Permethrin resistance in Aedes aegypti: Genomic variants that confer knockdown resistance, recovery, and death

Pyrethroids are one of the few classes of insecticides available to control Aedes aegypti, the major vector of dengue, chikungunya, and Zika viruses. Unfortunately, evolving mechanisms of pyrethroid resistance in mosquito populations threaten our ability to control disease outbreaks. Two common pyrethroid resistance mechanisms occur in Ae. aegypti: 1) knockdown resistance, which involves amino acid substitutions at the pyrethroid target site—the voltage-gated sodium channel (VGSC)—and 2) enhanced metabolism by detoxification enzymes. When a heterogeneous population of mosquitoes is exposed to pyrethroids, different responses occur. During exposure, a proportion of mosquitoes exhibit immediate knockdown, whereas others are not knocked-down and are designated knockdown resistant (kdr). When these individuals are removed from the source of insecticide, the knocked-down mosquitoes can either remain in this status and lead to dead or recover within a few hours. The proportion of these phenotypic responses is dependent on the pyrethroid concentration and the genetic background of the population tested. In this study, we sequenced and performed pairwise genome comparisons between kdr, recovered, and dead phenotypes in a pyrethroid-resistant colony from Tapachula, Mexico. We identified single-nucleotide polymorphisms (SNPs) associated with each phenotype and identified genes that are likely associated with the mechanisms of pyrethroid resistance, including detoxification, the cuticle, and insecticide target sites. We identified high association between kdr and mutations at VGSC and moderate association with additional insecticide target site, detoxification, and cuticle protein coding genes. Recovery was associated with cuticle proteins, the voltage-dependent calcium channel, and a different group of detoxification genes. We provide a list of detoxification genes under directional selection in this field-resistant population. Their functional roles in pyrethroid metabolism and their potential uses as genomic markers of resistance require validation.


Introduction
The mosquito Aedes aegypti is the primary urban vector of three globally important arboviral diseases-dengue, Zika, and chikungunya fever-for which vaccines and effective pharmaceuticals are still lacking. The only available strategy to suppress these arboviral outbreaks is to reduce vector populations. Control of Ae. aegypti is challenging and is further compromised by widespread pyrethroid resistance [1][2][3][4]. Heavy use of pyrethroid space sprays-due to their strong human safety profile-has created immense selection pressure for resistance [1][2][3][4]. This resistance is primarily under the control of the voltage-gated sodium channel (VGSC) and enhanced metabolism by detoxification enzymes.
Amino acid replacements at VGSC-the target site of pyrethroids-confer resistance to knockdown (kdr-mutations) [4,5]. Approximately 12 nonsynonymous substitutions in the VGSC gene (VGSC) are associated with pyrethroid resistance in Ae. aegypti. Three kdr-mutations are common in Latin America, including V410L, V1016I, and F1534C [6][7][8]. The role of V410L and F1534C to confer pyrethroid resistance was recently confirmed in electrophysiological assays [8,9]. Although V1016I alone had no effect on the channel sensitivity to permethrin or deltamethrin, it enhanced the F1534C-mediated resistance to both pyrethroids [10].
Additional mechanisms of pyrethroid resistance are conferred by enhanced insecticide metabolism (or sequestration) by three enzyme systems, the carboxyl/choline/esterases (CCE), glutathione-s-transferases (GST), and cytochrome P450 monooxygenases (CYP) [11]; reduced penetration of insecticides through the cuticle [12]; and behavioral avoidance [13]. The extent to which these different mechanisms contribute to the overall resistance phenotype seems to vary [5,6,14,15]. Previous studies in Ae. aegypti from Mexico showed two major quantitative trait loci (QTLs) controlling permethrin resistance [5]. One corresponded to the nonsynonymous mutation V1016I in the VGSC and the esterase CCEunk70. Additional QTLs contained several CYP genes of relatively minor effect. These results confirmed that target site insensitivity explained almost 60% of permethrin resistance but that other genes dispersed throughout the genome also contributed to the survival of mosquitoes following permethrin exposure.
A common methodology to diagnose pyrethroid resistance in mosquito populations is the bottle bioassay [16]. The major route of intoxication is through tarsal contact during exposure times that range from 30 min to 2 h, depending on the pyrethroid concentration. At the endpoint of the bioassay, two phenotypes are discriminated: knockdown or resistant mosquitoes. Intriguingly, once the knocked-down mosquitoes are removed from insecticide exposure, recovery rates from 20 to 60% have been reported in the literature [17]. Previous studies have shown that kdr mosquitoes often carry resistant homozygous genotypes at three loci in VGSC (V410L, V1016I, and F1534C). In contrast, recovered and dead mosquitoes often carry heterozygous or wild-type homozygous genotypes at these loci [5,7,15]. Two different studies showed that 42 and 32% of the recovered mosquitoes carry the V1016I heterozygous genotype [5,15], suggesting partial protection during recovery, however, the remaining recovered individuals (>58%) must rely in mechanisms other than kdr-mutations to recover.
In this study, we aim to identify genomic differences associated with kdr, recovered, and dead phenotypes in mosquitoes following permethrin exposure in the laboratory. We used a F 1 offspring of a pyrethroid-resistant Ae. aegypti field population from Tapachula, Mexico. Our hypothesis is that SNPs associated with kdr will occur at insecticide target site or cuticle genes, whereas recovery will be associated with genes linked to insecticide detoxification mechanisms. The identification of such SNPs will improve our understanding of the mechanisms of resistance associated with kdr, recovered, and dead in field pyrethroid-resistant mosquito populations.

Results
We exposed 401 mosquitoes for 1 h in bottles coated with 15 μg of permethrin. This permethrin concentration allowed us to discriminate three different phenotypes with significant sample size in a heterogeneous pyrethroid-resistant population from Tapachula. Fig 1 shows the three phenotypes discriminated by this concentration and time of exposure: 1) kdr (n = 58, 14.5% of total), 2) recovered (n = 130, 32.4% of total), and 3) dead (n = 213, 53.1% of total).
Six genomic libraries, consisting of two biological replicates of pooled mosquitoes (n = 25) exhibiting the kdr, recovered, or dead phenotypes were prepared. The cost of library enrichment using an exome-capture hybridization prevented us to process a third biological replicate of each phenotype. By using two replicates, we obtained between 112 and 129 million reads across the six pair-end libraries. Thus, sequencing coverage ranged from 196-fold to 288-fold. After removing repetitive DNA (coverage > 1000) and sites with fewer than 25 reads, we identified 30-35 million common sites among the three pairwise comparisons: 1) kdr vs recovered, 2) recovered vs dead, and 3) kdr vs dead. Between 1.69 and 2.3 million polymorphic sites (SNPs) were identified among the pairwise comparisons. Alternate nucleotides were defined as those differing from the reference genome. The frequency of the alternate nucleotide at each SNP was subjected to a contingency χ 2 analysis and then assigned a genetic association value (-log 10 (p value)), referred to as the "LOD". A Benjamini-Hochberg correction [18] for false discovery rate (FDR) was applied to SNPs at each chromosome separately using an α = 0.01. The LOD cutoff values ranged from 3.17 to 3.37 between the pairwise comparisons, resulting in 12,209 significant SNPs in the kdr vs recovered, 11,472 in the recovered vs dead, and 13,011 in the kdr vs dead comparison. The minimum and maximum LODs were 3.1 and 37, respectively. The mean LODs among SNPs ranged from 5.1 to 5.4, and the 95 quantiles, from 8.5 to 9.18 (Table 1). Approximately 55% of these SNPs occurred at intergenic sites, 7.2% at 3'-UTR sites, 7.5% at 5'-UTR sites, 36.9% at synonymous coding sites, 9.6% at nonsynonymous coding sites, 37.8% at intron sites, and 0.6% in noncoding RNA.
The mean LODs for SNPs belonging to three gene categories associated with insecticide resistance (cuticle, detoxification, and insecticide target sites) are shown in Table 1. For the cuticle category, the mean LODs were not significantly different between the phenotypes (F = 1.67, p value = 0.18). However, the mean LODs were significantly different for the target site category (F = 50.5, p value = 2e-16), in which the kdr vs recovered and kdr vs dead had mean LODs of 8.69 and 14.3, respectively, while the recovered vs dead comparison had a mean  Table 1. Mean and standard error (SE) of LOD values (-log 10 (p value)) assigned to SNPs differing between kdr, recovered and dead Aedes aegypti exposed to permethrin. The number of SNPs (N) and the LOD mean and SE for three categories of genes associated with insecticide resistance are shown separately. LOD of 5.19. In the detoxification category, we found significant differences between mean LODs (F = 6.76, p value = 0.001), with the kdr vs recovered and kdr vs dead explaining this difference.
The distributions of the SNPs by LOD and relative physical position across chromosomes for each of the pairwise comparisons are shown in Fig 2. Interestingly, the kdr vs recovered and the kdr vs dead comparisons showed a cluster of SNPs with high association values (> 95 quantiles) in chromosome 3 (Fig 2A and 2C). These clusters consist of SNPs in VGSC, the major pyrethroid target site. In contrast, the recovered vs dead comparison lacked this cluster of SNPs (Fig 2B), validating the role of VGSC mutations in kdr but not in recovery.

SNPs that differ between kdr and recovery
The mosquitoes included in this comparison survived the exposure to permethrin; however, some mosquitoes exhibited knockdown resistance at 1 h of exposure (kdr), whereas others survived by recovering from initial knockdown during the 4-h observation. A total of 12,209 significant SNPs resulted in this comparison (2242, 4520, and 5447 in chromosomes 1, 2, and 3, respectively). Fig 2A shows the distribution of SNPs by their LOD values and relative positions in the genome. We purposely highlighted the SNPs associated with insecticide target site, detoxification, and cuticle genes.
The 50 mosquitoes used to generate the kdr, recovered, and dead libraries were individually genotyped for V410L using an allele-specific PCR. Ninety percent of the kdr mosquitoes were resistant homozygotes (V410L/V410L), whereas 10% were heterozygotes (V410L/V410), and Distribution of SNPs associated with kdr, recovered, and dead Aedes aegypti exposed to permethrin. The relative physical position is based in Ae. aegypti AaegL5 assembly. LODs correspond to the-log 10 (p value) obtained in a chi square test that compared the proportion of the alternate allele between pairwise comparisons. A) kdr vs recovered and B) recovered vs dead and C) kdr vs dead. SNPs located in insecticide target site (yellow), detoxification (red), and cuticle genes (green) are highlighted. https://doi.org/10.1371/journal.pgen.1009606.g002

PLOS GENETICS
0% were wild-type homozygotes (V410/V410). In contrast, 8% of the recovered were resistant homozygotes, 80% were heterozygotes, and 12% were wild-type homozygotes. None of the dead mosquitoes were resistant homozygotes, 36% were heterozygotes, and 64% were wildtype homozygotes at this locus. A 3 x 3 contingency table showed significant differences between observed and expected genotypes at each phenotype (χ 2 = 168.8, df = 4 and p value = 1.8e-35). Fig 3 shows the Pearson residuals between the phenotypes and genotypes at locus V410L. Positive residuals in blue indicate a positive association between kdr and V410L/V410L and dead with V410/V410. Interestingly, recovery was positively associated with heterozygotes (V410L/V410).

Cuticle SNPs associated with kdr
Changes in the cuticle proteins can result in reduced penetration by insecticides. Approximately 88 SNPs were identified in this category. The eight nonsynonymous SNPs located in seven genes are shown in Table 3. The highest LODs occurred in cuticle protein CP14.6 (LOC5577605, LOD = 6.46), in an uncharacterized cuticle protein (LOC5572415, LOD = 6.02), and in cuticle protein 8 (LOC5565392, LOD = 5.89).

Detoxification SNPs associated with kdr
Two hundred ninety-five SNPs located in 91 genes differed significantly in the kdr vs recovered comparison. The 41 nonsynonymous SNPs were located in 31 detoxification genes Table 2. Nonsynonymous SNPs at insecticide target sites genes associated with kdr, recovered, or dead Aedes aegypti exposed to permethrin. Gene identification, chromosome (Ch), SNP site and amino acid residue (based in the AaegL5 genome assembly), LOD = -log 10 (p value vector base gene identification (from www.vectorbase. org). Voltage-gated calcium channel (VGCC), voltage-gated sodium channel (VGSC), gamma-aminobutyric acid type B receptor subunit 2 (GPRGBB3), acetylcholinesterase/hydrolase (ACE), and G protein-activated inward rectifier potassium channel 3 (KCNJ3).

Gene
Gene  In the kdr vs recovered pairwise comparison, SNPs in which the alternate nucleotide frequency was higher in the kdr group were assumed to be potentially beneficial in the presence of insecticides and therefore under directional selection. In contrast, those SNPs where the alternate nucleotide frequency was higher in the recovered group were suspected to be nonbeneficial for the kdr phenotype and under neutral or purifying selection.  Pearson residuals between three genotypes at V410L in VGSC and the three phenotypes: kdr, recovered, and dead in Aedes aegypti mosquitoes exposed to permethrin. Genotypes were determined by an allele-specific PCR to detect three genotypes: homozygote resistant = V410L/V410L, heterozygote = V410L/V410), and wild-type homozygote = V410/V410. Positive residuals in blue specify a positive association between the corresponding row and column variables. Negative residuals are in red, implying a negative association between the corresponding row and column variables.
https://doi.org/10.1371/journal.pgen.1009606.g003 was categorized as low (probably novel mutations) when the frequency ranged from 0 to 0.4, moderate from 0.4 to 0.8, and high from 0.8 to 1.0. Example of genes under directional selection for kdr include GSTD3, CCEae3O, CCEunk7O, and CCEjhe2O, whereas CYP4C51 had high frequency in the recovered group, suggesting that this SNP is under purifying selection. In the recovered vs dead group, SNPs in which the alternate nucleotide frequency was higher in the recovered group were assumed to provide a beneficial protection when no kdr-mutations are present and therefore under directional selection. Sites where the alternant SNP was higher in the dead group were assumed to be under purifying selection. Table 3. Nonsynonymous SNPs in cuticular protein genes associated with kdr, recovered, or dead Aedes aegypti exposed to permethrin. Gene identification, chromosome (Ch), SNP site, amino acid residue (based in the AaegL5 genome assembly), LOD = -log 10 (p value) and vector base identification (www.vectorbase.org).

Gene
Gene

SNPs that differ between recovered and dead
During exposure to permethrin, 86% (343 out 401) of the mosquitoes experienced knockdown. In this knockdown group, 37% recovered (130 out of 343), whereas 62% died (213 out of 343). When comparing the genomes of the recovery and dead phenotypes, we identified

Top nonsynonymous SNPs associated with recovery
The nonsynonymous SNPs associated with recovery are shown in S2 Table. In chromosome 1, we identified 217 nonsynonymous located in 161 genes. The SNP with the highest LOD was a Y309C in putative inositol mono-phosphatase 3 (LOC110677190, LOD = 12.4), followed by a R435Q in the uncharacterized LOC5564494 (LOD = 11.6), and a L505V in uncharacterized protein F54F2.9 (LOC5579009, LOD = 11.34). In chromosome 2, 388 nonsynonymous mutations in 279 genes were identified (S2 Table).

Insecticide-target site SNPs associated with recovery
We identified 88 SNPs in 22 genes in this category ( Table 2). Six nonsynonymous SNPs were found in four genes, including the acetylcholine receptor subunit alpha-like (LOC5575838, LOD = 6.7), GPRGBB3 (LOD = 7.6 and 7.1) and the voltage-dependent calcium channel type A subunit alpha-1 (VGCC, LOC5564339) with an LOD of 10.9. The S679T and V410L at VGSC had LODs of 3.8 and 4.6, respectively.

Cuticle SNPs associated with recovery
Approximately 112 SNPs located across 50 genes were identified in this category. The 12 nonsynonymous mutations were located in seven genes ( Table 3). The highest LOD occurred in the larval cuticle protein A2B (LOC5577372, LOD = 8.8 and 6.2) and endocuticle structural protein SgAbd-6 (LOC5570308, LOD = 5.8) in chromosome 3.

SNPs that differ between kdr and dead
This comparison identified SNPs that would be hypothetically selected if kdr were the only mechanism under selection, ignoring recovery as a mechanism of survival. A total of 12,381 SNPs were associated with kdr. In chromosome 1, the nonsynonymous mutations with the highest LOD values were uncharacterized LOC5576488, nose-resistant to fluoxetine protein 6 (LOC5576756), and protein msta (LOC5579017) with LOD values of 15.1, 13.1, and 11.1, respectively (S3 Table). In chromosome 2, vitamin K-dependent gamma-carboxylase (LOC110675715, LOD = 14.1) and a fatty acid synthase (LOC5573930, LOD = 13.3) had the highest association with kdr. In chromosome 3, the highest LODs occurred in the two nonsynonymous mutations on VGSC (S679T and V410L), with LODs of 26.1 and 25.7, respectively. These mutations were followed by two nonsynonymous mutations in a membrane-bound transcription factor site-2 protease (LOC5580236), with LODs of 24.02 and 17.5, respectively. Additional target site associations occurring in this comparison were located in GPRGBB3 (LOC5569525, LODs = 7.13 and 7.62) and ACE (LOC5570776, LOD = 5.3). Table 2 shows the seven genes with nonsynonymous mutations, with LODs ranging from 3.97 to 7.93.

Discussion
We investigated the genomic differences between the three phenotypes discriminated after exposure to permethrin, including kdr (14.5% of total), recovery (32.4% of total), and dead (53.1% of total). We used a permethrin concentration recommended by the CDC bottle bioassay, in which knocked-down mosquitoes are recorded after 30 to 60 min from the exposure, and a mortality rate is calculated. Although, the methodology does not recommend further observations, recovery rates of 20-60% of the original knocked-down mosquitoes are commonly observed within 4 h of exposure [15,19], suggesting that recovery might be an important phenomenon in overall survival. The biological significance of recovery in the field is not well understood. A common assumption is that knocked-down mosquitoes are likely to die due to desiccation and predation. However, under ideal environmental conditions, individuals with the recovery mechanism could be favored by sublethal exposure to insecticides provided by the poor penetration of space sprays into sheltered indoor microhabitats where the mosquito rests [20]. In this study, we assume that recovery likely contributes to overall insecticide resistance in the field. Understanding the genomic differences between mosquitoes that exhibit kdr from those that recover or die after knockdown will pinpoint possible mechanisms regulating these phenotypes.
Our genomic analysis resulted in the identification of thousands of significant SNPs associated with the phenotypes. We focused on nonsynonymous mutations under the assumption that these can confer protein changes that affect the phenotype. Our results showed that SNPs at VGSC had the highest association with kdr but not with recovery. Significant but moderate association occurred between recovery and SNPs in additional insecticide target site, detoxification, or cuticle genes.
The role of VGSC as the primary site of action for the neurotoxic effects of pyrethroids in mammals and insects has been demonstrated. However, the actions of pyrethroids on secondary targets also have been associated with toxicity, mostly with the choreoathetosis and salivation (CS) intoxication syndrome by pyrethroids type 2 [21,22]. In our study, three nonsynonymous mutations at VGSC (V410L and S679T) were highly associated with the kdr phenotype but not with recovery. Additionally, nonsynonymous SNPs uniquely associated with kdr were located at the ACE/hydrolase (target site of organophosphate and carbamates), the G protein-activated inward rectifier potassium channel 3 (target site of novel insecticides), and the gamma-aminobutyric acid type B receptor subunit 2 (possible target site of cyclodienes). The level of association of these SNPs with the phenotype was much lower than VGSC but still significant.
The association between nonsynonymous mutations at VGSC and different levels of pyrethroid resistance has been confirmed in Ae. aegypti [8,9]. Moreover, some mutations interact and are restricted to geographical distributions [4]. In this study, two nonsynonymous mutations were identified in the VGSC, including V410L and S679T (S723T following M. domestica). V410L has been associated with kdr in Ae. aegypti from Mexico and has evolved in close linkage disequilibrium with V1016I [7]. Interestingly, V410L reduces the binding of pyrethroids to VGSC in mosquitoes by itself [8], whereas V1016I does not [9,10]. Individual genotyping at V410L in the mosquitoes used for our libraries showed strong positive association between kdr with resistant homozygous genotypes, dead with wild-type homozygous genotypes, and recovered with heterozygous genotypes. Moreover, the resistance allele frequencies (q) calculated from individual genotyping for kdr (q = 0.95), recovered (q = 0.48), and dead (q = 0.18) did not differ significantly from the genome sequencing frequencies obtained by the pooled libraries for kdr (q = 0.93), recovered (q = 0.43), and dead (q = 0.15). These results and previous observations on kdr-mutations in Ae. aegypti from Mexico lead us to infer that V410L segregates as a recessive allele. For example, a QTL mapping confirmed the recessive nature of the V1016I to confer permethrin resistance in Ae. aegypti from Mexico [5]. In addition, the strong linkage disequilibrium between V1016I and V410L, in which 95% of the individuals collected in 2014-2016 had resistant genotypes at both loci [7] suggests that both V1016I and V410L segregate as recessive alleles.
Additional kdr-associated mutations at VGSC include V1016I, F1534C, and S679T. The latter has low fitness and can only survive when F1534C is present [15]. Although, F1534C alone confers seven-to 14-fold resistance to pyrethroids [6], the combination of F1534C with V1016I enhanced the levels of permethrin and deltamethrin resistance in electrophysiology assays [10]. In our study, we did not identify V1016I or F1534C in our libraries for different reasons. V1016I was inconsistent between the biological replicates and therefore was eliminated in the independence test. F1534C was eliminated because C1534 is approaching fixation in Tapachula, therefore, the pipeline did not identify polymorphisms at this locus. A third nonsynonymous mutation, S679T, was strongly associated with permethrin resistance in this study. This mutation was previously identified in a deltamethrin-resistant population from Merida, Mexico. This nonsynonymous mutation corresponds to residue V723 in M. domestica, and it is probably located at the linker IS6-IIS1 in the VGSC. The role of S679T in pyrethroid resistance is unknown [23].
SNPs associated with recovery were in the acetylcholine receptor alpha, GPRGBB3 and the VGCC. Interestingly, the P1661S mutation at VGCC had high association with recovery (LOD = 10.9). Although the direct binding of pyrethroids to these channels has not been demonstrated, different biochemical or electrophysiological effects in these channels indicate they might play a secondary role in the toxicity of pyrethroids [18,22]. For example, five type 2 pyrethroids and permethrin (type 1) were potent enhancers of both calcium uptake and glutamate neurotransmitter release in rat brain synaptosomes [22,24]. Although, biochemical and electrophysiological studies show direct effects of pyrethroids on calcium channel function in vitro, a causal connection between these effects and pyrethroid intoxication remains unclear [21].
GABA-receptors are a major target site of picrotoxinin, chlorinated cyclodienes, and phenylpyrazoles [21,25]. Although some evidence exists of pyrethroid and DDT effects on insect GABA responses [26], the relative low potency and incomplete stereospecificity of pyrethroids as GABA receptor antagonists in functional assays do not support a significant role in the production of pyrethroid intoxication [21]. Whether the presence of nonsynonymous mutations in the GPRGBB3 in permethrin-exposed survivors is a result of selection by pyrethroids or other classes of insecticides requires further research.
Metabolism of insecticides is a major mechanism of resistance in Ae. aegypti. A common model of metabolic resistance assumes that the survival of mosquitoes exposed to insecticides results from enhanced metabolism (or sequestration), preventing the insecticide from reaching its target site [27]. In addition to this barrier, survival will depend on target site mutations that prevent the insecticide from exerting its toxic effects. Following the assumptions of this model, our bioassay showed that 85% of the mosquitoes were knocked down by permethrin, suggesting that metabolism was not a major mechanism to prevent intoxication of the target site. Of the 15% resistant to knockdown (kdr), 90% were homozygous resistant for the V410L mutation. The identification of metabolic mechanisms that prevent the insecticide from reaching its target site-as the model assumes-would be found in heterozygous and wild-type homozygous individuals in the kdr group (8% of the kdr mosquitoes were heterozygotes; 0% were wild-type homozygotes). We did not use this group of mosquitoes for comparisons because small sample sizes prevented us from completing pools of 25 individuals. But future experiments that include this group could decipher this model of metabolism resistance.
By using the permethrin discriminating concentration of 15 μg, mosquito survival was mostly explained by recovery (32.4%) and then by kdr (14.5%) in the mosquito population from Tapachula. Whether these responses are dose-dependent need to be addressed in future studies. Interestingly, 80% of the recovered mosquitoes were heterozygotes for the V410L locus in VGSC. This result suggests that carrying a single V410L allele does not protect a mosquito against knockdown but favors recovery once the pyrethroid dissociates from the ion channel and is metabolized and excreted. Therefore, the rates of recovery would be conferred by metabolism and detoxification mechanisms occurring between 1 and 4 h after insecticide exposure. Approximately 41 nonsynonymous SNPs that differ in the kdr vs recovery comparison and 36 that differ in the recovery vs dead comparison were identified. Except for five SNPs that overlapped in both comparisons (CYP9J32, CCEjhe2o, CYP325U1, GPXH3, and CYP325R1), unique genes were associated with kdr and recovery, suggesting different genes might explain differences between these phenotypes. Genes associated with kdr included three esterases (CCEae3O, CCEaeB1, and CCEunk7O); seven CYP6; three CYP325; CYP4J; and four redox, two delta and one epsilon GSTs. Interestingly, CCEunk7O was previously associated with kdr following exposure to permethrin in a previous QTL mapping study in Ae. aegypti from Mexico [5]. Because CCEunk7O is close to VGSC in the AaegL5 genome assembly, this association may be the result of a genetic sweep, as was observed recently in a study where several genes close to VGSC were highly associated with deltamethrin resistance [23]. The same study also identified SNPs at CCEae30, CCEaeB1, CYP325, and CYP4J in association with deltamethrin resistance in Merida, Mexico [19]. Additionally, CYP4J and CYP325 genes have been upregulated in pyrethroid resistant Ae. aegypti from the United States, Mexico, Vietnam, and Thailand [11,[28][29][30]. Specifically, in Ae. aegypti from Mexico, the upregulation of CYP325G3, CYP4J13, CYP6NAE1, GSTE2, and other genes was selected for after one generation of artificial selection with permethrin [29].
Enhanced metabolism can result from two different mechanisms: gene overexpression or allele variants conferring higher catalytic properties to the enzymes coded by these genes. To date, most studies have focused on expression analysis between resistant and laboratory susceptible strains, but a few studies have identified allele variants associated with resistance. For example, in CYP6P9A and CYP6P9B from An. funestus, a single allele was associated with pyrethroid resistance [37]. In addition, allele variants have been associated with temephos resistance in Ae. aegypti from Brazil or Thailand [38]. Interestingly, one study compared the transcription and copy number in pyrethroid-resistant strains from different continents [39]. The study showed that detoxification genes differentially expressed in resistant populations were contained in genomic clusters affected by copy number variations associated with resistance. Positive correlation occurred in three CCEs (CCEae3A, CCEae4A, and CCEae6A), CYP9J21, CYP9J22, CYP6BB2, and CYP6P12 [39]. In this study, we did not test for expression analysis. However, previous studies evaluated the expression profiles in Mexican Ae. aegypti colonies artificially selected for permethrin resistance [29]. The study showed that resistance ratios increased between 60-and 165-fold among strains, but this increase in resistance was associated with only slight changes in CYP expression profiles (less than three-fold). Mostly, the highest differences seemed to be constitutive [29]. Whether SNPs found in specific comparisons are directly associated with more efficient metabolism of permethrin during the first hour (kdr) or following recovery after 4 h, requires further research.
Identifying the specific genes associated with resistance phenotypes is the first step in the development of genetic SNP markers of resistance. Our results show that several detoxification genes are associated with kdr and recovery. Some genetic markers might only be artifacts of genetic sweeps or might follow different evolutionary mechanisms of selection. The applicability of this study beyond southeastern Mexican populations requires further investigation. Two studies have used exome-wide sequencing to identify SNPs associated with pyrethroid resistance in Mexico. These included two sites located less than 10 km apart in the Yucatan Peninsula (Vergel and Viva Caucel) [19,33] and Tapachula (this study), which is~1000 km from Yucatan. Both studies have shown a strong association between kdr and mutations at VGSC. Selection of kdr-mutations across geographical locations can be explained by the uniform practices of insecticide application by the Mexican vector control programs. Also, by the high gene flow in southeastern Ae. aegypti populations recorded in a large-scale mitochondrial population genetics study [40] and a small-scale SNP study in the Yucatan Peninsula [41]. Currently, kdr-mutations are widespread in Ae. aegypti from Mexico and the United States. Whether the same SNPs at compensatory or detoxification genes associated with insecticide resistance become selected across Ae. aegypti populations in North America remains unknown. Recent genome-wide comparisons in Ae. aegypti populations from California show large differentiation between genetic clusters due to recent introductions and from multiple genetically diverged populations [42,43]. For example, Southern California were less genetically diverse that Northern California populations. This low diversity was likely a signature of bottlenecks caused by recent founder effects and/or vector control measures [43].
Recent studies have shown the importance of reduced insecticide penetration of the cuticle as a mechanism of resistance in mosquito vectors. In Anopheles gambiae, thickening of the cuticle was associated with reduced permeability to pyrethroids in resistant mosquitoes [44,45]. These studies have revealed that the basis of cuticular thickening is quantitative changes in the composition of the cuticle. Furthermore, the role of certain enzymes (CYP4G16, CYP4G17) in enhancing the biosynthesis of epicuticular hydrocarbons and the upregulation of cuticular genes in resistant strains has been demonstrated. So far, 293 cuticle proteins (CPR) have been characterized in Anopheles gambiae [46], and at least 300 are present in the recently annotated Ae. aegypti genome assembly. Our study identified nonsynonymous mutations in several CPR genes; however, three genes had high LODs associated with kdr (LOC5571160 and LOC5571167) and recovery (LOC5577598). Further studies are required to test how these qualitative changes in cuticle proteins are associated with pyrethroid resistance in Ae. aegypti field populations.

Conclusion
Pyrethroid resistance in Ae. aegypti threatens our ability to control arboviral diseases. Understanding the mechanisms of pyrethroid resistance and the interactions and evolution of that resistance will be necessary to develop diagnostic tools to support insecticide management strategies. Two phenotypes-kdr and recovery-are involved in pyrethroid survival. This study identified mutations at VGSC controlling kdr but not recovery. Additional target site SNPs in the VGCC and GABA receptor genes were associated with recovery. Additionally, some specific detoxification genes were uniquely associated with kdr or with recovery. Understanding the role of these genes in the metabolism of pyrethroid will increase our knowledge of the evolution of resistance mechanisms in the field.

Mosquito colony and bioassays
We used a field mosquito colony named Colinas, which was collected in 2017 from Tapachula, Chiapas, Mexico (N 14' 55" 43.6, W 92' 14" 58.8). Briefly, we collected larvae from patios of approximately 25 houses in this neighborhood. Approximately 1,000 larvae were transferred to the Insectary at Centro Regional de Investigaciones en Salud Publica. Emerged mosquitoes were identified as Ae. aegypti and bloodfed to produce the F 1 offspring. Mosquito F1-egg papers were shipped to Colorado State University, where we performed the bottle bioassay on emerged adult mosquitoes to establish the levels of permethrin resistance relative to the New Orleans susceptible reference strain. We exposed approximately 75 mosquitoes to five different concentrations of permethrin for 1 h. Then, mosquitoes were removed from the treated bottle and were kept on a holding cup to score the mortality at 24 hours. Permethrin concentrations ranged from 0.1-1.5 μg/bottle for New Orleans and 7-50 μg/bottle for Colinas. Knockdown and mortality was scored at 1 or 24 h of observation, respectively. Following a binomial regression model and using the IRMA quick calculator [47], we calculated the lethal concentration that kills 50% of the mosquitoes. The permethrin LC 50  For the genome-wide association study, we used a permethrin-discriminating concentration (LC50 = 15 μg) that allowed us to discriminate three different phenotypes. The interior walls of 10 Wheaton bottles (250 ml) were coated with 15 μg of permethrin (Sigma-Aldrich, St. Louis, Missouri) diluted in 1 ml acetone. Following the evaporation of acetone overnight, we performed the bioassay. Approximately 40 non-bloodfed female mosquitoes (3-4 days old) were aspirated into each permethrin-coated bottle. After the mosquitoes had been exposed for 1 h, we transferred active mosquitoes to a clean 1-quart cup for observation. The remaining knocked-down mosquitoes in the bottle were transferred to a second cup for observation. These cups were placed in an incubator for 4 h. At this time, we recorded our three phenotypes: 1) "knockdown-resistant" (kdr) refers to those mosquitoes still active after 1 h of insecticide exposure and that maintained activity 4 h after being transferred to clean cups; 2) "recovered" refers to those mosquitoes that were knocked-down at 1 h and became active again in the 4 h after being transferred to clean cups; and 3) "dead" those mosquitoes that were knocked down at 1 h and did not become active again in the 4 h after being transferred to clean cups (Fig 1). In this study, we compared the genome of pooled mosquitoes from the three phenotypic groups: kdr vs recovered, recovered vs dead, and kdr vs dead.

Sample pooling and library preparation
We constructed six genomic libraries (gDNA): two for the knockdown-resistant (kdr) group, two for the recovery group, and two for the dead group. Each library consisted of pools of 25 mosquitoes. Before the individual mosquitoes were pooled, gDNA from each mosquito was quantified using the Quant-IT Pico Green kit (Life Technologies, Thermo Fisher Scientific Inc.). Approximately 40 ng of each individual DNA sample was used for a final DNA pool of 1 μg. Pooled DNA was sheared and fragmented by sonication to obtain fragments between 300-500 bp (Covaris Ltd., Brighton, UK). We prepared one library for each of the six DNA pools following the Low Sample (LS) protocol from the Illumina TrueSeqDNA PCR-Free Sample preparation guide (Illumina, San Diego CA). Because 47% of the Ae. aegypti genome consists of transposable elements and other forms of repetitive DNA [48], we performed an exome-capture hybridization to enrich for coding sequences using custom SeqCap EZ Developer probes (NimbleGen, Roche). Probes covered protein coding sequences (not including UTRs) in the AaegL1.3 genebuild using the exonic coordinates specified previously [49]. In total, 26.7 Mb of the genome (2%) was targeted for enrichment. TruSeq libraries were hybridized to the probes using the xGen Lock Down recommendations (Integrated DNA Technologies). The targeted DNA was eluted and amplified (10-15 cycles). Pair-end sequencing was run in a flow cell of NovaSEQ S4 and performed by the Genomics Core University of Colorado Anschutz Medical Campus (Aurora, Colo.).
allele counts in each of the four libraries (e.g., kdr1, kdr2, dead1, dead2). Allele frequencies for the alternant allele were calculated for each phenotype, and a goodness of fit test identified SNPs with consistent proportions within replicates (p > 0.05). Then, we built contingency tables and calculated the heterogeneity χ 2 with n-1 degrees of freedom to compare the proportion of the alternate allele between the phenotypes (the probability derived from this analysis was -log 10 transformed to provide a "LOD" value). Additionally, we calculated the expected heterozygosity (H exp ) of each site where Hexp ¼ 1 À P n i¼1 pi 2 and n is the number of alternate nucleotides at a site. We applied a Benjamini-Hochberg correction for false discovery rate [18] for each chromosome separately (α = 0.01). SNPs with a LOD above the cutoff were considered significant. Each significant SNP was annotated with the position (intergenic region, gene, introns, 3'-UTR, 5'-UTR, synonymous or nonsynonymous site). This information is included in S4-S6 Tables. Finally, a Fortran Program was used to identify whether the alternate SNP conferred a nonsynonymous mutation.
Individual genotyping of the V410L at VGSC was performed in the 50 individuals included in the libraries by using melting curve analysis of the allele-specific PCR methodology described in Saavedra-Rodriguez et al. 2018 [18].
Supporting information S1 Table. List of nonsynonymous SNPs differing between kdr and recovered Ae. aegypti mosquitoes exposed to permethrin. The gene ID, site (based in AaegL5 genome assembly), alternate nucleotide frequency in kdr, alternate nucleotide frequency in recovered, LOD (log 10 (p value)), heterozygosity in kdr, heterozygosity in recovered, total heterozygosity, amino acid substitution, functional category, vector base aliases and gene description. (XLSX) S2 Table. List of nonsynonymous SNPs differing between recovered and dead Ae. aegypti mosquitoes exposed to permethrin. The gene ID, site (based in AaegL5 genome assembly), alternate nucleotide frequency in recovered, alternate nucleotide frequency in dead, LOD (log 10 (p value)), heterozygosity in recovered, heterozygosity in dead, total heterozygosity, amino acid substitution, functional category, vector base aliases and gene description. (XLSX) S3 Table. List of nonsynonymous SNPs differing between kdr and dead Ae. aegypti mosquitoes exposed to permethrin. The gene ID, site (based in AaegL5 genome assembly), alternate nucleotide frequency in kdr, alternate nucleotide frequency in dead, LOD (-log 10 (p value)), heterozygosity in kdr, heterozygosity in dead, total heterozygosity, amino acid substitution, functional category, vector base aliases and gene description. (XLSX) S4 Table. List of SNPs differing between kdr and recovered Ae. aegypti mosquitoes exposed to permethrin. The gene ID, site (based in AaegL5 genome assembly), alternate nucleotide frequency in kdr, alternate nucleotide frequency in recovered, LOD (-log 10 (p value)), heterozygosity in kdr, heterozygosity in recovered, total heterozygosity, reference nucleotide, alternate nucleotide, gene orientation in chromosome, SNP position in gene, amino acid substitution classification (synonymous vs replacement), alternate nucleotide position in codon, amino acid residue, reference codon, reference amino acid, alternate codon, alternate aminoacid, chromosome, vector base identification, gene description and insecticide resistance category (target, cuticle, detoxification, other). (CSV) S5 Table. List of SNPs differing between recovered and dead Ae. aegypti mosquitoes exposed to permethrin. The gene ID, site (based in AaegL5 genome assembly), alternate nucleotide frequency in recovered, alternate nucleotide frequency in dead, LOD (-log 10 (p value)), heterozygosity in recovered, heterozygosity in dead, total heterozygosity, reference nucleotide, alternate nucleotide, gene orientation in chromosome, SNP position in gene, amino acid substitution classification (synonymous vs replacement), alternate nucleotide position in codon, amino acid residue, reference codon, reference amino acid, alternate codon, alternate aminoacid, chromosome, vector base identification, gene description and insecticide resistance category (target, cuticle, detoxification, other). (CSV) S6 Table. List of SNPs differing between kdr and dead Ae. aegypti mosquitoes exposed to permethrin. The gene ID, site (based in AaegL5 genome assembly), alternate nucleotide frequency in kdr, alternate nucleotide frequency in dead, LOD (-log 10 (p value)), heterozygosity in kdr, heterozygosity in dead, total heterozygosity, reference nucleotide, alternate nucleotide, gene orientation in chromosome, SNP position in gene, amino acid substitution classification (synonymous vs replacement), alternate nucleotide position in codon, amino acid residue, reference codon, reference amino acid, alternate codon, alternate amino acid, chromosome, vector base identification, gene description and insecticide resistance category (target, cuticle, detoxification, other). (CSV)