Loss of pyrethroid resistance in newly established laboratory colonies of Aedes aegypti

Background Resistance to pyrethroid insecticides in Aedes aegypti has become widespread after almost two decades of the frequent use of these pesticides to reduce arbovirus transmission. Despite this resistance, pyrethroids continue to be used because they are relatively inexpensive and have low human toxicity. Resistance management has been proposed as a way to retain the use of pyrethroids in natural populations. A key component of resistance management is the assumption that negative fitness is associated with resistance alleles such that resistance alleles will decline in frequency when the insecticides are removed. At least three studies in Ae. aegypti have demonstrated a decrease in pyrethroid resistance once the insecticide has been removed. Methods/Principal findings The present study aims to evaluate variation in the loss of pyrethroid resistance among newly established laboratory populations of Ae. aegypti from Mexico. Eight field collections were maintained for up to eight generations, and we recorded changes in the frequencies of the mutations at the V1,016I locus and at the F1,534C locus in the voltage-gated sodium channel gene (VGSC). I1,016 and C1,534 confer resistance. We also examined resistance ratios (RR) with type 1 and 2 pyrethroids. Conclusions/Significance We demonstrate that, in general, the frequency of the Ae. aegypti pyrethroid-resistance alleles I1,016 and C1,534 decline when they are freed from pyrethroid pressure in the laboratory. However, the pattern of decline is strain dependent. In agreement with earlier studies, the RR was positively correlated with the frequencies of the resistance allele I1,016 and showed significant protection against permethrin, and deltamethrin, whereas F1,534C showed protection against permethrin but not against deltamethrin.


Introduction
After almost two decades of frequent pyrethroid use for Aedes aegypti (L.) control, widespread resistance now exists [1,2]. Despite this, pyrethroid use continues because these insecticides are relatively inexpensive and have low human toxicity. Resistance management has been proposed as a way to preserve the effectiveness of pyrethroids for Ae. aegypti control programs [2]. A key component of resistance management is the assumption that negative fitness is associated with resistance alleles such that the resistance alleles will decline in frequency when the insecticides are removed.
Laboratory strains of Ae. aegypti have shown a decrease in resistance once a pyrethroid is removed, thereby suggesting a fitness cost is associated with resistance. To date, three laboratory studies have evaluated the loss of pyrethroid resistance in Ae. aegypti. In Taiwan, a permethrin-resistant laboratory strain was maintained for 47 generations under permethrin pressure. Following 15 generations without exposure, a significant decrease was observed in the permethrin resistance ratio (RR) and resistance alleles in the Voltage-Gated Sodium Channel gene (VGSC) [3]. In Brazil, after 15 generations, the frequency of I1,016 decreased from 0.75 to 0.20 [4]. A study in Mexico showed a significant increase in the proportion of knockeddown mosquitoes ten generations after removal of pyrethroid exposure [5] but without a decrease in the frequency of I1,016 and C1,534 mutations in the VGSC. All three studies were done in laboratory cages.
The present study aims to evaluate the loss of pyrethroid resistance from eight collections of Ae. aegypti (six field collections from or near the city of Merida and two collections from Tapachula and Acapulco from southern Mexico). These collections were maintained without pyrethroid pressure for eight consecutive generations; during this time, we recorded changes in the frequencies of two mutations-VGSC I1,016 and C1,534-and analyzed the resistance ratios (RR) for permethrin (pyrethroid type 1) and deltamethrin (pyrethroid type 2).

Aedes aegypti field populations
In 2014, larvae were collected from eight public sites in Mexico. These larvae were reared to adults and then identified as Ae. aegypti. The adult females were then blood fed (citrated sheep blood-Colorado Serum Co., Denver Colorado) using an artificial membrane feeder, and the eggs were collected for shipment. The GPS coordinates and name abbreviations of the collection sites appear in Table 1. Eggs from Yucatan were collected from three urban sites in Merida and from three villages near Merida. Two additional collections were taken from Tapachula and Acapulco, in the states of Chiapas and Guerrero, respectively.

Establishment and maintenance of field populations
F 1 eggs were sent to Colorado State University. Egg papers were placed in a container with 2 L of tap water to promote development and hatching. Larvae were fed 2 mL of 10% (w/v) liver powder solution every other day. We transferred pupae to cups inside plastic cages for adult emergence. Larvae and mosquitoes were maintained in an incubator at 27-28˚C, 70-80% humidity and a photoperiod of 12 h light:12 h dark. Adults were fed with raisins and allowed access to tap water. Females were offered citrated sheep blood in artificial membrane feeders, every four days, to obtain eggs. Females laid their eggs on moistened filter papers. The eggs were allowed to mature for 48 h before being allowed to partially dry at room temperature; they were then stored in sealed plastic bags. Each collection was split at the F 1 larval stage into three groups of 500 to act as biological replicates. In each subsequent generation, �250 adult ♀ and 250 adult ♂ were used to maintain each of the three replicates.

Resistance ratio
We evaluated the changes in the resistance levels by determining the LC 50 for permethrin and deltamethrin at each of the eight sites in the F 3 , F 5 , and F 8 generations. Resistance ratios were obtained by dividing the LC 50 calculated for each site by the LC 50 calculated for the New Orleans Ae. aegypti (NO) susceptible reference strain. This colony was originally collected in New Orleans, Louisiana, by the Centers of Disease Control and Prevention and donated by Dr. William Brogdon. Our laboratory has maintained this colony since 2005 free of insecticide exposure. None of the knockdown resistance (kdr) alleles (I1,016 nor C1,534) are present in NO, and pyrethroid susceptibility is routinely confirmed by bottle bioassays. Commonly, the LC 50 for permethrin ranges between 0.4 and 0.6 ug/bottle, and for deltamethrin, between 0.09 and 0.15 ug/bottle. We examined the relationship between the RR and the frequencies of the I1,016 and C1,534 alleles. Analyses were performed in the F 3 , F 6 , and F 8 generations and then for all three generations combined. We used PROC CORR in SAS 9.4 to calculate Pearson's correlation coefficient and to test for significance. Genotyping V1,016I and F1,534C DNA was extracted at each generation (F 1 -F 8 ) from individual mosquitoes by the salt extraction method [6] and resuspended in 180 μL of TE buffer. To identify allelic variation, we used the allele-specific polymerase chain reaction (asPCR), followed by generation of a melting curve (CFX-96 BioRad), to identify the genotypes [7][8][9]. In each of the eight generations, we analyzed three replicates of 50 adult mosquitoes (~25♀ and 25♂) for each of the eight collection sites. Sample sizes were kept intentionally large to minimize the founder effect and genetic drift.

Allele frequencies and linkage disequilibrium analysis
We estimated allele frequencies in each of the eight generations from the genotypic frequencies (resistant allele frequency) = (2 x resistant homozygote + heterozygote) / (2 x sample size). Allele frequencies were compared among the replicates using a 2x3 contingency χ 2 test. Win-BUGS 2.0 [10] with 10 6 iterations was used to calculate the 95% high-density intervals (HDI 95%) around the allele and genotype frequencies. We used ggplot2 in R-3.5.1 to graph the data. In addition, we used LINKDIS [11] and χ 2 tests to calculate the pairwise linkage disequilibrium coefficients (R ij ) between alleles at loci 1,016 and 1,534 [12] according to the following: where Ci = Hobs (i)-p i 2 , Hobs (i) is the observed frequency of the i homozygotes, and N is the sample size, and p i and p j are the frequencies of the alleles at locus i = 1,016 and locus j = 1,534, respectively. T ij is the number of times that allele i and allele j occur in the same individual. A χ 2 test was performed to determine if significant disequilibrium exists among all alleles at 1,016 and 1,534. The statistic was calculated and summed over all two-allele interactions as follows:

Resistance declines in absence of pyrethroids
Permethrin and deltamethrin LC 50 was determined at generations F 3 , F 6 and F 8 (S1 Table). These values were compared against the NO susceptible colony to obtain the resistance ratios (RR). Fig 1 and S1 Table show the resistance ratios (RR) obtained for both pyrethroids at each of the generations evaluated. Permethrin resistance in the F 3 varied from 15-fold to 60-fold in six of the eight sites ( Fig 1A). The most resistant colony was Mer3 (60-fold) followed by Acp (45-fold) and Mer3 (39-fold). The remaining colonies had RR between 15-fold to 21-fold. At the F 6 generation, RR ranged between 3-fold to 50-fold. The most resistant sites were Mer3 (50 fold) and Ac (42-fold). Acp, Dz, Mer1 and Mer2 had RR from 17-fold to 25-fold. At the F 8, permethrin RR varied from 2-fold to 28-fold. In general, permethrin RR declined after eight generations in five sites (Ac, Acp, Co, Mer1 and Mer3) whereas two sites did not show changes in RR (Dz and Tap). Deltamethrin RR ranged from 18-fold to 108-fold in generation F 3 ( Fig 1B). The most resistant colonies were Co, Acp and Tap with RR higher than 86-fold. Ac and Dz had RR of 32-fold and 18-fold, respectively. Unfortunately, we did not have enough material to process Mer1, Mer2 and Mer3 at this generation. In generations F 6 and F 8 , RR ranged from 3-fold to 27-fold. In general, we observed a decline in RR in five sites (Ac, Acp, Co, Dz and Tap). Mer1 and Mer 2 showed a slight decrease from the F 6 to F 8 .

Allele frequency of I1,016 declines in the absence of pyrethroids
We determined the V1,016I and C1,534 genotypes of 9,563 mosquitoes from eight sites in southern Mexico (Table 1) over eight generations in absence of pyrethroids (F 1 -F 8 ). The genotype counts and allele frequencies appear in S2 and S3 Tables, respectively. Allele and genotype frequencies were analyzed separately for each locus. We show a general decline of the resistance alleles in Figs 2 and 4; and the different tendencies of each of the genotypes in Figs 3 and 5. Fig 2A plots the I1,016 frequencies in all three replicates in the eight different collections over eight generations. Initial allele frequencies in the F 1 varied from 0.35 in Dz to 0.83 in Acp, most of the sites had frequencies ranging between 0.59 to 0.73. In general, the I1,016 allele frequencies were statistically uniform among the three replicates, with five exceptions that are indicated with an asterisk (AcF 2 , CoF 2 , CoF 8 , DzF 2 , and Mer3F 7 ) (Fig 2A). Fig 2B plots the mean I1,016 frequencies among all three replicates and the 95% high-density intervals (HDI 95%). The Pearson correlation coefficients between the I1,016 frequencies and the generation number and the associated significance appear at the bottom of each graph. The correlation between the I1,016 allele frequency and generation number was negative in each of the eight collections, though statistical significance varied from highly significant in four out of the eight sites (Ac, Acp, Co, and Mer2) ( Fig 2B) to not significant in sites Dz, Mer1, and Tap. The frequency of I1,016 in Co, Mer3, and Tap declined initially and then, surprisingly, increased in the last 3-4 generations. In general, I1,016 declined in frequency over eight generations; however, the rate and pattern of the decline varied greatly among the collections.

Genotype frequencies at 1,016 indicate lower fitness of the resistanthomozygote
To understand the fitness of each genotype, we examined the fate of the three genotypes separately over time (Fig 3). Because the knockdown resistance (kdr) is a recessive trait (only homozygote resistant genotypes survive exposure to pyrethroids), we expect a strong negative fitness cost for homozygote-resistant individuals (II 1,016 ) in the absence of pyrethroids. We show in Fig 3C that a decline in II 1,016 occurs over eight generations, and Table 2 shows negative correlation coefficients for all sites (r = -0.5452, P < 0.0001), confirming a negative fitness for the II 1,016. Additionally, we show that the susceptible-homozygote or wild-type genotype (VV 1,016 ) increased in frequency over the generations with positive correlation coefficients (r = 0.4497, p = 0.0002). This suggests that the susceptible homozygote genotype has greater fitness in the absence of insecticides. Interestingly, we did not observe changes in the frequency of the heterozygotes over time (r = 0.1462, P = 0.2489), indicating that no negative fitness cost occurs for VI 1,016 in laboratory conditions. Therefore, the resistant allele will prevail in heterozygous individuals for several generations in the absence of insecticides.

Allele frequency of C1,534 declines in the absence of pyrethroids
We obtained the genotypes and allele frequencies for F1,534C in each of the 9,563 mosquitoes for which the V1,016I genotype frequencies had been determined (S2 and S3 Tables). Fig 4A shows the resistant allele C1,534 frequencies in the eight collection sites over eight generations. Initial allele frequencies in the F 1 ranged from 0.67 in Dz to 0.98 in Acp and Tap. The C1,534 frequencies differed between replicates in 15 of the 64 comparisons. The correlation coefficients between the C1,534 allele frequencies and generation number and their significance are at the base of each of the graphs in Fig 4B. The correlation coefficients between the C1,534 frequency and generation number were negative and significant in four out of the eight collection sites for the F 1 to F 8 generations (Acp, Co, Mer1, and Mer2) (Fig 4B). The correlation was not significant for sites Ac, Dz, Mer3, and Tap. In general, C1,534 declined in frequency over the eight generations; however, the rate and pattern of the decline varied greatly among the collections.

Genotype frequencies at 1,534 indicate lower fitness of the resistanthomozygote
The genotype frequencies at F1,534C are provided in Fig 5, and the correlation coefficients between genotypes and generations are provided in Table 3. Five correlation coefficients for FF 1,534 were positive, with two being significant; five of the correlation coefficients for FC 1,534 were positive, with three being significant; and five of the CC 1,534 were negative, with four being significant. Across all collection sites, a positive correlation (r = 0.3399, p = 0.006) existed between the homozygote susceptible FF 1,534 and the generation number. The frequencies of FC 1,534 did not change significantly over the generations (r = 0.2078, P = 0.0994), and as expected for a genotype that confers a lower fitness, the CC 1,534 genotypic frequencies decreased significantly over the generations (r = -0.3024, P = 0.0152), suggesting a fitness cost associated with the resistant-homozygote genotype.

Linkage disequilibrium
We performed pairwise linkage disequilibrium analyses between the alleles in V1,016I and F1,534C. Table 4 lists the linkage disequilibrium coefficients R ij , χ 2 , and the probability value obtained between pairwise loci. R ij is distributed from -1.00 (mutations occur on opposite chromosomes-trans) to 0.00 (mutations occur independently) to 1.00 (both mutations on the same chromosome-cis); therefore, R ij provides a standardized measure of the disequilibrium.
Alleles segregated in 57 of the 64 collections (Table 4). Forty-seven of the 57 collections exhibited significant linkage disequilibrium, with the R ij values ranging between 0.15 and 0.85 among the collections. This result suggests that the resistance alleles I1,016 and C1,534 occur together more often than expected by independent, random segregation. In general, and in agreement with earlier studies [12,13], alleles in V1,016I and F1,534C were in linkage disequilibrium.
The four possible haplotypes resulting from V1,016I and F1,534C are shown in Fig 6. Table 5 displays the correlation and the significance between the haplotype frequencies and generations. The frequency of the susceptible V 1,016 /F 1,534 haplotype increased over the generations (r = 0.323, p = 0.009). The frequency of the V 1,016 /C 1,534 haplotype remained relatively constant (r = 0.140, p = 0.271). The frequency of the resistant I 1,016 /C 1,534 haplotype decreased over time (r = -0.516, p <0.0001) across all collection sites. The frequencies of the I 1,016 /F 1,534 haplotypes were consistently low (r = 0.091, p = 0.474) across generations. This same trend has been noted in two previous studies [5] [13], suggesting that low fitness may occur in any mosquito in which I1,016 co-occurs with F1,534. Temporal analysis of di-locus genotypes Fig 7 shows the frequency of the nine di-locus genotype combinations (three genotypes at two loci), and Table 6 lists the correlation between frequencies of each di-locus genotype and the generation number without pyrethroid exposure. We estimated the frequencies of the nine genotype combinations in 9,563 mosquitoes. However, the di-locus coefficient among the nine graphs was only significant for the wild type susceptible VV 1,016 /FF 1,534 and the dual resistant II 1,016 /CC 1,534. The correlation between generation number and VV 1,016 /FF 1,534 was positive (r = 0.3376, P = 0.0064), indicating an increase in susceptible di-locus genotypes in the absence of insecticides. The correlation of II 1,016 /CC 1,534 was negative (r = -0.5465, P < 0.0001), indicating a decline. The low frequencies of VI 1,016 /FF 1,534 , II 1,016 /FF 1,534 , and II 1,016 /FC 1,534 in Fig  7B, 7C and 7F are consistent with the hypothesis that low fitness may occur in any mosquito in which I1,016 co-occurs with F1,534. However, Fig 7E identifies

Association between resistance alleles and resistance ratios
Pearson correlation coefficients and their significance were calculated between the frequencies of I1,016 or C1,534 alleles and resistance ratios for permethrin or deltamethrin as determined   Table 7 indicates that all correlations were positive in F 3 and F 8 , but none were significant. All correlations were positive in F 6 , and three were significant. When the results from all three generations were combined, all correlations were positive, and three were significant. Thus, in general, a weak but consistently positive correlation existed between the frequency of the I1,016 or C1,534 alleles and resistance ratios as determined by bioassay (Fig 8). This correlation also was noted in an earlier study in which I1,016 showed significant protection against permethrin and deltamethrin, whereas F1,534C showed protection against permethrin but not against deltamethrin [14,15]. The expression of C1,534 in Xenopus oocytes and exposure to both pyrethroids demonstrated that the resistant amino acid substitution in C1,534 is sensitive to permethrin but not to deltamethrin, which is consistent with our results.

Discussion
In this study, we established eight colonies of Ae. aegypti from the field and maintained them in a pyrethroid-free environment in the laboratory for eight generations. We demonstrated that, in general, the frequency of the Ae. aegypti pyrethroid-resistance alleles I1,016 and C1,534 decline when released from pyrethroid pressure in the laboratory (Figs 2 and 4, Tables  3 and 5). However, the pattern of decline appeared to be strain dependent, with some strains showing a steady rate of decline (Ac, Acp, and Mer2 in Fig 2; Acp and Mer2 in Fig 4), some showing a shallow decline (Mer1 and Tap in Fig 2; Acp, Mer2, and Co in Fig 4), and others displaying no net change (Dz in Fig 2; Ac, Mer3, and Tap in Fig 4).
A more surprising result was that, in Co and Mer3, the frequencies of I1,016 increased following a precipitous drop. Likewise, the frequencies of C1,534 in Mer1, Mer2, and Co increased after a drop. Such a change might occur if deleterious or lethal recessive mutations are linked to the susceptible allele that became homozygous through continuous inbreeding. However, this theory fails to explain why I1,016 increased simultaneously in all three replicates of Co.
We do not suggest, nor do we know whether the selection pressure in indoor cage studies is the same as or even correlated with outdoor selection pressure. This study only indicates that the loss of pyrethroid resistance is unlikely to follow a smooth linear or exponential decline for any one of a number of reasons. Nor should we expect the decline in resistance to be consistent among collections. Epistatic interactions between alleles may cause nonlinear trends in allele frequencies. Much depends on the genetic background of each population; some populations could take much longer to lose resistance, while others may do so much more rapidly.
A variety of possibilities exist to explain this variance in the gene frequency trajectories among the collections, but this heterogeneity is unlikely to have arisen from small sample sizes. We analyzed three replicates of 50 adult mosquitoes for each of the eight collection sites.

PLOS NEGLECTED TROPICAL DISEASES
Furthermore, 500 adults were used to generate the next generation of eggs for each replicate. The 95% HDI remained narrow in all graphs in Fig 2B and Fig 4B. Initial conditions may affect the shape of the curve. For example, a curve that begins with initial frequencies close to 1 (Mer3 and Tap in Fig 4A) would begin to decline much later than a curve that begins at 0.6 (Mer1 and Mer2 in Fig 2A). Metabolic resistance may account for much of this heterogeneity. A quantitative trait loci mapping study (QTL) [9] reported that 58.6% of the variation in knockdown could be accounted for by I1,016, but that a number of different QTL located throughout the genome accounted for the remainder of the variation. Saavedra-Rodriguez et al. used the "Aedes Detox" microarray [16] and showed an inverse relationship between I1,016 frequencies and the number of differentially transcribed metabolic genes [17]. Table 5 displays the correlation and the significance between haplotype frequencies and generation. The low frequencies of VI 1,016 /FF 1,534 , II 1,016 /FF 1,534 , and II 1,016 /FC 1,534 noted in this (Fig 7B, 7C and 7F) and in two previous studies [5] [12] suggest that low fitness may occur in any mosquito in which I1,016 co-occurs with F1,534. This finding, and because C1,534 historically appears before I1,016, suggests that the evolution of the mutations was sequential. If Deltamethrin RR appears to be correlated with the I1,016 allele frequency but not with the C1,534 allele frequencies, whereas permethrin RR correlates with both allele frequencies ( Table 7). We speculated that other resistance mechanisms are present that could drive resistance to deltamethrin, such as detoxifying enzymes or other mutations in VGSC that have not been identified. However, permethrin resistance appears to be driven heavily by both resistance alleles.
Interestingly, the I1,016 mutation was functionally expressed in Xenopus oocytes and did not show an alteration in the sodium channel sensitivity to either of these pyrethroids [14]. Our data indicate a positive correlation of I1,016 allele frequency with the RR of both pyrethroids. Clearly, a greater understanding is required concerning the role that I1,016 plays in the resistance of Ae. aegypti. This includes the possibility that V410L [15] may be the residue to which the pyrethroid actually binds. Since many field populations are resistant to pyrethroids, and few insecticides are available on the market, negative fitness is beneficial for vector control, providing an opportunity for the use of alternative insecticides as pyrethroid resistance is lost. Pyrethroids could then be "saved" to control susceptible populations during a disease outbreak. However, the results in this study suggest that even if dominant wild-type alleles (V1,016 and F1,534) increase in a population in the absence of insecticides, recessive resistance alleles will be hidden and maintained in heterozygotes. Therefore, the initial susceptibility of populations (before the introduction of a pyrethroid) will never be recovered. Furthermore, recent work in Sao Paulo State in Brazil showed that resistance alleles persist in natural populations for at least 11 years [1].
Supporting information S1