Skip to main content
Advertisement
  • Loading metrics

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

Abstract

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.

Author summary

Dengue, Zika, and chikungunya are viral diseases transmitted by the mosquito Aedes aegypti. Unfortunately, no vaccines or effective treatments are available, and public health strategies rely on the suppression of mosquito populations to reduce the impact of disease outbreaks. Pyrethroids are a common class of insecticides used to suppress mosquito populations. Unfortunately, mosquitoes have developed resistance to pyrethroids in several regions of the world, threatening disease control efforts. Since few classes of insecticides are available for public health, we have to find strategies that prevent or prolong the use of pyrethroids. Therefore, we need to understand mechanisms of resistance and develop methods to quantify these mechanisms in field populations before they become a threat to vector control strategies. Quantifying the frequency of mutations associated with target site resistance provides information about the presence of this mechanism in mosquito populations; however, we still need to understand additional mechanisms of resistance (e.g., metabolism, cuticle penetration) that play an important role in the survival of mosquitoes. In this study, we sequenced the genome of resistant (kdr), recovered, and dead mosquitoes after insecticide exposure. Our aim is to identify genomic variants associated with specific mechanisms of resistance. The results will improve the identification of specific genomic markers associated with resistance in the Aedes aegypti from Southeastern Mexico.

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 [14]. Heavy use of pyrethroid space sprays—due to their strong human safety profile—has created immense selection pressure for resistance [14]. 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 [68]. 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 F1 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).

thumbnail
Fig 1. Bioassay to differentiate three phenotypes in Aedes aegypti exposed to permethrin (15 ug/bottle) for 1 h.

Total number of mosquitoes used in bioassays are shown. Pooled libraries were prepared using 25 individual mosquitoes from each phenotypic group.

https://doi.org/10.1371/journal.pgen.1009606.g001

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 (-log10(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.

thumbnail
Table 1. Mean and standard error (SE) of LOD values (-log10(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.

https://doi.org/10.1371/journal.pgen.1009606.t001

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 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.

thumbnail
Fig 2. 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–log10(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

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.

Top nonsynonymous SNPs associated with kdr

In this section, we describe nonsynonymous mutations while assuming that these mutations confer changes in protein activity or functionality, therefore, making these proteins more likely to be subject to selection. S1 Table shows the list of nonsynonymous SNPs associated with kdr. On chromosome 1, 195 nonsynonymous mutations were located in 149 genes. The SNPs with the highest LODs were D765E (LOD = 22.15) and N772K (LOD = 10.84) in the ER degradation-enhancing alpha-mannosidase (LOC5566778). These were followed by P194S (LOD = 10.72) in the bis(5’-adenosyl)-triphosphatase ENPP4 (LOC5572000). In chromosome 2, 463 nonsynonymous mutations in 295 genes were identified (S1 Table). The SNP with the highest LOD was R577S (LOD = 13.34) in a zinc finger protein 883 (LOC5564538), followed by Q773P (LOD = 12.06) in an activating transcription factor 7-interacting protein (LOC5569161) and D875E (LOD = 11.79) in the fatty acid synthase (LOC5573930). In chromosome 3, 517 nonsynonymous mutations in 336 genes were significantly associated with kdr (S1 Table). The SNPs with the highest LODs were S679T (LOD = 13.7) and V408L (LOD = 12.25) at VGSC (LOC5567355). These SNPs correspond to loci V723 and V410L following the M. domestica annotation. Additional SNPs were H1711P (LOD = 12.08) and K1391R (LOD = 10.76), located in uncharacterized genes LOC5571908 and LOC5572722, respectively.

Insecticide target-site SNPs associated with kdr

We identified 192 significant SNPs in 23 genes coding for insecticide target sites; 70 SNPs were found only in VGSC. Table 2 shows the eight nonsynonymous SNPs identified in four genes, including acetylcholinesterase/hydrolase (ACE, LOC5570776), gamma-aminobutyric acid type B receptor subunit 2 (GPRGBB3, LOC5569525), G protein-activated inward rectifier potassium channel (KCNJ3, LOC5571228), and VGSC (LOC5567355). The highest LODs occurred in S679T and V410L in the VGSC (LOD = 13.7 and 12.25), followed by GPRGBB3 (LOD = 8.68 and 5.5) and ACE (LOD = 5.29 and 4.22).

thumbnail
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 = -log10(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).

https://doi.org/10.1371/journal.pgen.1009606.t002

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 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 wild-type 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).

thumbnail
Fig 3. 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

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).

thumbnail
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 = -log10(p value) and vector base identification (www.vectorbase.org).

https://doi.org/10.1371/journal.pgen.1009606.t003

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 4). In chromosome 1, the highest association occurred in CYP6AL3 (LOD = 7.57), GSTD3 (LOD = 7.27), and CYP4C51 (LOD = 5.18). In chromosome 2, the highest association occurred in CYP4AR2 (LOD = 9.4), CYP6N15 (LOD = 5.3), and aldo-keto reductase (LOD = 5.2). In chromosome 3, the highest association occurred in CYP285A (LOD = 5.29), CYP325U1 (LODs = 5.27, 5.22), and in the esterases CCEunk7O (LOD = 4.8) and CCEbe1O (LOD = 4.25).

thumbnail
Table 4. Nonsynonymous SNPs in detoxification genes associated with kdr, recovered, or dead Aedes aegypti exposed to permethrin.

The gene identification, site, amino acid residue (based in the AaegL5 genome assembly), frequency of alternative nucleotide for each phenotype, the associated LOD values, chromosome (Ch) and vector base identification (www.vectorbase.org) are shown. * indicates SNPs occurring in two or more pairwise comparisons.

https://doi.org/10.1371/journal.pgen.1009606.t004

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. Fig 4 shows the detoxification genes under directional and purifying selection. The frequency of the alternate allele 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.

thumbnail
Fig 4. Detoxification genes under purifying or directional selection in kdr, recovered, and dead Aedes aegypti exposed to permethrin.

We conducted three pairwise comparisons between phenotypes. The kdr vs recovery comparison involved two different resistant phenotypes. The recovery vs dead and, kdr vs dead comparisons involve one resistant and the susceptible phenotype (dead). If the frequency of the alternant SNP was higher in the resistant phenotype, we labeled the gene as directional selection. If the frequency of the alternant SNP was higher in the susceptible group, we labeled the SNP as purifying. The frequency was categorized as low (probably novel mutations) when the frequency ranged between 0 to 0.4, moderate from 0.4–0.8 and high 0.8 to 1.0. The detoxification genes with the highest LOD values for each comparison are highlighted in bold.

https://doi.org/10.1371/journal.pgen.1009606.g004

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 11,472 SNPs that differed significantly (2,713, 4,293 and 4,466 on chromosomes 1, 2, and 3, respectively).

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). The highest LOD was a N349D in the zinc finger protein 583-like gene (LOC110677088, LOD = 12.47). The following highest SNPs consisted of L44F in transmembrane protein 231 (LOC5574513, LOD = 12.45) and V5M in uncharacterized LOC23687531 (LOD = 12.06). In chromosome 3, 403 replacements in 279 genes were associated with recovery. The SNP with the highest association was a I145T in a zinc transporter ZIP1 (LOC5578569, LOD = 13.82), followed by a P4L in general transcription factor IIF subunit 1 (LOC5563781, LOD = 12.73) and S45P in a general odorant-binding protein 45-like (LOC5566895, LOD = 11.2).

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.

Detoxification SNPs associated with recovery

A total of 280 SNPs categorized in 103 detoxification-associated genes were associated with recovery. The 36 nonsynonymous SNPs were located in 22 genes, with LODs ranging from 3.3 to 6.6 (Table 4). In chromosome 1, the highest association occurred in CYP6P12v2 (LOC5576391, LOD = 6.6) and CYP6CC1 (LOC5565579, LOD = 6.4). In chromosome 2, the highest association with recovery occurred in CYP6M6 (LOC5571537, LOD = 6.34) and CYP6-A1 (LOC23687976, LOD = 4.9 and 4.6). In chromosome 3, the highest association occurred in CYP325U1 (LOC23687635, LOD = 5.9 and 5.74) and GPXH3 (LOC5578481, LOD = 5.72). Fig 4 shows that CYP6CC1, CYP6P12v2, CYP9J32, and CYP6M6 had SNPs under purifying selection, whereas CYP325U1, CCEglt2G, CYP4D24, GSTS-1, and GPXH3 were under moderate to high directional selection. Novel SNPs under directional selection occurred in CYP9J15 and CYP9J19.

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.

In the cuticle category, approximately 138 SNPs located across 49 genes were identified. The 13 nonsynonymous mutations were located in twelve genes (Table 3). The highest LOD occurred in the adult cuticle protein 1 (LOC5573913, LOD = 6.58), endocuticle SgAbd-6 (LOC5570308, LOD = 5.8) and cuticle protein CP14.6 (LOC5575766, LOD = 5.23).

In the detoxification category, 55 SNPs in 29 genes had nonsynonymous mutations with LODs ranging from 3.7 to 7.9. In chromosome 1, the highest association was CYP6CC1 (LOC5565579, LOD = 7.2, 7.1, and 6.9), HPX8C (LOC5564684, LOD = 7.0), and GSTD3 (LOC5568347, LOD = 6.5). In chromosome 2, the highest association was CYP6N6 (LOC5571528, LOD = 6.1), CYP4J13 (LOC5565336, LOD = 5.2), and CYP6N13 (LOC5571524, LOD = 5.1). In chromosome 3, the highest association was CCEunk7O (LOC5571034, LOD = 7.9), CCEglt4h (LOC5567220, LOD = 6.5) and a probable cytochrome LOC23687786 (LOD = 5.7). SNPs under high directional selection occurred in GSTD3, CYP4J13, CYP9-f2, CCEglt4H, and aldo-keto reductase. SNPs under purifying selection included CYP6CC1, CYP6BB2, CYP9J20, CYP6P12, and CYP325V1 (Fig 4).

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, 2830]. 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].

Specific genes associated with recovery included five CYP9J (-2, 15, 20, 26, and 29), seven CYP6 (including CYP6P12, CYP6BB2, CYP6Z9, and CYP6M6), and four delta and two theta GSTs. Similarly, the overexpression of several of these genes has been reported in insecticide resistant Ae. aegypti from Malaysia [31], Laos [32], Mexico [33] and Thailand [11, 29]. Moreover, the functional metabolism of pyrethroids has been demonstrated for CYP9J26, CYP9J28, CYP9J32, and CYP6BB2 in Ae. aegypti [34, 35] and CYP6P12 in Ae. albopictus [36].

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.

Methods

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 F1 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 LC50 for Colinas was 37.3-fold higher than that of the New Orleans susceptible strain. The LC50 values were 21.25 μg/bottle (95% CI = 18.5–24.4 μg) for Colinas and 0.56 μg/bottle (95% CI = 0.51–0.56 μg) for New Orleans.

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.).

Analysis pipeline

Our analysis compared the frequency of the alternate allele at each polymorphic site between the two different phenotypes. The alternate allele consisted of an allele not present in the reference AaegL5 genome assembly. In this study, we performed three pairwise comparisons: 1) kdr vs recovered, 2) recovered vs dead, and 3) kdr vs dead.

All sequencing data generated under this project is available at the National Center for Biotechnology Information (NCBI) Sequence Read Archive, Bioproject PRJNA731165. Colinas kdr replicates 1 and 2, are submitted as SRR14609309 and SRR14609308, respectively. Colinas recovered replicates 1 and 2, are SRR14609307 and SRR14609306, respectively. Colinas dead replicates 1 and 2, are SRR14609305 and SRR14609304.The raw sequence files (*.fastq) for each pair-ended gDNA library were aligned to a custom reference physical map generated from the assembly AaegL5 [50] using the package GSNAP (Genomic Short-read Nucleotide Alignment) [51]. We chose a minimum coverage of 25 and minimum base quality of 30.

A series of R scripts were used to split the genomic sites in the three chromosomes. A list of common sites within biological replicates and common sites between the phenotypes was then generated. The following script selected polymorphic sites (SNPs) and generated tables for 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 -log10 transformed to provide a “LOD” value). Additionally, we calculated the expected heterozygosity (Hexp) of each site where 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 S4S6 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 (–log10(p value)), heterozygosity in kdr, heterozygosity in recovered, total heterozygosity, amino acid substitution, functional category, vector base aliases and gene description.

https://doi.org/10.1371/journal.pgen.1009606.s001

(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 (–log10(p value)), heterozygosity in recovered, heterozygosity in dead, total heterozygosity, amino acid substitution, functional category, vector base aliases and gene description.

https://doi.org/10.1371/journal.pgen.1009606.s002

(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 (–log10(p value)), heterozygosity in kdr, heterozygosity in dead, total heterozygosity, amino acid substitution, functional category, vector base aliases and gene description.

https://doi.org/10.1371/journal.pgen.1009606.s003

(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 (–log10(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).

https://doi.org/10.1371/journal.pgen.1009606.s004

(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 (–log10(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).

https://doi.org/10.1371/journal.pgen.1009606.s005

(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 (–log10(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).

https://doi.org/10.1371/journal.pgen.1009606.s006

(CSV)

Acknowledgments

We thank Amaury Rodriguez-Penilla, Asha Gupta, and India Hilty for their help to process samples at Colorado State University.

References

  1. 1. Aponte A, Penilla P, Dzul-Manzanilla F, Che-Mendoza A, D. Lopez A, Solis F, et al. The pyrethroid resistance status and mechanisms in Aedes aegypti from the Guerrero state, Mexico. Pesticide Biochemistry and Physiology 2013;107:226–34.
  2. 2. Flores AE, Ponce G, Silva BG, Gutierrez SM, Bobadilla C, Lopez B, et al. Wide spread cross resistance to pyrethroids in Aedes aegypti (Diptera: Culicidae) from Veracruz state Mexico. J Econ Entomol. 2013;106(2):959–69. pmid:23786088
  3. 3. Marcombe S, Mathieu RB, Pocquet N, Riaz MA, Poupardin R, Selior S, et al. Insecticide resistance in the dengue vector Aedes aegypti from Martinique: distribution, mechanisms and relations with environmental factors. PLoS One. 2012;7(2):e30989. pmid:22363529
  4. 4. Moyes CL, Vontas J, Martins AJ, Ng LC, Koou SY, Dusfour I, et al. Contemporary status of insecticide resistance in the major Aedes vectors of arboviruses infecting humans. PLoS Negl Trop Dis. 2017;11(7):e0005625. pmid:28727779
  5. 5. Saavedra-Rodriguez K, Strode C, Flores Suarez A, Fernandez Salas I, Ranson H, Hemingway J, et al. Quantitative trait loci mapping of genome regions controlling permethrin resistance in the mosquito Aedes aegypti. Genetics. 2008;180(2):1137–52. pmid:18723882
  6. 6. Fan Y, Scott JG. The F1534C voltage-sensitive sodium channel mutation confers 7- to 16-fold resistance to pyrethroid insecticides in Aedes aegypti. Pest Manag Sci. 2020;76(6):2251–9. pmid:31981401
  7. 7. Saavedra-Rodriguez K, Maloof FV, Campbell CL, Garcia-Rejon J, Lenhart A, Penilla P, et al. Parallel evolution of vgsc mutations at domains IS6, IIS6 and IIIS6 in pyrethroid resistant Aedes aegypti from Mexico. Sci Rep. 2018;8(1):6747. pmid:29712956
  8. 8. Haddi K, Tome HVV, Du Y, Valbon WR, Nomura Y, Martins GF, et al. Detection of a new pyrethroid resistance mutation (V410L) in the sodium channel of Aedes aegypti: a potential challenge for mosquito control. Sci Rep. 2017;7:46549. pmid:28422157
  9. 9. Du Y, Nomura Y, Satar G, Hu Z, Nauen R, He SY, et al. Molecular evidence for dual pyrethroid-receptor sites on a mosquito sodium channel. Proc Natl Acad Sci U S A. 2013;110(29):11785–90. pmid:23821746
  10. 10. Chen M, Du Y, Wu S, Nomura Y, Zhu G, Zhorov BS, et al. Molecular evidence of sequential evolution of DDT- and pyrethroid-resistant sodium channel in Aedes aegypti. PLoS Negl Trop Dis. 2019;13(6):e0007432. pmid:31158225
  11. 11. Strode C, Wondji CS, David JP, Hawkes NJ, Lumjuan N, Nelson DR, et al. Genomic analysis of detoxification genes in the mosquito Aedes aegypti. Insect Biochem Mol Biol. 2008;38(1):113–23. pmid:18070670
  12. 12. Balabanidou V, Grigoraki L, Vontas J. Insect cuticle: a critical determinant of insecticide resistance. Curr Opin Insect Sci. 2018;27:68–74. pmid:30025637
  13. 13. Zalucki MP, Furlong MJ. Behavior as a mechanism of insecticide resistance: evaluation of the evidence. Current Opinion in Insect Science. 2017;21:19–25. pmid:28822484
  14. 14. Donnelly MJ, Corbel V, Weetman D, Wilding CS, Williamson MS, Black WCt. Does kdr genotype predict insecticide-resistance phenotype in mosquitoes? Trends Parasitol. 2009;25(5):213–9. pmid:19369117
  15. 15. Vera-Maloof FZ, Saavedra-Rodriguez K, Elizondo-Quiroga AE, Lozano-Fuentes S, Black Iv WC. Coevolution of the Ile1,016 and Cys1,534 Mutations in the Voltage Gated Sodium Channel Gene of Aedes aegypti in Mexico. PLoS Negl Trop Dis. 2015;9(12):e0004263. pmid:26658798
  16. 16. CDC. CDC Bottle Bioassay. (https://wwwcdcgov/parasites/education_training/lab/bottlebioassayhtml). 2013.
  17. 17. Kuri-Morales PA, Correa-Morales F, Gonzalez-Acosta C, Moreno-Garcia M, Santos-Luna R, Roman-Perez S, et al. Insecticide susceptibility status in Mexican populations of Stegomyia aegypti (= Aedes aegypti): a nationwide assessment. Med Vet Entomol. 2018;32(2):162–74. pmid:29165810
  18. 18. Benjamini Y, Hochberg Y. Controlling the False Discovery Rate—a Practical and Powerful Approach to Multiple Testing. Journal of the Royal Statistical Society Series B-Methodological. 1995;57(1):289–300.
  19. 19. Saavedra-Rodriguez K, Campbell CL, Lenhart A, Penilla P, Lozano-Fuentes S, Black WCt. Exome-wide association of deltamethrin resistance in Aedes aegypti from Mexico. Insect Mol Biol. 2019;28(5):591–604. pmid:30758862
  20. 20. Perich MJ, Davila G, Turner A, Garcia A, Nelson M. Behavior of resting Aedes aegypti (Culicidae: Diptera) and its relation to ultra-low volume adulticide efficacy in Panama City, Panama. J Med Entomol. 2000;37(4):541–6. pmid:10916294
  21. 21. Soderlund DM. State-Dependent Modification of Voltage-Gated Sodium Channels by Pyrethroids. Pestic Biochem Physiol. 2010;97(2):78–86. pmid:20652092
  22. 22. Breckenridge CB, Holden L, Sturgess N, Weiner M, Sheets L, Sargent D, et al. Evidence for a separate mechanism of toxicity for the Type I and the Type II pyrethroid insecticides. Neurotoxicology. 2009;30 Suppl 1:S17–31. pmid:19766671
  23. 23. Rinkevich FD, Du Y, Dong K. Diversity and Convergence of Sodium Channel Mutations Involved in Resistance to Pyrethroids. Pestic Biochem Physiol. 2013;106(3):93–100. pmid:24019556
  24. 24. Clark JM, Symington SB. Neurotoxic implications of the agonistic action of CS-syndrome pyrethroids on the N-type Ca(v)2.2 calcium channel. Pest Manag Sci. 2008;64(6):628–38. pmid:18383452
  25. 25. Casida JE, Durkin KA. Novel GABA receptor pesticide targets. Pestic Biochem Physiol. 2015;121:22–30. pmid:26047108
  26. 26. Lummis SC. GABA receptors in insects. Comp Biochem Physiol C. 1990;95(1):1–8. pmid:1971549
  27. 27. Hemingway J, Hawkes NJ, McCarroll L, Ranson H. The molecular basis of insecticide resistance in mosquitoes. Insect Biochem Mol Biol. 2004;34(7):653–65. pmid:15242706
  28. 28. Reid WR, Thornton A, Pridgeon JW, Becnel JJ, Tang F, Estep A, et al. Transcriptional analysis of four family 4 P450s in a Puerto Rico strain of Aedes aegypti (Diptera: Culicidae) compared with an Orlando strain and their possible functional roles in permethrin resistance. J Med Entomol. 2014;51(3):605–15. pmid:24897853
  29. 29. Saavedra-Rodriguez K, Suarez AF, Salas IF, Strode C, Ranson H, Hemingway J, et al. Transcription of detoxification genes after permethrin selection in the mosquito Aedes aegypti. Insect Mol Biol. 2012;21(1):61–77. pmid:22032702
  30. 30. Strode C, de Melo-Santos M, Magalhaes T, Araujo A, Ayres C. Expression profile of genes during resistance reversal in a temephos selected strain of the dengue vector, Aedes aegypti. PLoS One. 2012;7(8):e39439. pmid:22870187
  31. 31. Ishak IH, Kamgang B, Ibrahim SS, Riveron JM, Irving H, Wondji CS. Pyrethroid Resistance in Malaysian Populations of Dengue Vector Aedes aegypti Is Mediated by CYP9 Family of Cytochrome P450 Genes. PLoS Negl Trop Dis. 2017;11(1):e0005302. pmid:28114328
  32. 32. Marcombe S, Fustec B, Cattel J, Chonephetsarath S, Thammavong P, Phommavanh N, et al. Distribution of insecticide resistance and mechanisms involved in the arbovirus vector Aedes aegypti in Laos and implication for vector control. PLoS Negl Trop Dis. 2019;13(12):e0007852. pmid:31830027
  33. 33. Campbell CL, Saavedra-Rodriguez K, Kubik TD, Lenhart A, Lozano-Fuentes S, Black WCt. Vgsc-interacting proteins are genetically associated with pyrethroid resistance in Aedes aegypti. PLoS One. 2019;14(1):e0211497. pmid:30695054
  34. 34. Stevenson BJ, Pignatelli P, Nikou D, Paine MJ. Pinpointing P450s associated with pyrethroid metabolism in the dengue vector, Aedes aegypti: developing new tools to combat insecticide resistance. PLoS Negl Trop Dis. 2012;6(3):e1595. pmid:22479665
  35. 35. Kasai S, Komagata O, Itokawa K, Shono T, Ng LC, Kobayashi M, et al. Mechanisms of pyrethroid resistance in the dengue mosquito vector, Aedes aegypti: target site insensitivity, penetration, and metabolism. PLoS Negl Trop Dis. 2014;8(6):e2948. pmid:24945250
  36. 36. Ishak IH, Riveron JM, Ibrahim SS, Stott R, Longbottom J, Irving H, et al. The Cytochrome P450 gene CYP6P12 confers pyrethroid resistance in kdr-free Malaysian populations of the dengue vector Aedes albopictus. Sci Rep. 2016;6:24707. pmid:27094778
  37. 37. Riveron JM, Irving H, Ndula M, Barnes KG, Ibrahim SS, Paine MJ, et al. Directionally selected cytochrome P450 alleles are driving the spread of pyrethroid resistance in the major malaria vector Anopheles funestus. Proc Natl Acad Sci U S A. 2013;110(1):252–7. pmid:23248325
  38. 38. Helvecio E, Romao TP, de Carvalho-Leandro D, de Oliveira IF, Cavalcanti A, Reimer L, et al. Polymorphisms in GSTE2 is associated with temephos resistance in Aedes aegypti. Pestic Biochem Physiol. 2020;165:104464. pmid:32359546
  39. 39. Faucon F, Gaude T, Dusfour I, Navratil V, Corbel V, Juntarajumnong W, et al. In the hunt for genomic markers of metabolic resistance to pyrethroids in the mosquito Aedes aegypti: An integrated next-generation sequencing approach. PLoS Negl Trop Dis. 2017;11(4):e0005526. pmid:28379969
  40. 40. Gorrochotegui-Escalante N, Gomez-Machorro C, Lozano-Fuentes S, Fernandez-Salas L, De Lourdes Munoz M, Farfan-Ale JA, et al. Breeding structure of Aedes aegypti populations in Mexico varies by region. Am J Trop Med Hyg. 2002;66(2):213–22. pmid:12135296
  41. 41. Saavedra-Rodriguez K, Beaty M, Lozano-Fuentes S, Denham S, Garcia-Rejon J, Reyes-Solis G, et al. Local evolution of pyrethroid resistance offsets gene flow among Aedes aegypti collections in Yucatan State, Mexico. Am J Trop Med Hyg. 2015;92(1):201–9. pmid:25371186
  42. 42. Lee Y, Schmidt H, Collier TC, Conner WR, Hanemaaijer MJ, Slatkin M, et al. Genome-wide divergence among invasive populations of Aedes aegypti in California. BMC Genomics. 2019;20(1):204. pmid:30866822
  43. 43. Pless E, Gloria-Soria A, Evans BR, Kramer V, Bolling BG, Tabachnick WJ, et al. Multiple introductions of the dengue vector, Aedes aegypti, into California. PLoS Negl Trop Dis. 2017;11(8):e0005718. pmid:28796789
  44. 44. Balabanidou V, Kampouraki A, MacLean M, Blomquist GJ, Tittiger C, Juarez MP, et al. Cytochrome P450 associated with insecticide resistance catalyzes cuticular hydrocarbon production in Anopheles gambiae. Proc Natl Acad Sci U S A. 2016;113(33):9268–73. pmid:27439866
  45. 45. Yahouedo GA, Chandre F, Rossignol M, Ginibre C, Balabanidou V, Mendez NGA, et al. Contributions of cuticle permeability and enzyme detoxification to pyrethroid resistance in the major malaria vector Anopheles gambiae. Sci Rep. 2017;7(1):11091. pmid:28894186
  46. 46. Zhou Y, Badgett MJ, Orlando R, Willis JH. Proteomics reveals localization of cuticular proteins in Anopheles gambiae. Insect Biochem Mol Biol. 2019;104:91–105. pmid:30278207
  47. 47. Lozano-Fuentes S, Saavedra-Rodriguez K, Black WCt, Eisen L. QCal: a software application for the calculation of dose-response curves in insecticide resistance bioassays. J Am Mosq Control Assoc. 2012;28(1):59–61. pmid:22533088
  48. 48. Nene V, Wortman JR, Lawson D, Haas B, Kodira C, Tu ZJ, et al. Genome sequence of Aedes aegypti, a major arbovirus vector. Science. 2007;316(5832):1718–23. pmid:17510324
  49. 49. Juneja P, Ariani CV, Ho YS, Akorli J, Palmer WJ, Pain A, et al. Exome and transcriptome sequencing of Aedes aegypti identifies a locus that confers resistance to Brugia malayi and alters the immune response. PLoS Pathog. 2015;11(3):e1004765. pmid:25815506
  50. 50. Matthews BJ, Dudchenko O, Kingan S, Koren S, Antoshechkin I, Crawford JE, et al. Improved Aedes aegypti mosquito reference genome assembly enables biological discovery and vector control. Nature. 2018:(in press).
  51. 51. Wu TD, Nacu S. Fast and SNP-tolerant detection of complex variants and splicing in short reads. Bioinformatics. 2010;26(7):873–81. pmid:20147302