Impact of deltamethrin selection on kdr mutations and insecticide detoxifying enzymes in Aedes aegypti from Mexico

Background Insecticide resistance is a serious problem for vector control programmes worldwide. Resistance is commonly attributed to mutations at the insecticide’s target site or increased activity of detoxification enzymes. Methods We determined the knockdown concentration (KC50) and lethal concentration (LC50) of deltamethrin in six natural populations of adult Aedes aegypti from southeastern Mexico. These populations were then selected over five generations using the LC50 from the preceding generation that underwent selection, and the heritability of deltamethrin resistance was quantified. For each generation, we also determined the frequency of the kdr alleles L410, I1016 and C1534, and the levels of activity of three enzyme families (α- and β-esterases, mixed-function oxidases and glutathione S-transferases (GST)) associated with insecticide detoxification. Results There was an increase in KC50 and LC50 values in the subsequent generations of selection with deltamethrin (FS5vs FS0). According to the resistance ratios (RRs), we detected increases in LC50 ranging from 1.5 to 5.6 times the values of the parental generation and in KC50 ranging from 1.3–3.8 times the values of the parental generation. Triple homozygous mutant individuals (tri-locus, LL/II/CC) were present in the parental generations and increased in frequency after selection. The frequency of L410 increased from 1.18-fold to 2.63-fold after selection with deltamethrin (FS5vs FS0) in the populations analyzed; for I1016 an increase between 1.19-fold to 2.79-fold was observed, and C1534 was fixed in all populations after deltamethrin selection. Enzymatic activity varied significantly over the generations of selection. However, only α- esterase activity remained elevated in multiple populations after five generations of deltamethrin selection. We observed an increase in the mean activity levels of GSTs in two of the six populations analyzed. Conclusions The high levels of resistance and their association with high frequencies of kdr mutations (V410L, V1016I and F1534C) obtained through artificial selection, suggest an important role of these mutations in conferring resistance to deltamethrin. We highlight the need to implement strategies that involve the monitoring of kdr frequencies in insecticide resistance monitoring and management programmes.


Background
Aedes aegypti (L.) is an arbovirus vector of great public health importance, and is the principal vector transmitting dengue, chikungunya and Zika viruses in the Americas [1][2][3]. Currently, there are no widely licensed vaccines or specific treatments for any of these diseases, rendering vector control the principal strategy for preventing their transmission. The reduction or elimination of larval habitats and chemical control with the use of larvicides and adulticides are the main strategies to reduce vector-human contact and break the transmission cycles of these diseases [4]. The insecticides used for such purposes worldwide are mainly from the following chemical groups: organophosphates, organochlorines, carbamates and pyrethroids [5,6]. Pyrethroids are often the first line of defense to break the transmission of these diseases, and they are the most commonly used worldwide due to their low toxicity to mammals including humans [7,8].
In Mexico, the application of pyrethroid insecticides is regulated by the Ministry of Health, and their use dates back to 1999 when the use of DDT was banned in the health sector. At this time, permethrin became the insecticide of choice for space spraying and deltamethrin for indoor residual spraying [9]. In 2011, the use of permethrin was discontinued; however, deltamethrin remained on the list of insecticides approved for use as an adulticide in indoor residual spraying. In 2012 other pyrethroids were added to the list; however, deltamethrin remained, and in 2016 its use for space spraying applications was also approved, a practice that remains to date [10]. The use of pyrethroids in Mexico has been favored for more than 20 years, which led to widespread resistance in mosquito populations [11][12][13].
Insecticide resistance is one of the main obstacles to the success of vector control programmes [12,14,15]. Two of the principal mechanisms that cause resistance in mosquitoes are alterations at the target site of the insecticide, including knockdown resistance (kdr) mutations on the voltage-gated sodium channel gene (vgsc), and increased metabolic activity [16,17]. Targetsite resistance arising from kdr mutations can cause resistance to both DDT and pyrethroids, since both types of insecticides target sodium channels at the axon level [18].
Since the discovery of the L1014F mutation in Musca domestica and since the determination of its role in insecticide resistance [19], more than 50 mutations have been described for the vgsc gene in different pests and disease vectors [20]. Thirteen of these sequence changes are commonly found in Ae. aegypti in ten different positions covering five regions of the vgsc gene (IIS4-5, IIS5-6, IIS-6, IIIS6 and IVS5), identified at positions 410 (V → L), 419 (V → L), 923 (G → I), 982 (L → W), 989 (S → P), 1011 (I → M or V), 1016 (V → G or I), 1520 (T → I), 1534 (F → C or L), and 1763 (D → Y) [21][22][23][24][25][26][27][28][29]. It should be noted that only five of the mutations described for this species have been functionally associated with sensitivity to pyrethroids and DDT by heterologous expression and electrophysiology assays; only S989P, I1011M, V1016G, F1534C and V410L have been confirmed [20,27,30,31]. Similarly, it is important to note that mutations are usually associated with a specific geographic sector, for example, mutation V1016G has been described in Southeast Asia but not in the Americas [32].
In Mexico, mutations F1534C, V1016I and V410L have been reported in populations of Ae. aegypti in different states, sites where their temporal distribution indicates the absence of mutant alleles in 2000, with a gradual increase in these mutations during the period of 2002-2008, up to fixation of the mutant alleles in 2016, showing an allelic frequency for the three simultaneous changes of 0.47 [30]. Of these mutations, the change from phenylalanine to cysteine at position 1534 was described initially by [25], who associated its occurrence with pyrethroid resistance; however, confirmation by heterologous expression and electrophysiology has demonstrated that this change reduces sensitivity of VGSC to type I but not type II pyrethroids [33]. The presence of this mutation has been linked to the development of low resistance levels and its association with other gene changes to the establishment of a more resistant phenotype. Such is the case of its association with V1016I, described by Saavedra et al. [22], where its participation in the reduction of sensitivity of VGSC has not been identified when appearing as a single mutation [31]. An analysis of the frequency of V1016I and F1534C in 24 populations of Ae. aegypti in Mexico showed that the most likely evolutionary process is the appearance of F1534C and then V1016I, this being corroborated by the low fitness of haplotype I1016/F1534 [34] and the increase in insecticide resistance expressed in Xenopus oocytes when inducing the two mutations [31]. Temporal distribution studies of this mutation in Mexico identified an increase in the frequency of the I1016 mutant allele during the period of 1996 to 2000 with a frequency of 0.04% up to an increase to 33.2% in 2007 to 2009 [35].
The V410L mutation has been recently described by Haddi et al. [27], who determined their association with resistance to type I and II pyrethroids. This mutation has been studied in Mexico in populations of Ae. aegypti collected between 2000 and 2016 where its appearance was reported in 2002 in heterozygous individuals collected in Coatzacoalcos, Veracruz, and where the increase in the frequency of L410 was also documented as being as high as 0.9 in populations collected in 2014 in Merida, Yucatan [30]. A similar pattern was found for 26 populations of Ae. aegypti collected in the eastern part of Mexico, showing a frequency of the L410 allele close to the one in populations from Minatitlan and Jose Cardel in Veracruz and Cancun in Quintana Roo (0.99, 0.97 and 0.93, respectively), but the average values, as well as the interval of frequencies observed in populations ranged from 0.3 to 0.99 [36].
Metabolic resistance to pyrethroids is mediated mainly by glutathione S-transferases (GSTs), esterases and mixed-function oxidases (MFOs) [37][38][39]. In Mexico, Flores et al. [14] previously reported high levels of αand β-esterase activity related to permethrin selection in populations of Ae. aegypti from Baja California Norte and Sur. The increased activity of these enzymes has also been reported in populations from the state of Sonora in a study of permethrin resistance and in the state of Veracruz in a study of chlorpyrifos resistance [40,41]. A further study reported high levels of GST activity in the state of Guerrero, which was found to be related to resistance to DDT [12].
Given the background levels of pyrethroid and DDT resistance in many parts of Mexico, the extent to which ongoing insecticide pressure selects for specific mechanisms of resistance is not well understood. To better understand how insecticide exposure selects for resistance we evaluate the effect of laboratory selection with deltamethrin on the underlying resistance mechanisms in populations of Ae. aegypti from southeastern Mexico, over five generations of selection.

Biological material
Eggs and larvae of Ae. aegypti were collected during 2014 from six sites: Merida, Progreso, Hunucma, and Hoctun in the State of Yucatan; and Agua Dulce and Jose Cardel, in the State of Veracruz (Table 1, Fig. 1). The criteria for site collection were: high incidence of dengue and frequent application of insecticides for its control; high infestation index of homes with Ae. aegypti; and previous registration of insecticide resistance. The collection of the entomological material was carried out in peridomestic sites. In each, entomological inspections were performed for the collection of stages of Ae. aegypti in at least 20 potential breeding places per site, such as containers, tires and flowerpots. The biological material obtained from the field was reared under insectary conditions at 26 ± 2 °C and 70-80% relative humidity with a 12h:12h (light:dark) photoperiod. The adults obtained from field-collected material were allowed to intermate, and the eggs resulting were designated F S0 (without previous selection). The New Orleans strain (NO) was used as a susceptible reference in the study, this strain was originally obtained from the CDC (Atlanta, GA, USA) and has been maintained since 2002.

Bioassays
The bioassays consisted of exposing 20-25 non-bloodfed, 2-3-day-old F S0 female mosquitoes to different concentrations of the pyrethroid deltamethrin (> 98% purity; Chem Service, West Chester, PA, USA) based on the bottle bioassay methodology described by the CDC [42]. Each bioassay consisted of bottles containing different concentrations of deltamethrin resulting in mortalities between 9-90%, with at least three replicates for each concentration, and an untreated control bottle. The numbers of knocked-down mosquitoes were recorded at 1 h. After 1 h of exposure, all mosquitoes were gently transferred to a recovery container without insecticide and were offered a cotton ball soaked in a sugar solution. Mortality was recorded at 24 h. Both the bottles and recovery containers were kept at 24 ± 2 °C and 70% RH.
Rates of kdr and 24 h recovery were analyzed using a logistic regression model, QCal (https ://sourc eforg e.net/ proje cts/irmap roj/files /Qcal/) [43]. We determined the KC 50 (concentration causing 50% knockdown) after 1 h of exposure. The LC 50 (concentration causing 50% mortality) was estimated from mortality data 24 h after recovery. The confidence intervals were calculated using a significance level of α = 0.05. The mortalities were corrected according to Abbott's formula [44] when mortality was observed in the control bottles. The resistance ratios (RR) were calculated by dividing the KC 50 or LC 50 by the KC 50 or LC 50 of the susceptible New Orleans (NO) Ae. aegypti reference strain. The magnitude of resistance was classified as high (RR > 10-fold), moderate (RR between 5-10-fold) or low (RR < 5-fold) according to the criteria proposed by Mazarri & Georghiou [45].

Selection with deltamethrin
Selection cohorts were generated according to the methodology described by Saavedra et al. [46]. The selection was carried out exposing 650-1000 F S1 males and females (1:4) for 1 h to the LC 50 obtained for the previous generation (F S0 ) by the method described above. The survivors at 24 h were transferred to mosquito cages for breeding the subsequent generation (F S2 ). The dose response parameters were calculated again for the F S1 cohort and were used to select the next generation (F S2 ). This procedure was repeated for each generation of selection until F S5 . The New Orleans strain was tested at the same time as the F S0 and was not tested in each generation of selection.
KC 50 and LC 50 values were compared between generations to monitor changes in resistance. The KC 50 and LC 50 values were considered significantly different if their 95% confidence intervals did not overlap.

Realized heritability
Resistance risk assessment was made by calculating realized heritability values (h 2 ) in all selection cohorts for the knockdown (KC 50 ) and lethal (LC 50 ) parameters as described by Tabashnik [47]. The heritability index (h 2 = R/S) is calculated as the ratio of the response to selection (R) to the selection differential (S) according to the artificial selection technique of Falconer [48] and is related to the additive genetic variance for a trait. A low h 2 predicts no additive genetic variance for a trait and a poor or very slow response to artificial selection, while a high h 2 predicts a large additive genetic variance at one or a few loci that govern a trait and predicts a rapid response to artificial selection.
The selection differential is expressed as the product of selection intensity (i) and phenotypic standard deviation (S= iσ). The response to selection was determined for each generation as the difference in population means between subsequent generations with the probit analysis: R = log (Final KC 50 or LC 50 ) − log (Initial KC 50 or LC 50 )/n, where final KC 50 or LC 50 area the values of the offspring after n generations of selection and initial KC 50 or LC 50 are the values of the parental generation before n generations of selection. The selection intensity was calculated from the proportion of surviving individuals in the examined population. The difference between KC 50 or LC 50 was calculated on a logarithmic scale because the logarithm of tolerance was assumed to be normally distributed. The phenotypic standard deviation at t-th generation (σ t ) was obtained as the inverse of the regression slope: σ t = 1/b t . The parameters (R, i and σ) were determined at every generation and the realized heritability was estimated as the regression coefficient of cumulative responses on cumulative selection differentials [48].

Molecular assays
DNA was extracted from individual mosquitoes by the technique described by Coen et al. [49], and the DNA pellet was resuspended in 30 μl of ultrapure molecular grade water (Corning Cellgro TM , Manassas, VA, USA). The quantity and quality of DNA were determined using a NanoDrop spectrophotometer 2000 (Thermo Fisher Scientific, Woonsocket, RI, USA).
The kdr alleles 1534C, 1016I and 410L were detected by endpoint PCR using a Bio-Rad (Hercules, CA, USA) T100TM thermocycler using the primers for each of the mutations described by Saavedra et al. [22,30] and Yanola et al. [25].
DNA from the New Orleans susceptible strain was used as a negative control for the kdr assays, and previously genotyped individuals were used as positive controls.
The frequencies of the alleles were determined for each population at each generation. We verified that the populations at the parental generation (F S0 ) were in Hardy-Weinberg equilibrium by means of Chi-square analysis. In addition, Wright's F IS inbreeding coefficient was estimated, along with Wald's correction [51,52].

α-and β-esterases
To measure the activity levels of α-and β-esterases, 100 μl of either α-or β-naphthyl acetate (Sigma-Aldrich, St. Louis, MO,USA) dissolved in acetone (CTR Scientific, Monterrey, N.L., Mexico), and phosphate buffer (KPO 4 , pH 7.2) was added to each well and incubated for 20 min at room temperature. Subsequently, 100 μl of fast blue (tetrazotized O-dianisidine (Sigma-Aldrich) dissolved in distilled water) was added to each well and the microplates incubated for another 4 min, before absorbance was read at 540 nm using a spectrophotometer (ASYS Hitech GmbH, Eugendorf, Austria).
The values obtained after subtracting the initial reading (T0) from the reading at 10 min (T10) were used for the statistical analyses.
Total protein content was determined by the Bradford method [56], which was used to correct the activity values for all of the enzymes evaluated [57]. Positive and negative controls were included for MFOs and esterases. The same volume of homogenate which was used in the respective assays was used in running the controls. For the α-and β-esterases, α-and β -naphthyl acetate solution was used as positive controls, respectively. Cytochrome-C (Merck, Darmstadt, Germany) solution was the positive control for the MFO assay.

Calibration curves
The calibration curves were generated using control solutions for each enzyme family. α-or β-naphthol (Sigma-Aldrich) was used for α-esterases and β-esterases, with a concentration range of 0.3 to 5 μg/μl and 0.5 to 4 μg/μl, respectively. Cytochrome C concentrations from 0.002 to 0.65 μg/μl were used for MFOs, and 0.1 to 6.5 μg/μl BSA solutions were used for total protein determination.
The absorbance values obtained for each enzyme were used to calculate the mean absorbance values per mosquito, which was then converted to enzyme activity after accounting for the homogenization volume, total protein content and the activity unit considered for each enzyme. Each value was multiplied by a conversion factor that was obtained from the calibration curves [57,58].
The mean activity values underwent an analysis of variance (ANOVA, P < 0.05) and Tukey's multiple comparison of means (P < 0.05). Normality of the variance was verified by the Bartlett test. Statistical analyses were carried out with GraphPad Prism v.7 (GraphPad Software, Inc, Version 6.01, La Jolla, CA, USA; https ://www.graph pad.com). Activity levels for all enzymes were determined in the parental population and in each deltamethrin-selected generation (F S0 -F S5 ).
The activity values corresponding to the 99th percentile of the New Orleans reference strain were calculated for each enzyme family. The average enzymatic activity was classified as non-altered (NA) when < 15% of individuals did not exceed the 99th percentile of the reference strain, incipiently altered (IA) if 15-50% of the individuals exceeded the 99th percentile, and altered (A) when this percentage exceeded 50%, according to the criteria established by Montella et al. [59]. Additionally, the LC 50 values and mean enzyme activity values of the parental generation and each selected generation (F S1 -F S5 ) underwent linear regression analysis. Correlation coefficients were calculated to determine the degree of association between the two variables.

Bioassays
Deltamethrin KC 50 and LC 50 were determined for the parental generation and in each selected generation (F S1 -F S5 ) as well as for the insecticide-susceptible New Orleans (NO) reference strain. Although some variation was observed, KC 50 and LC 50 generally increased with respect to the parental generation (F S5 vs F S0 ) (Tables 2, 3, Fig. 2).
When comparing the LC 50 values obtained from the parental generations to the NO strain, five of the populations (Merida, Progreso, Hoctun, Agua Dulce and Jose Cardel) were highly resistant to deltamethrin with an RR LC50 of 18-134-fold, and only one population (Hunucma) showed moderate resistance (RR LC50 of 6-fold) ( Table 3).
With respect to knockdown in the parental generation, RR KC50 was high in two populations, Merida with RR KC50 of 20-fold and Jose Cardel with RR KC50 of 14-fold. Three populations (Hoctun, Progreso and Agua Dulce) showed moderate resistance with RR KC50 ranging between 5-9fold, and one population (Hunucma) had a low RR KC50 of 3-fold (Table 2).
When the populations were selected with deltamethrin, the increase in LC 50 was not significant for F S1 and F S2 for the Merida population (based on overlapping 95% CIs). However, there was a significant increase in LC 50 from F S1 to F S5 in all other populations of Ae. aegypti (Table 3). For KC 50 , the increases were significant between all generations of selection (F S1 -F S5 ) for the Progreso, Agua Dulce and Jose Cardel populations but were not significant for the Merida and Hunucma populations for F S1 -F S2 and the Hoctun population for F S1 -F S3 ( Table 2).
When comparing the LC 50 and KC 50 values for the last deltamethrin-selected generation (F S5 ) vs the parental generation (F S0 ), we observed an increase of ~1.5-5.5 times for LC 50 and ~1.4-2.7 times for KC 50 in all populations (Fig. 2).
The results revealed that deltamethrin resistance increased in response to selection pressure (Tables 2, 3, Fig. 2). The knockdown resistance ratio (RR KC50 ) value increased from 20-fold to 38-fold from F S0 -F S5 in the Merida population, from 6-fold to 23-fold for Progreso, from 3-fold to 8-fold for Hunucma, from 5-fold to 10-fold in Hoctun, and slight increases were observed in the Agua Dulce and Jose Cardel populations with 9-12-fold and 14-22-fold differences, respectively. Lethal concentrations also increased between F S0 and F S5 , with the RR LC50 increasing in the Hunucma populations from 6-fold to 34-  32-52-fold and 41-71-fold, respectively). When analyzing the relative increases in RR KC50 and RR LC50 during the selection, the Hunucma population showed the highest single-generation increase (> 0.3-fold) from F S1 vs F S2 for KC 50 ; Hunucma showed also the highest increase in RR LC50 (2.83) in the first selected generation relative to the parental generation (F S1 vs F S0 ).

Realized heritability
The heritability of resistance (h 2 ) estimated over the generations of deltamethrin selection was highest in the Hoctun and Merida populations, with values of 0.948 and 0.900, respectively, for LC 50 . In the case of KC 50 , the highest values were obtained for the Progreso population (0.779) and Merida populations (0.519). The response to the selection (R) was also high for Merida and Progreso, for both KC 50 and LC 50 . The selection differential (S) was high in the Progreso population for LC 50 , but in the case of KC 50 , the highest values were obtained for Hoctun and Hunucma. The number of generations required for a 10-fold increase in LC 50 , (G), is the reciprocal of the response to selection (R) [47]. Thus, for the Progreso,  A total of 20 combinations of tri-locus genotypes were detected among the populations in the parental generation (F S0 ). Figure 3 shows the frequency of each of the 20 tri-locus genotype combinations. Triple homozygote resistant genotype (LL 410 /II 1016 /CC 1534 ) occurred in all parental (F S0 ) populations of Ae. aegypti with the highest frequency for the Merida population (0.63) and the lowest for Hoctun (0.06), and their frequency increased as generations were subsequently selected (Fig. 4).
Interestingly, the triple heterozygous genotype (VL 410 / VI 1016 /FC 1534 ) was only found in the populations of Hunucma, Hoctu and Jose Cardel, in which the triple wild-type genotype (VV 410 /VI 1016 /FC 1534 ) was also found.
In general, the frequencies of the three kdr alleles increased over the selected generations. The frequencies of L410 from the parental generation to the last selected generation increased from 0.78 to 0.92 in the Merida population, 0.65 to 0.92 in Progreso, 0.36 to 0.95  The frequencies of the C1534 allele increased from 0.59 to 1 in the Hoctun population, 0.64 to 1 in Hunucma, 0.89 to 1 in Jose Cardel and 0.97 to 1 in Progreso; the C1534 allele was fixed at 1.0 in the Merida and Agua Dulce populations in all generations. This allele reached fixation after the first selection with deltamethrin (F S1 ) in the population from Progreso and Jose Cardel, in the third selected generation for Hoctun and in the fourth selected generation for Hunucma.
The L410 and I1016 alleles were selected in almost the same proportion from F S0 to the F S5 in all populations   Figure 6 shows that the increase in the frequencies of the L410 alleles in response to the selection with deltamethrin correlates significantly (P < 0.05) with the increase in the frequencies of the I106 allele by 98% in the population of Merida, 96% in Hunucma, 95% in Hoctun and 95% in Agua Dulce. Additionally, in the populations of Hoctun and Hunucma we found that the allelic frequency of C1534 increased significantly (P < 0.05) with the increase in allele frequency of I1016 and L410 throughout the selection in the populations of Hunucma and Hoctun (Fig. 7). The frequency of the mutant alleles for each mutation was correlated with the KC 50 and LC 50 values across the generations of selection (Table 9). KC 50 values were significantly correlated with L410 frequencies in the Merida and Hoctun populations (P < 0.05); and with the frequencies of I1016 in the Merida, Hoctun and Jose Cardel (P < 0.05). Meanwhile, the frequencies of L410 and I1016 were significantly correlated with LC 50 values in the Merida and Hunucma populations as well as C1534 correlated with LC 50 values in Hunucma and Jose Cardel (P < 0.05). The presence and frequency of triple homozygous mutant individuals showed a significant correlation with KC 50 only in the Merida population and LC 50 in the Hunucma population (P < 0.05). Table 10 shows the mean enzyme activity levels detected in the biochemical assays in the parental generation and all deltamethrin-selected generations, where at least 30 individuals per population were analyzed. When comparing values from the parental generations with the 99th percentile of the New Orleans reference strain, five of the six populations analyzed (Merida, Progreso, Hunucma, Jose Cardel and Agua Dulce) showed altered enzyme activity for α-esterases (˃ 50% of the individuals exceeded the threshold of resistance), and Hoctun showed incipient altered enzyme activity (39% of the individuals exceeded the threshold). Merida, Agua Dulce and Jose Cardel showed altered enzyme activity for β-esterases (˃ 50%), and Progreso, Hunucma and Hoctun showed incipient enzyme alteration with percentages ranging between 16-43%. Altered activity of MFOs was detected in Progreso, Hoctun, Agua Dulce and José Cardel in 69-95% of the specimens analyzed. Only the Agua Dulce population showed altered enzyme activity for GSTs (71%) in the parental generation.

Biochemical assays
The results show variations in enzyme activity across populations over the generations post-selection. However, when comparing F S5 vs F S0 , α-esterases remained altered in the Progreso, Hunucma and Agua Dulce; β-esterases in Progreso and Agua Dulce, and an increase in the mean activity levels of GSTs was observed only in the Hoctun and Agua Dulce populations (P < 0.0001) ( Table 10). GSTs were the only enzyme family for which a significant correlation was detected between level of activity and LC 50 values, but only in the F S4 generation (r = 0.83, P < 0.05).

Discussion
Target-site insensitivity and increased metabolic activity are key mechanisms of pyrethroid resistance in mosquito populations worldwide. The control of adult mosquitoes in Mexico has been carried out mainly with pyrethroids since 1999, and deltamethrin has been in use since the beginning of 2000 [9]. In the present study, five of the populations analyzed were found to be highly resistant to the pyrethroid deltamethrin in the parental generation according to the LC 50 (RRLC 50 > 10-fold), and only one population showed moderate resistance (RRCL 50 5-10-fold). Likewise, high RRKC 50 values were detected in two populations, three showed moderate resistance, and one showed low resistance. These results suggest that deltamethrin resistance is well-established in these populations.
In response to further deltamethrin selection, both LC 50 and KC 50 increased over the generations, with the concentrations capable of knocking down or killing 50% of the population at least doubling by F S5 compared to F S0 . The highest heritability values (h 2 ≥ 0.90) for LC 50 were detected in the populations from Merida and Hoctun, followed by the Progreso population (0.61), suggesting high values of additive genetic variance. The highest value of h 2 for KC 50 was detected in the Progreso population (h 2 > 0.70), followed by the Merida population (0.52). A high h 2 value predicts a rapid response to selection, while a low value suggests a slow response to artificial selection [46].
The frequencies of the L410, I1016 and C1534 alleles increased with deltamethrin selection in all populations. When analyzing the frequencies in the last selected generation among the different populations, the frequencies of the L410 allele ranged between 0.80-0.98; this allele is able to significantly decrease the sensitivity of the sodium channel for both permethrin and deltamethrin [27]. The I1016 allele frequencies were 0.88-0.96 in F S5 and the C1534 allele reached fixation by F S5 in all populations analyzed. These results suggest that deltamethrin can select for C1534 more rapidly than I1016, corroborating the findings of Alvarez et al. [60].
The co-occurrence of I1016 and C1534 and high levels of pyrethroid resistance has been reported in other countries [61][62][63][64][65][66]. In Mexico, the co-occurrence of these two mutations and pyrethroid resistance has been also documented [34,50]. Vera-Maloof et al. [34] originally suggested that in order to present resistance to pyrethroids, the sequential evolution of both mutations was necessary. They considered it unlikely that I1016 had evolved independently due to the low fitness exhibited by haplotype I1016/F1534; instead, they hypothesized that C1534 first evolved by conferring a low level of resistance individually and that I1016 arose from haplotype V1016/C1534 and selected quickly due to the high level of resistance conferred by the double mutant haplotype. The sequential selection of F1534C and V1016I was later confirm by Chen et al. [31]. The V410L mutation was reported for the first time in the sodium channels of mosquitoes by Haddi et al. [27]. They demonstrated that this mutation reduced the sensitivity of the vgsc gene to both type I (i.e. permethrin) and type II pyrethroids (i.e. deltamethrin). The presence of this mutation has been demonstrated in Ae. aegypti from Mexico for more than 16 years and it has been shown that L410 is in greater linkage disequilibrium with I1016 than with C1534. Our results show that L410 and I1016 responded to the selection with deltamethrin in a similar way in populations where C1534 was fixed; or in Hoctun and Hunucma where the increase in C1534 was associated with the increase in L410 and I1016. This is consistent with the proposed sequential model where both V410L and V1016I might have occurred independently on a C1534 haplotype and then became cis to C1534 by recombination. Or alternatively and considering that the F1534C mutation was fixed in practically four of the six basal populations (F S0 ), the three mutations arose independently at very low frequencies, and then by two recombination events, came to occur in a cis arrangement [30].
Increased activity of mixed-function oxidases (MFOs) and esterase activity have been associated with pyrethroid resistance in Ae. aegypti [38,67]. Our results detected incipiently altered levels of MFO activity in four populations. Altered levels of α-esterase activity were detected in five populations, and at least four populations showed altered levels of β-esterase activity. However, there was no association between increased activities of these enzymes and phenotypic resistance, suggesting that these enzymes are not strongly associated with the metabolism of deltamethrin. However, contrary to the findings of Son et al. [68], who reported increased levels of MFOs in response to deltamethrin selection, we found that MFO activity decreased with successive selection. Given that the frequencies of the kdr mutations increased significantly with deltamethrin selection, the rapid selection of these mutations could have conditioned the activity of the enzymes in relation to deltamethrin resistance.
One criticism in our methodology is that we compare the enzymatic activity of each generation of deltamethrin in each of the populations with respect to the enzymatic activity of the susceptible New Orleans strain, when the most appropriate comparison should be with respect to the same population, in the same generation, without selection with deltamethrin.
In the metabolic resistance, three families of detoxifying enzymes, oxidases, esterases and glutathione transferases are mainly involved, all associated to confer  Table 9 Pearsonʼs correlation coefficients (r) (P-values in parenthesis) between the resistant allelic frequencies for L410, I1016, C1534 and with the frequencies of tri-locus genotypes vs KC 50  resistance through overexpression, increased insecticide metabolism or a greater affinity for the chemical [6,17]. An example of this contribution is reported by Lumjuan et al. [69] who determined an overexpression of GSTe2, GSTe5, GSTe7 and GSTE5-5 of epsilon glutathione transferase in individuals with resistance to DDT and pyrethroids. This fact is directly associated with the reduction in the susceptibility of these individuals to chemical control. Similarly, using quantitative trait loci (QTL) it was shown that temephos resistance is related to a QTL on chromosome two where a carboxylesterase cluster occurs; subsequent studies showed that there is an increase in the expression of these enzymes and therefore a possible intervention as a resistance mechanism [70]. Similar results were reported by Saavedra et al. [71], where the significant participation of an esterase marker discovered by QTL, CCEunk7o, and its relationship to resistance development in conjunction with mutations in the voltage-gated sodium channel was demonstrated. Alongside gene expression studies, several investigations focus on studying the levels of detoxifying enzymatic activity. Proof of this is the increase in activity suffered by multiple-function oxidases, which includes cytochrome P450, in individuals selected with deltamethrin for 15 generations compared to the unselected susceptible group [60]; however, as for this mechanism to be in itself the one that confers insecticide resistance is questionable since other mechanisms were described for these populations.
To verify that the enzymatic activity has a direct effect on the decrease in susceptibility, synergists are used, chemical compounds that inhibit detoxifying enzymes. By inhibiting the enzymatic activity of multiple function oxidases, esterases and glutathione transferases directly involved as a resistance mechanism, proof of this is reported by Bharati & Saha [72] where an increase in enzymatic activities of multiple-function oxidases (CYP450) and carboxylesterases was identified as a mechanism of resistance, by inhibiting such enzymes, the susceptibility to the pyrethroid deltamethrin and the carbamate propoxur was recovered in populations of Ae. aegypti from Bangladesh. A similar pattern was demonstrated in populations of Ae. aegypti resistant to DDT and pyrethroids from Selangor, Malaysia, where the inhibition of multiple function oxidases and glutathione-S transferases through an assay with synergists correlated with an increase in the activity of such enzymes as a possible mechanism of resistance [73].
Identifying the direct and exclusive participation of metabolic resistance to insecticides is a complex process due to the simultaneous presence of mechanisms such as target site mutations and elevated levels of detoxifying enzymes. Despite this, the occurrence of Ae. albopictus populations free of insecticide-resistant target site mutations has been reported, demonstrating the direct and exclusive participation of cytochrome P450 genes such as CYP6P12 and the overexpression of cuticular genes as the main resistance mechanisms [74].

Conclusions
The high levels of resistance associated with high frequencies of kdr alleles L410, I1016 and C1534 obtained through artificial selection suggest the important role of these mutations in resistance to deltamethrin. The role of metabolic resistance was less clear, as the activity levels of key enzyme groups appeared to increase or decrease without clearly relating to the selection with deltamethrin. These findings highlight the importance of monitoring the susceptibility status of mosquito populations before choosing insecticide products for vector control, especially in areas where resistance may already be present and further selection could quickly result in fixation of kdr mutations.