Mechanisms associated with pyrethroid resistance in populations of Aedes aegypti (Diptera: Culicidae) from the Caribbean coast of Colombia

Aedes aegypti is the main vector of dengue, chikungunya, and Zika viruses, which are of great public health importance in Colombia. Aedes control strategies in Colombia rely heavily on the use of organophosphate and pyrethroid insecticides, providing constant selection pressure and the emergence of resistant populations. In recent years, insecticide use has increased due to the increased incidence of dengue and recent introductions of chikungunya and Zika. In the present study, pyrethroid resistance was studied across six populations of Ae. aegypti from the Caribbean coast of Colombia. Susceptibility to λ-cyhalothrin, deltamethrin, and permethrin was assessed, and resistance intensity was determined. Activity levels of enzymes associated with resistance were measured, and the frequencies of three kdr alleles (V1016I, F1534C, V410L) were calculated. Results showed variations in pyrethroid susceptibility across Ae. aegypti populations and altered enzyme activity levels were detected. The kdr alleles were detected in all populations, with high variations in frequencies: V1016I (frequency ranging from 0.15–0.70), F1534C (range 0.94–1.00), and V410L (range 0.05–0.72). In assays of phenotyped individuals, associations were observed between the presence of V1016I, F1534C, and V410L alleles and resistance to the evaluated pyrethroids, as well as between the VI1016/CC1534/VL410 tri-locus genotype and λ-cyhalothrin and permethrin resistance. The results of the present study contribute to the knowledge of the mechanisms underlying the resistance to key pyrethroids used to control Ae. aegypti along the Caribbean coast of Colombia.

The present study builds upon earlier work by further investigating the intensity and spatial extent of pyrethroid resistance in Ae. aegypti along the Caribbean coast of Colombia and links the frequency of kdr alleles and tri-locus kdr haplotypes to insecticide resistant phenotypes. To further understand the mechanisms of resistance, we also analyzed the activity levels of key detoxification enzyme groups.

Materials and methods
The insectary that was used to rear mosquitoes belonged to the Public Health Laboratory of the Department of Atlántico. Mosquitoes were reared using their existing protocols which include the routine use of feeding mosquito colonies on mice. This is standard practice in government insectaries in Colombia and is overseen by the Colombian National Institute of Health's National Public Health Laboratory Network.  (Fig 1).

Ae. aegypti collections
The study was undertaken in collaboration with and under the supervision of the Entomology Laboratory of the Colombian National Institute of Health, which coordinates the Entomology Laboratory Network in Colombia. This network for entomological surveillance includes the monitoring of insecticide resistance in insects of public health importance, to guide decision-making in the control of vector-borne diseases at the national level. The project also had the support of the local health secretaries in each of the departments selected for sampling sites.
The study sites were selected taking into account the national insecticide susceptibility baseline records. Municipalities with a susceptibility baseline were selected in order to identify variations in susceptibility over time and expand knowledge regarding enzymatic and molecular resistance mechanisms. Sites without an existing baseline were also selected in order to generate baseline data regarding the susceptibility to insecticides in these areas.
Immature Aedes were collected from habitats including tanks, pools, plastic/metallic cans, tires, animal water dishes, and flower vases located around houses. The specimens were reared to adults and maintained under controlled conditions of temperature (28˚C ± 2˚C), relative humidity (60% ± 10%), and photoperiod (12 h light:12 h dark) in the Public Health Laboratory of the department of Atlántico.
Upon emergence, male mosquitoes were fed with 10% sugar solution, and the females were fed with mouse (Mus musculus) blood every third day to obtain eggs of the F1 generation. Eggs were stored in an airtight plastic container, until they were hatched to obtain the adult mosquitoes used in the bioassays.

Bioassays
Insecticide bioassays were performed following the methodologies described by the CDC [33] and WHO [34]. The pyrethroid insecticides and their concentrations were as follows: λ-cyhalothrin ( For each population, 20-25 F1 generation, 3-to 5-day-old, unfed female Ae. aegypti were used in each bioassay replicate; as a control, the susceptible Rockefeller laboratory Ae. aegypti strain was used. Each bioassay consisted of four replicates per insecticide for each population. The diagnostic time post-exposure was 30 min for the CDC bioassays and 24h for the WHO bioassays. Upon the completion of the diagnostic time, the living and dead specimens were classified as phenotypically resistant (R) or susceptible (S), and individually stored in 0.5 mL tubes with a hole in the lid and desiccated in tightly sealed bags containing silica gel. The bags containing the tubes were stored at −80˚C for the subsequent detection of the V1016I, F1534C, and V410L kdr alleles.
In populations where resistance was detected via the CDC bioassay, resistance intensity was determined by conducting additional bioassays employing 2X the original insecticide concentration [34].

Biochemical assays
Biochemical assays were conducted on F1 generation adults. One day post-emergence, 40 unfed female Ae. aegypti from each population were preserved at −80˚C until processing. Individuals from the susceptible Rockefeller strain were used as controls. Mosquitoes were homogenized individually in 30 μL of distilled water for 5-10 seconds with an electric macerator and an additional 270 μL of distilled water was added for a final volume of 300 μL. Subsequently, samples were centrifuged at 12,000 rpm for 60 seconds and aliquoted in triplicate in 96 well microplates: 10 μL for α, β, pNPA-esterases; 15 μL for GST; 20 μL for mixed-function oxidases (MFO); and 25 μL for iAChE. For the tests of MFO and acetylcholinesterase, the samples were transferred without being centrifuged. Enzyme activity levels were determined according to the methodology described by Valle et al. [35], which measures the optical densities at predetermined wavelengths to estimate the activity levels of MFO, iAChE, esterases, and GSTs. Total protein concentration was also determined for each individual mosquito to correct for differences in body sizes [36]. Results were read using an ELISA plate reader (Multiskan TM -Thermo Fisher Scientific1).

Detection of kdr alleles
Real-time PCR was used to identify the V1016I, F1534C, and V410L kdr mutations. To estimate the allele frequencies in natural populations, 40-50 Ae. aegypti parental (F0) mosquitoes from each population were analyzed. To estimate associations between genotype and phenotype, all phenotypically resistant (R) and 30 randomly selected susceptible (S) individuals were analyzed per insecticide per population.
DNA was extracted from individual mosquitoes using the Quanta Biosciences Extracta TM Kit. Individual mosquitoes were placed in sterile 0.2 mL tubes and 25 μL extraction buffer was added to each tube, followed by an incubation at 95˚C for 30 min in a C1000 Bio-Rad CFX 96 Touch TM Real-Time System thermocycler. At the end of the incubation, 25 μL of stabilization buffer was added. DNA was quantified using a NanoDrop TM 2000/2000c spectrophotometer (ThermoFisher Scientific).
PCR reactions were performed in a Bio-Rad C1000 CFX96 Real-Time System thermocycler. Genotype was determined by analyzing the melting curves of the PCR products. The V1016I mutation was amplified following the methodology described by Saavedra-Rodríguez et al. [37], using a final reaction volume of 20 μL, containing 6 μL of ddH 2 O, 10 μL of iQ TM SYBR1 Green Supermix (Bio-Rad), 1 μL of each of the V1016f, I1016f, and I1016r primers (Table 1), and 1 μL of DNA template. The cycling conditions were as follows: an initial denaturation at 95˚C for 3 min followed by 40 cycles of: 95˚C for 10 s, 60˚C for 10 s, and 72˚C for 30 s; and a final extension at 95˚C for 10 s. The melting curves were determined by a denaturation gradient from 65˚C to 95˚C with an increase of 0.2˚C every 10 seconds.
The F1534C mutation was detected following the methodology described by Yanola et al. [38], using a final reaction volume of 20 μL comprised of 7.15 μL of ddH 2 O, 9 μL of iQ TM SYBR1 Green Supermix (Bio-Rad), 0.6 μL of each of the F1534f, C1534f, and F1534r primers (Table 1), 0.65 μL of the C1534f primers, and 2 μL of DNA template. The cycling conditions were as follows: an initial denaturation at 95˚C for 3 min followed by 37 cycles of: 95˚C for Table 1. Primer sequences used for detecting kdr alleles.

Mutation
Primer 5´-TTCTTCCTCGGCGGCCTCTT-3h ttps://doi.org/10.1371/journal.pone.0228695.t001 10 s, 57˚C for 30 s, and 72˚C for 30 s; and a final extension at 95˚C for 10 s. The melting curves were determined by a denaturation gradient from 65˚C to 95˚C with an increase of 0.5˚C every 5 s. The V410L mutation was detected following the methodology described by Haddi et al. [39], using a final reaction volume of 21 μL comprised of 9.6 μL of ddH 2 O, 10 μL of iQ TM SYBR1 Green Supermix (Bio-Rad), 0.1 μL of each of the L410f, V410f, and L410r primers (Table 1), 0.2 μL of the L410r primer, and 1.0 μL of DNA template. The cycling conditions were as follows: an initial denaturation at 95˚C for 3 min followed by 39 cycles of: 95˚C for 10 s, 60˚C for 10 s, and 72˚C for 30 s; and a final extension at 95˚C for 10 s. The melting curves were determined by a denaturation gradient from 65˚C to 95˚C with an increase of 0.2˚C every 10 s.
Each mosquito was analyzed in duplicate. For all assays for each mutation, three positive controls were included: a wild-type homozygote, a homozygote mutant, and a heterozygote. All assays also included a negative control consisting of master mix without DNA template.

Data analysis
Bioassays. Mortality was scored at the diagnostic time per insecticide per population. Populations were categorized according to the WHO criteria [34], whereby 98%-100% mortality indicates susceptibility, 90%-97% suggests resistance is developing and <90% mortality indicates resistance.
Biochemical assays. Absorbance values were entered into Excel databases to calculate the mean and standard deviation for each mosquito. To express the absorbance values in terms of enzymatic activity, data regarding the homogenate volume of each mosquito, total protein content for each mosquito, and units of activity for each enzyme group were calculated according to the protocol described by Valle et al. [35]. The cutoff value for the susceptible Rockefeller strain was determined based on the 99 th percentile of absorbance, and the percentage of individuals from the field strains with activity levels that exceeded this cutoff value were classified according to the criteria proposed by Montella et al. [40]: <15% unaltered, 15%-50% altered, and >50% highly altered.
After determining the activity levels for each enzyme group, an analysis of variance was performed, followed by Tukey's multiple comparison test, with the significance level set at p � 0.05, to identify populations with any statistically significant differences as compared to the Rockefeller reference strain.
Allelic and genotypic frequencies of the V1016I, F1534C, and V410L mutations. Results were obtained using Bio-Rad's Precision Melt Analysis Software TM and were interpreted as follows. For the V1016I mutation, a melting peak at 77˚C corresponded to a mutant homozygote (I/I), a peak at 82˚C corresponded to a wild-type homozygote (V/V), and peaks at both 77˚C and 82˚C corresponded to a heterozygote (V/I). For the F1534C mutation, a peak at 82˚C corresponded to a mutant homozygote (C/C), a peak at 78˚C corresponded to a wildtype homozygote (F/F), and peaks at both 78˚C and 82˚C corresponded to a heterozygote (F/ C). For the V410L mutation, a peak at 80˚C corresponded to a mutant homozygote (L/L), a peak at 83˚C corresponded to a wild-type homozygote (V/V), and peaks at both 80˚C and 83˚C corresponded to a heterozygote (V/L).
n mosquitoes with the genotype to be calculated total n mosquitoes analyzed ð2Þ The Hardy-Weinberg principle was tested, as shown in Eq (3) where p is the number of wild-type homozygotes, pq is the frequency of heterozygotes, and q is the frequency of mutant homozygotes. Expected wild-type V 1016 /V 1016 , F 1534 /F 1534 , V 410 /V 410 homozygotes = p 2 (n). Expected V 1016 /I 1016 , F 1534 /C 1534 , V 410 /L 410 heterozygotes = 2pq (n). Expected mutant I 1016 /I 1016 , C 1534 /C 1534 , L 410 /L 410 homozygotes = q 2 (n). The Chi square test was used to determine whether the populations were in Hardy-Weinberg equilibrium, as shown in Eq (4): f 0: Frecuency observed value . f e: Frecuency expected value . If the calculated value of χ 2 was < tabulated χ 2 (1 gl) = 3.84 and P < 0.05, the H 0 that the study population was in Hardy-Weinberg equilibrium was accepted; otherwise, if the calculated χ 2 was � tabulated χ 2 , the H a that the study population was not in Hardy-Weinberg equilibrium was accepted.
In addition, the coefficient of endogamy was calculated using Eq (5) as follows: where, H obs is the number of observed heterozygotes and H exp is the number of expected heterozygotes; if F IS was significantly higher than 0, an excess of homozygotes was considered, and if F IS was significantly less than 0, an excess of heterozygotes was considered in the population, with a significance of P < 0.05. In addition, the frequencies of tri-locus genotypes were determined in the study populations. Association of kdr mutations with pyrethroid resistance. The association between resistant and susceptible phenotypes and their kdr genotypes was tested using contingency tables, and the relationship between phenotype and tri-locus genotype was tested using the statistical software programs OpenEpi version 3.0 (https://www.openepi.com/TwobyTwo/TwobyTwo. htm) and GraphPad Prism version 8.1.

Bioassays
A total of 1732 adult female Ae. aegypti were tested in WHO bioassays for susceptibility to λcyhalothrin (n = 564), deltamethrin (n = 586), and permethrin (n = 582). Resistance to λ-cyhalothrin and permethrin was detected in all six evaluated populations. Resistance was most frequent in Montería with 43.3% mortality to λ-cyhalothrin and 24.0% mortality to permethrin. Cartagena was the least resistant, with mortalities of 86.4% to λ-cyhalothrin and 77.6% to permethrin. Susceptibility to deltamethrin was observed in the populations from Juan de Acosta (98% mortality) and Barranquilla (100% mortality), and possible development of resistance was detected in Valledupar (96.8% mortality) and Montería (93.2% mortality). The populations from Cartagena (87.9% mortality) and Chiriguaná (86.0% mortality) were found to be resistant to deltamethrin (Fig 2).
In populations where resistance intensity was assessed, 100% mortality at the diagnostic time was observed after exposure to twice the concentration (2X) of the recommended diagnostic dose for λ-cyhalothrin (20 μg/bottle) and permethrin (30 μg/bottle) ( Table 2).

PLOS ONE
(44%), Montería (34%) and Chiriguaná (4%) (Fig 4E). AChE activity remained unaltered in all populations evaluated ( Fig 4F). Overall, significant differences were observed between the mean activity levels of most enzyme groups between the field populations and the Rockefeller reference strain (p < 0.05) (Fig 4).  frequency of 0.05. For the other populations, the frequencies of the L410 allele ranged between 0.12 and 0.32 (Table 3). For loci 1016 and 1534, all genotypes were found to be in Hardy-Weinberg equilibrium. In the case of locus 410, the genotypes of most populations, except Valledupar, were in Hardy-Weinberg equilibrium (p < 0.05). When determining the inbreeding coefficients (F IS ) for I1016, values < 0 were obtained for the populations from Barranquilla and Valledupar due to an excess of heterozygotes, in contrast to the populations from Cartagena, Chiriguaná, Juan de Acosta, and Montería, where values > 0 were recorded due to a deficiency of heterozygotes. For C1534, a generalized excess of heterozygotes was observed, with the exception of Barranquilla, where a deficiency of heterozygotes was observed. Similarly, for L410, the populations from Barranquilla, Cartagena, and Valledupar showed a deficiency of heterozygotes, in contrast to Chiriguaná, Juan de Acosta, and Montería, where an excess of heterozygotes was detected (Table 3).

kdr allele frequencies
Of the 27 combinations of tri-locus genotypes, 13 combinations were detected in 281 mosquitoes collected from the six evaluated populations. The triple homozygous wild-type genotype (VV 1016 / FF 1534 /VV 410 ) was detected only in the populations from Barranquilla and Juan de Acosta, with frequencies of 0.08 and 0.04, respectively, whereas the triple homozygous mutant genotype (II 1016 /CC 1534 /LL 410 ) was present in all populations except Valledupar, with frequencies between 0.02 (Barranquilla) and 0.49 (Montería). Similarly, the triple heterozygous  5).

Association of kdr alleles with phenotypic resistance to pyrethroids
Based on the results obtained with the mosquitoes exposed to insecticides in the WHO bioassays, a significant association (p < 0.05) was identified between the mutant kdr alleles 1016I, 1534C, and 410L and resistance to λ-cyhalothrin in the populations from Juan de Acosta, Montería, and Valledupar. Similarly, an association was observed between the 1534C allele and resistance to deltamethrin in the populations of Chiriguaná, Montería, and Valledupar and between the 1016I and 410L alleles and resistance to deltamethrin in the population of Montería. A significant association (p < 0.05) was also detected between the 1016I, 1534C, and 410L alleles and resistance to permethrin in the populations from Chiriguaná, Montería, and Valledupar; between the 1534C allele and resistance to permethrin in Barranquilla, Cartagena, and Juan de Acosta; and between the 410L allele and permethrin resistance in Juan de Acosta (Tables 4-6).  Less association was detected between kdr alleles and the observed phenotype in the CDC bioassays. A significant association (p < 0.05) between the 1534C allele and resistance to λcyhalothrin was detected in the population from Barranquilla and between the 1016I and 410L alleles and resistance to permethrin in the population from Montería. Despite the resistance to

PLOS ONE
pyrethroids detected with the CDC bioassays in the populations from Chiriguaná and Juan de Acosta, no significant associations were detected between kdr alleles and resistant phenotypes in these populations (Tables 7 and 8).  susceptible to deltamethrin (p < 0.05). Heterozygotes at both loci 1016 and 410 in the presence of CC1534 were significantly more likely to be resistant to λ-cyhalothrin and permethrin (p < 0.05) and susceptible to deltamethrin (p < 0.05) ( Table 9). From the CDC bioassays, 15 combinations of tri-locus genotypes were observed in 465 mosquitoes assayed with λ-cyhalothrin and permethrin in Barranquilla, Juan de Acosta, and Montería. Similar to the WHO bioassays, the most common haplotypes were VI 1016 /CC 1534 / VL 410 (n = 161, 34.6%) and VV 1016 /CC 1534 /VV 410 (n = 117, 25.2%). Wild-type double homozygotes at loci 1016 and 410 in the presence of CC1534/FC1534 were significantly more likely to be phenotypically susceptible to λ-cyhalothrin and permethrin (p < 0.05) ( Table 10).

Discussion
In Colombia, the use of pyrethroids for the control of Ae. aegypti is a fairly recent phenomenon. Among the pyrethroids, λ-cyhalothrin and deltamethrin have most commonly been used to control Ae. aegypti in Colombia. However, resistance to λ-cyhalothrin has been more commonly reported than resistance to deltamethrin in Colombia, as demonstrated by results from previous studies [7, 9-11, 13, 26] as well as those obtained in the present study. In the findings presented here, we detected resistance to permethrin and λ-cyhalothrin in all populations and varying degrees of susceptibility to deltamethrin. This heterogeneity of resistance patterns within the pyrethroid class suggests that diverse mechanisms are contributing to these phenotypes.
Resistance to DDT is widespread in Colombia owing to the application of this organochlorine compound for more than five decades in the country [20]. DDT and pyrethroids share the mode of action consisting of delayed sodium channel closure and membrane repolarization [41]. The modification of this target site due to the presence of kdr mutations on the para gene can lead to cross-resistance to both DDT and pyrethroids. As such, the high prevalence of kdr alleles detected in our study may also be linked to previous selection pressures caused by DDT [42,43]. Our findings of resistance to λ-cyhalothrin in populations of Ae. aegypti in the Caribbean region of Colombia are consistent with those reported by Maestre et al. [13] for the populations of Barranquilla and Montería, Granada et al. [15] for the population of Riohacha, and Atencia et al. [31] for the population of Sincelejo. Likewise, Maestre et al. [13] reported resistance to λ-cyhalothrin in Valledupar and moderate resistance in Cartagena, while in the present work we detected resistance using the WHO bioassay but susceptibility using the CDC bioassay in both populations. However, other studies in Colombia have found resistance to λcyhalothrin using both the WHO and CDC techniques in the departments of Cundinamarca, Caquetá, Meta, Guaviare, Santander, Chocó, Antioquia, Putumayo and Casanare [9,10,12]. Ocampo et al. [11] found susceptibility to this insecticide using the CDC technique in the departments of Cauca, Nariño, Valle del Cauca and Huila.    Our findings of resistance to permethrin were consistent with those reported by Maestre et al. [13] for Barranquilla and Montería; however, for Cartagena and Valledupar, Maestre et al. had previously reported susceptibility, while we observed resistance using the WHO bioassay but susceptibility using the CDC bioassay. For this same insecticide Ardila et al. [12] encountered resistance in the department of Casanare, as did Fonseca et al. [10] in the departments of Chocó Antioquia and Putumayo using both bioassay methodologies. The results of the intensity bioassays for both permethrin and λ-cyhalothrin for the populations that had shown resistance using the CDC methodology resulted in 100% mortality when the diagnostic dose was doubled, suggesting that while resistance was present in the populations, it had not yet reached a high level of intensity.
For deltamethrin, Maestre et al. [13] had reported resistance in Barranquilla; however, we found susceptibility using both bioassay methodologies. In Montería and Valledupar, Maestre et al. [13] reported resistance to deltamethrin in both populations, while we found susceptibility with the CDC bioassay and indications that resistance was developing with the WHO bioassay. In Cartagena, Maestre et al. [13] had reported moderate resistance to deltamethrin, while we observed resistance with the WHO bioassay and susceptibility with the CDC bioassay. Other studies carried out in Colombia showed resistance in Ae. aegypti to deltamethrin in the departments of Cundinamarca and Santander and susceptibility in the departments of Caquetá and Meta [9], as well as susceptibility in the departments of Cauca, Nariño, Valle del Cauca, Huila [11], Casanare [12] and Caldas [14] using the two bioassay methodologies.
In Colombia, both the CDC bioassay methodology and the WHO bioassay methodology have been used for insecticide susceptibility studies in adult Ae. aegypti. Typically, using both techniques, resistance to DDT has been observed in all Ae. aegypti populations evaluated in the country, together with variable susceptibility to pyrethroids and susceptibility to organophosphates in most populations [7]. In the present study, some discrepancies were observed between the results obtained with the WHO and CDC bioassay methodologies, indicating that the two techniques may not always provide consistent results. In studies by Aizoun et al. [44] and Fonseca et al. [20], WHO and CDC bioassays were compared to determine the susceptibility of Anopheles gambiae to deltamethrin and Anopheles nuñeztovari to fenitrothion. Both studies reported susceptibility when using the WHO bioassay and resistance when using the CDC bioassay. The authors observed that the exposure time of the mosquitoes to the insecticide (diagnostic time) was considerably shorter in the case of the CDC bioassay, which could have led to an overestimation of resistance; although in fact the opposite was observed in our study. Despite the shorter exposure time in the CDC bioassay, populations that were classified as resistant in the WHO bioassay were classified as susceptible in the CDC bioassay. This could potentially be explained due to the mechanisms underlying the resistance; for example, resistance that is primarily caused by kdr would likely result in populations that are not quickly knocked down and thus scored as 'resistant' at 30 minutes. However, if the main mechanisms of resistance are metabolic, mosquitoes may initially be knocked down but could recover over time as their detoxification enzymes metabolize the insecticide. Indeed, our biochemical assay data suggest that elevated enzymatic activity is present in the populations that were studied.
Most previous studies regarding enzymatic activity have been conducted on Ae. aegypti populations from other regions of Colombia where alterations were detected, mainly in MFOs and nonspecific esterases, in populations from Antioquia, Chocó, Putumayo, Cauca, Valle del Cauca, Nariño, Huila, Santander, Meta, and Casanare [9][10][11][12]. The one previous study conducted in the Caribbean region of Colombia reported altered α-esterases and MFOs in Ae. aegypti from Valledupar, MFOs in Cienaga, and GSTs in Sincelejo. In that previous study, no alterations in enzyme activity were detected in Cartagena, Montería, Barranquilla, San Juan, Puerto Colombia, and Soledad, [13]. Our results are consistent with the finding of highly altered MFOs in Valledupar, and we also detected altered β-esterases in that same population. We also detected highly altered α-esterases, β-esterases, MFOs and GSTs in Montería; altered β-esterases and GSTs in Barranquilla; and altered GSTs in Cartagena. Additionally, in the present study we detected altered pNPA-esterases in the population of Juan de Acosta.
There are no studies in Colombia incriminating insensitive acetylcholinesterase as a mechanism associated with resistance to organophosphates and carbamates in Ae. aegypti. A study by Grisales et al. [24] reported resistance to temephos in the population of Ae. aegypti from Cúcuta (RR: 15X) without evidence of insensitive acetylcholinesterase, although they did detect esterase and oxidase-based mechanisms.
Kdr mutations are important mechanisms involved in DDT and pyrethroid resistance. In Colombia, the first kdr mutation reported in populations of Ae. aegypti was V1016I, which was identified in populations from Puerto Colombia, Soledad, Barranquilla, Valledupar, San Juan, Sincelejo, Montería, Cienaga and Cartagena, which are all located in the Caribbean region. In that initial report, the V1016I mutation showed frequencies ranging between 0.07 and 0.35; the lowest frequency was found in the Cienaga population and the highest was found in Soledad, Montería, and Barranquilla, with frequencies of 0.35, 0.33, and 0.32, respectively [13]. The highest frequency of 1016I that we detected in the present study was in Montería, with a frequency of 0.70, showing a large increase in the frequency in this population from what was originally reported by Maestre et al. [13]. In addition, an increase in the frequency of 1016I from 0.09 to 0.16 was detected in Cartagena and a reduced frequency was detected in Barranquilla and Valledupar, from 0.32 and 0.27, respectively, to 0.15 in both populations. V1016I had also previously been reported in Quindío at low levels of frequency (0.02-0.05) [25].
The F1534C mutation was first detected in Colombia in the department of Sincelejo (Sucre), in the Caribbean region [31]. It had also previously been reported in Ae. aegypti populations from Puerto Colombia, Soledad, Barranquilla, Valledupar, San Juan, Sincelejo, Montería, Cienaga and Cartagena with frequencies ranging between 0.74 and 0.88. When compared with the results reported previously, we observed increased frequencies of 1534C, having risen in Barranquilla from 0.74 to 0.76, in Cartagena from 0.86 to 0.97, in Montería from 0.88 to 1.00, and in Valledupar from 0.82 to 0.94. These increases are likely attributable to the constant pressure exerted by pyrethroid insecticides, which were heavily applied during the period between the two studies for the control of dengue, chikungunya, and Zika. Although there are no previous studies reporting this mutation in Juan de Acosta and Chiriguaná, these populations also showed high frequencies (0.76 and 0.95, respectively). Moreover, high frequencies of 1534C have been reported in other areas of Colombia, including Villavicencio, Riohacha, and Bello, with frequencies of 0.63, 0.71, and 0.56, respectively [15]. In these latter three populations, the V410L mutation was also identified in Colombia for the first time, with frequencies of 0.46, 0.30, and 0.06, respectively. It is noteworthy that in that study, Ae. aegypti from Bello were susceptible to λ-cyhalothrin, whereas those from Riohacha and Villavicencio were resistant. In these latter two populations, the researchers detected a positive association between V410L and V1016I and resistance to λ-cyhalothrin. In the present study, the V410L mutation was detected for the first time in the study populations, with frequencies ranging between 0.05 in Valledupar and 0.72 in Montería. The frequencies of the V1016I mutation were very similar to those of the V410L mutation in all the evaluated populations; this result is consistent with the findings reported by Granada et al. [15] for Ae. aegypti in Bello, Villavicencio, and Riohacha.
Haddi et al. [39] reported the presence of the V410L mutation in resistant Ae. aegypti in Brazil and observed that this mutation, either alone or in combination with the F1534C mutation, was strongly associated with increased the resistance to type I and II pyrethroids. This is consistent with the results of the present study, where the 1534C and 410L alleles were associated with resistance to permethrin in the population of Juan de Acosta. The 1016I, 1534C, and 410L alleles were all associated with resistance to permethrin in the Chiriguaná, Montería, and Valledupar populations based on phenotyping by the WHO bioassay. In addition, F1534C was associated with resistance to deltamethrin in Chiriguaná, Valledupar, and Montería; V1016I and V410L were also associated with deltamethrin resistance in the case of the latter population. Similarly, an association was found between all three mutations and resistance to λ-cyhalothrin in Valledupar, Montería, and Juan de Acosta. This last result is consistent with the results of the study by Maestre et al. [32] which detected a significant positive correlation between the frequency of the 1016I allele and resistance to permethrin, λ-cyhalothrin, and cyfluthrin. However, no significant correlation was observed in that same study between 1534C and resistance to any pyrethroids [32].
Recent studies conducted in Mexico proposed three sequential models to explain the evolution of the V1016I, F1534C, and V410L mutations. The first model suggests that F1534C appeared first, providing low resistance levels, followed by the appearance of V1016I, which provided higher levels of resistance. The second model challenges the first model and proposes that V410L and V1016I occurred independently on a C1534 haplotype followed by cis conversion by recombination. Finally, a third model assumes that the three mutations appeared independently at low frequencies and that two recombination events rearranged them in a cis configuration [52]. Considering these previous models and the results obtained in the present investigation, it is possible to hypothesize that the appearance of V410L and V1016I did not occur independently because their allelic frequencies were so similar and they almost always appeared together.
Regarding the 1016/1534/410 phenotype-genotype association, a relationship between the VI 1016 /CC 1534 /VL 410 genotype and resistance to λ-cyhalothrin and permethrin was detected in the present study. These results are consistent with the study conducted by Haddi et al. [39] in a pyrethroid-resistant Ae. aegypti strain from Brazil, where V410L alone or in combination with F1534C was shown to reduce sodium channel sensitivity to type I (permethrin) and type II pyrethroids (λ-cyhalothrin and deltamethrin). In addition, these results further support the notion that the presence of VI 1016 and VL 410 heterozygotes is sufficient to confer resistance to deltamethrin [52]. These findings suggest that the interactions of multiple mutations play a role in the response of Ae. aegypti sodium channels to insecticides [53].
In Colombia, previous studies have identified the frequency of V1016I, F1534C and V410L in populations of Ae. aegypti and have correlated those frequencies with the results of outcomes obtained through the CDC bioassay technique for populations from Puerto Colombia, Soledad and Barranquilla (Atlántico); Valledupar and San Juan (Cesar); Sincelejo (Sucre); Montería (Córdoba); Cienaga (Magdalena); Cartagena (Bolívar) [32]; Villavicencio (Meta); Riohacha (La Guajira); Bello (Antioquia) [15]; Giron (Santander); Buga, Palmira, Yumbo, Cali (Valle del Cauca) and Medellin (Antioquia) [54]. Unlike the aforementioned studies, our study evaluated for the first time for Colombia the direct association between phenotype and genotype in individual mosquitoes phenotyped through CDC and WHO bioassays. Furthermore, the present study established for the first time for Colombia differences between the phenotypes observed in CDC and WHO bioassays for pyrethroids and trilocus kdr haplotypes present in populations of Ae. aegypti from Montería (Cordoba); Cartagena (Bolívar); Juan de Acosta, Barranquilla (Atlántico); Chiriguaná and Valledupar (Cesar). This contributes to knowledge about the role and co-occurrence of these mutations in Ae. aegypti from the Caribbean region and how they relate to phenotypic resistance, providing further evidence to guide the selection of insecticides to be used in Ae. aegypti control.
A key strength of the present study is that it expands the knowledge base regarding the susceptibility status of Ae. aegypti to pyrethroid insecticides in Colombia. Additionally, it provides new information regarding the frequency and distribution of kdr mutations and a detailed analysis of phenotype-genotype associations. An increased understanding of the role of the mechanisms involved in resistance will contribute to improved resistance surveillance strategies which can better guide control programs for the selection of insecticides for the control of Ae. aegypti. An important limitation of the current study is the lack of synergist bioassays, which would permit an estimation of the relative contribution of metabolic mechanisms as compared to kdr in conferring phenotypic resistance.
Regardless, the results presented here provide important input for territorial and national entities in vector control decision-making. In addition, the detailed information on resistance mechanisms provide deeper insight into the types of resistance mitigation and management strategies might be most effective in the populations evaluated.

Conclusions
Variability was observed in pyrethroid susceptibility using the WHO and CDC bioassay methodologies, highlighting the importance of using a consistent methodology to routinely screen populations for susceptibility. The altered activity levels of β-esterases, α-esterases, MFOs, and GSTs suggest that metabolic resistance may be important in these populations. The kdr mutations V1016I, F1534C, and V410L were detected in all populations, with 1534C being nearly fixed in all except two populations. Finally, associations were observed between the F1534C mutation and resistance to permethrin in all populations, the F1534C mutation with resistance to deltamethrin in Chiriguaná, Montería, and Valledupar, and the V1016I, F1534C, and V410L mutations and resistance to λ-cyhalothrin in Juan de Acosta, Valledupar, and Montería.