Parallel evolution of vgsc mutations at domains IS6, IIS6 and IIIS6 in pyrethroid resistant Aedes aegypti from Mexico

Aedes aegypti is the primary urban mosquito vector of viruses causing dengue, Zika and chikungunya fevers –for which vaccines and effective pharmaceuticals are still lacking. Current strategies to suppress arbovirus outbreaks include removal of larval-breeding sites and insecticide treatment of larval and adult populations. Insecticidal control of Ae. aegypti is challenging, due to a recent rapid global increase in knockdown-resistance (kdr) to pyrethroid insecticides. Widespread, heavy use of pyrethroid space-sprays has created an immense selection pressure for kdr, which is primarily under the control of the voltage-gated sodium channel gene (vgsc). To date, eleven replacements in vgsc have been discovered, published and shown to be associated with pyrethroid resistance to varying degrees. In Mexico, F1,534C and V1,016I have co-evolved in the last 16 years across Ae. aegypti populations. Recently, a novel replacement V410L was identified in Brazil and its effect on vgsc was confirmed by electrophysiology. Herein, we screened V410L in 25 Ae. aegypti historical collections from Mexico, the first heterozygote appeared in 2002 and frequencies have increased in the last 16 years alongside V1,016I and F1,534C. Knowledge of the specific vgsc replacements and their interaction to confer resistance is essential to predict and to develop strategies for resistance management.

Pyrethroids are the most common class of insecticides used to suppress adult populations of Aedes aegypti, the principal vector of dengue, chikungunya, Zika and yellow fever viruses. The lack of vaccines for most of these arboviruses results in a tremendous reliance on chemical control. Unfortunately, intense application of pyrethroids has resulted in pyrethroid resistance in Ae. aegypti populations worldwide [1][2][3][4][5] . A major mechanism underlying pyrethroid resistance is known as knockdown resistance (kdr), which is caused by mutations in the voltage-gated sodium channel (vgsc) 6 . Knowledge of the specific vgsc mutations that confer resistance is essential to predict the rise of pyrethroid resistance and to develop strategies for resistance management.
Globally, eleven vgsc mutations have been identified in Ae. aegypti and in most cases, have been linked to conferring some degree of pyrethroid resistance. These include G923V, L982W, I1,011 M and V1,016 G first identified in 2003 7 , I1,011 V and V1,016I in 2007 8 , D1,763Y in 2009 9 , S989P and F1,534C in 2010 10,11 , T1,520I in 2015 12 and V410L in 2017 13 . These kdr mutations are usually confined to specific geographic areas and co-occurrence of certain mutations is a common phenomenon sometimes associated with higher levels of phenotypic resistance 14 .
In Mexico, pyrethroid-resistant Ae. aegypti populations carry at least two vgsc mutations, V1,016I which is linked to permethrin survival 8 and F1,534C. Interestingly, F1,534C reduces permethrin binding to vgsc channels whereas V1,016I has no effect in pyrethroid binding 14  Recently, a novel mutation V410L in domain IS6 was identified in a pyrethroid resistant laboratory strain of Ae. aegypti from Brazil 13 . Alone or in combination with F1,534C, V410L reduced the sensitivity of mosquito sodium channels expressed in Xenopus oocytes to both type I (eg. permethrin) and II pyrethroids (eg. deltamethrin). Interestingly, authors did not detect this mutation in field populations from Pernambuco, Brazil, concluding the novel mutation may not yet be widespread in nature. In contrast, we identified high frequencies of V410L alongside V1,016I and F1,534C in a genome-wide association study of Ae. aegypti from Mexico. To extend this observation, we screened the frequencies of V410L in different temporal and geographical Ae. aegypti collections made over the last 16 years in Mexico. We compared the evolution and linkage disequilibrium of V410L alongside the V1,016I and F1,534C replacements, which occur in different domains IS6, IIS6 and IIIS6, respectively (Fig. 1). In addition, we determined the phenotype/genotype association in a field population exposed to permethrin and deltamethrin.   19 , L410 also increased in frequency from 2000 to 2016, noting that at each of these time points L410 and I1,016 alleles changed frequencies in parallel from 0.00 to 0.71 (Fig. 2).

V410L in
Pairwise linkage disequilibrium analyses were performed between locus 410-1,016, 410-1,534 and 1,016-1,534. Table 2 shows the linkage disequilibrium coefficients (R ij ), χ 2 and associated probabilities obtained between pairwise loci. Fourteen out of 25 collections had alleles segregating at loci 410-1,016 with R ij values ranging between 0.53-0.99 among collections; overall R ij was 0.96 (p = 0.0001). For loci 410-1,534, only five collections had mutant alleles segregating, and four were in significant linkage disequilibrium with R ij ranging from 0.33 to 0.99; the overall R ij coefficient was 0.76 (p = 0.0001). At loci 1,016-1,534, segregating alleles from six collections were in significant linkage disequilibrium, with R ij coefficients ranging from 0.31 to 0.99; overall R ij coefficient was 0.76 (p = 0.0001).
Association of V410L with pyrethroid resistance. L410 was present at a frequency of 0.69 in Viva Caucel mosquitoes used for our phenotype/genotype association study. We used a dose of permethrin or deltamethrin to discriminate three phenotypes: knockdown resistant, recovered and dead. The susceptible genotype at locus 410 was VV 410 , heterozygote was VL 410 and resistant was LL 410 . Table 3 shows the outcomes of mosquitoes carrying a specific genotype in terms of response to pyrethroid treatment. For permethrin, 53% of the resistant homozygotes (LL 410 ) were knockdown resistant, 40% recovered and only 7% died. In the heterozygote group (VL 410 ) 4% was knockdown resistant, 28% recovered and 64% died following permethrin exposure. For deltamethrin, 72% of the resistant homozygotes (LL 410 ) were knockdown resistant whereas 23% recovered and only 5% died. Note that the phenotype outcome is very similar between genotypes at loci 410 and 1,016. The most striking difference at locus 1,534, is that more than 92% of the heterozygotes died after exposure to permethrin or deltamethrin. Across all analyses, strong correlations were detected between phenotype and genotype.
Association of tri-locus genotypes with pyrethroid resistance. Because our results indicated that L410, I1,016 and C1,534 do not occur independently, we analyzed the phenotype outcome by tri-locus genotype combinations. In the Viva Caucel population, 13 tri-locus genotype combinations were identified. Figure 4 shows the distribution of knockdown resistant, recovered and dead mosquitoes for the eight most common tri-locus genotype combinations. The presence of resistant alleles in the tri-locus homozygote genotype is strongly associated with knockdown resistance and recovery for both permethrin and deltamethrin (Fig. 4a,b) 1,534 ) was associated with the dead phenotype for permethrin exposure but was associated with knockdown resistance and recovery in the deltamethrin exposure group.

Discussion
Different replacements at residue V410 have been reported in the vgsc of several pyrethroid resistant insect species. Specifically, V410L was associated with deltamethrin resistance in the common bed bug Cimex lectularis 29 . However, replacements V410M in the tobacco budworm 30    A recent study found that V410L and F1,534C occurred without V1,016I in a deltamethrin resistant laboratory strain originated from Rio de Janeiro 13 . In contrast, we show large linkage disequilibrium between L410 and I1,016, except among very few individuals collected early in 2002-2005; however, this genotype combination was no longer detected in following years. Whether L410 remains independent of I1,016 in Brazil will provide evidence of the mutations arising independently at a local level in the sequential model. However, I1,016 and C1,534 are already widespread in several regions of Brazil, and due to high migration rates among Ae. aegypti populations, we would expect I1,016 and C1,534 to recombine with L410-C1,534 in future years. However, an alternative scenario is that as in Mexico, L410 is already present at high frequencies in Brazilian collections previously genotyped with I1,016-C1,534 but simply has not yet been detected.
The selection of the triple homozygote resistant genotypes detected in our data suggests higher fitness of this genotype in the presence of pyrethroids. The role of L410 and C1,534 in conferring pyrethroid resistance was determined in Haddi et al. 2016 13 . L410 alone or in combination with C1,534 confers high levels of resistance,   however, it remains to be seen if it is fit in the field. We found 4 out of 1,176 individuals with L410 and C1,534 occurring independently (in the absence of I1,016), and this genotype became extinct in Mexican populations. In our phenotype and genotype studies, the triple homozygote resistant individuals had better survival (knockdown resistance and recovery) following either permethrin or deltamethrin exposure. One particular genotype, VL 410 /VI 1016 /CC 1,534 had different outcome depending of the specific pyrethroid, with this genotype associated with dead in mosquitoes following permethrin exposure. In contrast, this genotype was mostly associated with survival (knockdown and recovery) in mosquitoes exposed to deltamethrin. Apparently, the presence of heterozygotes at loci 410 and 1,016 was sufficient for deltamethrin survival. F1,534C is located in the PYR-1 receptor site and is responsible for reducing vgsc sensitivity to permethrin. Although residue 1,016 is located in the PYR-1 receptor site, only the V1,016G replacement occurring in Ae. aegypti from Asia reduces vgsc whilst the V1,016I replacement found in the Americas does not 16 . In contrast, V410L is located in DIS6 but does not form part of the PYR-2 receptor site. It has been suggested that the reduction of sodium channel sensitivity to permethrin and deltamethrin by V410L might result from changes in the gating properties of vgsc without inhibiting molecule docking 13 . Because pyrethroids prefer to bind to sodium channels in the open state, kdr mutations that deter the open state would counteract the pyrethroid effects 13 . In recent structural models, pyrethroids make multiple contacts with helices IIL45, IIS5, IIS6, and IIIS6, as well as IL45, IS5, IS6, and IIS6 that would maintain vgsc in an open state 17,32 . Simultaneous binding of pyrethroids to both PYR-1 and PYR-2 is thought to effectively prolong the opening of vgsc 14 . It is possible that co-occurrence of V410L and V1,016I, although in different receptor sites, provide fitness advantages in the presence of pyrethroids, thus favoring co-selection. The interaction of both mutations in electrophysiology experiments will address if the co-occurrence of these mutations is compensatory or synergistic in the presence of pyrethroids. Allele frequencies and linkage disequilibrium analysis. V410L allele frequencies were calculated from genotypic frequencies following Garcia et al. 18 . The 95% high density intervals (HDI 95%) around these frequencies were calculated in WINBugs2.0 following 1,000,000 iterations. Departure from Hardy-Weinberg expectations were expressed as Wright's inbreeding coefficient (F IS ) and a χ 2 test was used to test the hypothesis F IS = 0 (d.f. = 1). Pairwise linkage disequilibrium between alleles at loci 410 and 1,016 or between loci 410 and 1,534 was calculated using LINKDIS following Vera-Maloof et al. 19 . Composite disequilibrium between resistant alleles was tested and a χ 2 test determined if significant disequilibrium existed among alleles at both loci.

Detection
V410L association with pyrethroid resistance. To determine phenotype/genotype associations, we used the Viva Caucel (long −89.71827, lat 20.99827) collection from Yucatan, Mexico made in 2011. First, we calculated the permethrin and deltamethrin (Chem Service) lethal concentration that killed 50% of the mosquitoes (LC 50 ) in 3-4 days old adults of the F 3 generation. The insecticide treatment consisted of a 1 h exposure in an insecticide coated bottle, transfer of mosquitoes to recovery chambers and mortality scored at 24 h after treatment. We assessed the levels of permethrin and deltamethrin resistance in Viva Caucel relative to the New Orleans (NO) reference strain. The permethrin LC 50 was 47.9-fold higher than NO (26.5 µg vs 0.55 µg, respectively) and the deltamethrin LC 50 was 47.6-fold higher than NO (10.49 µg vs 0.22 µg).
Once the LC 50 was calculated, six to 10 groups of 50 mosquitoes at a time were exposed to 25 µg of permethrin or 3 µg of deltamethrin coated in 250 mL Wheaton bottles during 1 h. Immediately following exposure, active (knockdown resistant) and inactive mosquitoes were transferred to separate containers. To ensure correct categorization, we phenotyped mosquitoes 4 h after treatment. We observed activity and if the mosquitoes were capable SCIENTIFIC REPORtS | (2018) 8:6747 | DOI:10.1038/s41598-018-25222-0 of flight, they were scored as 'knockdown resistant' . In the inactive group we separated the newly recovered mosquitoes from the inactive mosquitoes and scored them as 'recovered' and 'dead' , respectively. Supplementary Table 3 shows the total number of mosquitoes exposed to permethrin and deltamethrin and the distribution between the three phenotypic categories.
A subsample of mosquitoes from each group was individually frozen; DNA was isolated by the salt extraction method 33 and resuspended in 150 µl of TE buffer (10 mM Tris-HCl, 1 mM EDTA pH 8.0). For the Viva Caucel mosquitoes exposed to permethrin we randomly selected 95 knockdown-resistant, 95 recovered and 95 dead mosquitoes. For deltamethrin we randomly selected 111 knockdown-resistant, 67 recovered and 92 dead mosquitoes. We conducted genotyping at locus 410 using the V410L melting curve system described above. For V1,016I and F1,534C genotypes, we used previously described methodologies 8,11 . A contingency table was used to test for association between the phenotypes (knockdown resistant, recovered and dead) and genotypes (mutant homozygote, wild type homozygote, and heterozygote) at each locus separately (410, 1,016 and 1,534) and for the 27 (3 × 3 × 3) tri-locus genotype combinations.
Disclaimer. The findings and conclusions in this paper are those of the authors and do not necessarily represent the official position of the Centers for Disease Control and Prevention.