Exome‐wide association of deltamethrin resistance in Aedes aegypti from Mexico

Abstract Aedes aegypti is the major vector of a number of arboviruses that cause disease in humans. Without vaccines or pharmaceuticals, pyrethroid insecticides remain the major tool for public health protection. Pyrethroid resistance is now widespread. Replacement substitutions in the voltage‐gated sodium channel (vgsc) that reduce the stability of pyrethroid binding account for most of the resistance, but metabolic mechanisms also inactivate pyrethroids. High‐throughput sequencing and the A. aegypti L5 annotated physical map has allowed interrogation of the exome for genes and single‐nucleotide polymorphisms associated with pyrethroid resistance. We exposed females of A. aegypti from Mexico to a deltamethrin discriminating dose to designate them as resistant (active after 1 h) or susceptible (knocked down with no recovery after 4 h). The vgsc on chromosome 3 had the highest association, followed by genes proximal to vgsc. We identified potential detoxification genes located singly (eg HPX8C) or within clusters in chromosome 2 [three esterase clusters, two of cytochrome P450 monooxygenases (CYP)] and chromosome 3 (one cluster of 16 CYP325 and seven CYP9 genes). Deltamethrin resistance in A. aegypti is associated with mutations in the vgsc gene and a large assortment of genes.


Introduction
Insecticide resistance in Aedes aegypti mosquito populations threatens our ability to control arboviral diseases such as dengue, chikungunya, and Zika. During the last decade, intense global reliance on pyrethroid insecticides in vector control programmes worldwide has resulted in widespread pyrethroid resistance in A. aegypti populations (da-Cunha et al. , 2005;Harris et al. , 2010;Marcombe et al. , 2012;Aponte et al. , 2013;Flores et al. , 2013). Assessment and characterization of pyrethroid resistance mechanisms are important to the establishment of insecticide resistance management plans as insecticidal spraying efforts increase.
Various mechanisms of resistance have been reported in insects, including the insensitivity of the target site, increased metabolism, reduced penetration throughout the cuticle and changes in behaviour that usually result in avoidance of insecticide exposure. Although insecticide resistance is multifactorial, most research has focused on target site insensitivity and metabolic mechanisms separately, whereas reduced penetration and behavioural changes are just starting to be understood (Zalucki and Furlong, 2017;Balabanidou et al. , 2018).
Pyrethroids are neurotoxins that target the voltagegated sodium channel (vgsc ) in insect neurons, resulting in a disrupted action potential that culminates with rapid knockdown and mortality. Specific amino acid substitutions at vgsc result in target insensitivity in resistant insects. Around 57 nonsynonymous mutations in the vgsc have been associated with channel insensitivity (Rinkevich et al. , 2013); interestingly, convergence of knockdown resistance (kdr) mutations across different insect genera validate the important role of this mechanism. In A. aegypti , 11 nonsynonymous substitutions from different areas of the world have been reported (Rinkevich © Saavedra-Rodriguez et al. , 2018). Currently, five have been confirmed to contribute to insecticide resistance by electrophysiology Hirata et al. , 2014;Haddi et al. , 2017). Interestingly, co-existence of kdr mutations at different loci is common in insect species; furthermore, in A. aegypti these kdr mutations seem to be geographically distributed, conferring different levels of resistance (Hirata et al. , 2014;Saavedra-Rodriguez et al. , 2018). The three loci examined in this paper are 410, 1016 and 1534 corresponding to the position of amino acid substitutions. A valine at locus 410 (V410) confers susceptibility, whereas leucine (L410) confers resistance. A valine at locus 1016 (V1016) confers susceptibility, whereas isoleucine (I1016) confers resistance. A phenylalanine at locus 1534 (F1534) confers susceptibility, whereas cysteine (C1534) confers resistance.
Increased metabolism can result from overexpression (by amplification or cis-trans regulation) of insecticide metabolizing enzymes or by enzyme variants that provide higher or insecticide-specific catalytic properties (Hemingway and Karunaratne, 1998). Four enzyme families are associated with insecticide metabolism: cytochrome P450 monooxygenases (CYP), carboxyl/choline esterases (CCE), glutathione-S -transferases (GSTs) and reduction/oxidation (Red/ox) enzymes. Identification of detoxification genes in A. aegypti confirm the high redundancy of genes potentially involved in insecticide metabolism, describing 160 CYP, 44 CCE and 26 GSTs (Strode et al. , 2008). With improved microarray technology, whole transcriptomes had been interrogated to identify candidate genes conferring insecticide resistance in field or in artificially selected A. aegypti strains Marcombe et al. , 2009a;Bingham et al. , 2011;Bariami et al. , 2012;Grisales et al. , 2013;Kasai et al. , 2014). Unsurprisingly, genes and their expression patterns varied among regions of the globe. Nonetheless, a few common genes were successfully confirmed to produce insecticide-metabolizing proteins (Stevenson et al. , 2012).
As next-generation sequencing becomes more affordable, whole genomes or transcriptomes are interrogated for association with specific phenotypic traits. Few studies using this technology in mosquitoes have assessed the transcriptional changes and replacements identified in permethrin-resistant A. aegypti (David et al. , 2014;Faucon et al. , 2015). A recent complete partially annotated physical map of the A. aegypti genome (Matthews et al. , 2018) will allow us for the first time to perform genome-wide association studies. In the present study, we use an exome-enrichment genomic DNA sequencing approach (Juneja et al. , 2015) to identify singlenucleotide polymorphisms (SNPs) associated with deltamethrin resistance. We analysed an A. aegypti collection from Viva Caucel from Yucatán State in southern Mexico that was subjected to permethrin selection in the field from 1998 to 2010. Our hypothesis was that replacements in detoxification genes that have significantly different frequencies between resistant and susceptible mosquitoes are indicative of their potential involvement in deltamethrin resistance or are due to their physical proximity to vgsc on chromosome 3. The expected vgsc replacement polymorphisms were identified in the dataset and also allowed us to identify genome regions surrounding vgsc that are likely to have been involved in a genetic sweep.

Resistance phenotype
Deltamethrin resistance was phenotyped by exposing 390 A. aegypti adults from a field resistant strain collected in Viva Caucel to bottles coated with 3.0 µg of deltamethrin during 1 h. This concentration discriminated three phenotypes: (1) 149 mosquitoes were 'knockdown resistant' (38%), (2) 148 were 'susceptible' (38%) and (3) 93 'recovered' (24%) (Fig. 1). In this study, we compared pooled mosquitoes from the 'knockdown resistant' and 'susceptible' groups ( Fig.  1). We excluded the 'recovered' group to remove the uncertainty of assigning the incorrect phenotype. Also, we assumed that more contrasting phenotypes would allow us to get more contrasting SNPs underlying deltamethrin resistance.
A Benjamini-Hochberg correction (Benjamini and Hochberg, 1995) for false discovery rate (FDR) was applied to chromosomes 1, 2 and 3 separately (α = 0.01). A total of 14 532 SNPs (0.206% of the targeted exome) surpassed the FDR thresholds of 3.39, 3.52 and 3.22 for chromosomes 1, 2 and 3 respectively. These SNPs occurred in the 3′-untranslated region (UTR; n = 321), 5′-UTR (n = 443), coding sequence (n = 10 359; of which 2614 were nonsynonymous and 7745 were synonymous substitutions), intergenic regions (n = 40), noncoding RNA (n = 118) and introns (n = 3251). Fig. 2 shows the distribution of SNPs with −log 10 (prob) values above the FDR along chromosomes 1, 2 and 3. We provide a complete list of significant SNPs in Tables S2, S3 and S4. The approach used in this paper assumes that significant SNPs are independently distributed (in linkage equilibrium). This assumption is invalid, but we have no understanding of linkage disequilibrium blocks or haplotype structure in this or other mosquito species. Instead, we focus our results and discussion on those genes with the largest −log10(prob) values.
Our results showed that V410L located at IIS6 in vgsc had a high association [−log 10 (prob) = 17] with deltamethrin resistance. This mutation was recently associated with pyrethroid resistance in A. aegypti (Haddi et al. , 2017;Saavedra-Rodriguez et al. , 2018). To determine whether the library sequences accurately estimated the population of V410L allele frequencies, we performed a melting curve PCR SNP assay to genotype the individuals used to construct the four libraries (Saavedra-Rodriguez et al. , 2018). Fig. 3 shows the comparison of V410L allele frequencies obtained by library sequencing and by melting curve PCR SNP genotyping. Allele frequencies were not significantly different between replicates in any of the comparisons. Individual genotyping shows that the frequency of the resistant allele L410 was 0.88 and 0.94 in the knockdown-resistant individuals used in the replicates and 0.28 for the susceptible individuals used in the replicates (Fig. 3). In the sequencing libraries, L410 had a −log 10 (prob) = 17; frequencies in the knockdown-resistant libraries for replicates 1 and 2 were 0.88 and 0.73 respectively. In the susceptible individuals, frequencies were 0.10 and 0.3 in replicates 1 and 2 respectively. Heterogeneity χ 2 tests showed that the allele frequencies were similar regardless of the detection technique used, suggesting Testing for a genetic sweep by vgsc Because 45% of the 168 SNPs with high −log 10 (prob) > 10 were located in close proximity to vgsc (Fig. 2C), we analysed a 1.4 Mbp region upstream/downstream of this gene. By observing SNPs with −log 10 (prob) > 10 in this region, the sweep extended from a TATA-box binding protein upstream to vgsc to a vanin-2 protein downstream to vgsc . Genes in this sweep included a ubiquitin protease (OTU domain-containing protein 7B), a cytoskeleton binding protein (moesin/ezrin/radixin), a protein that forms a structural framework for protein-protein interactions (leucine-rich repeat-containing protein 47), a ribosomal biogenesis protein (large subunit GTPase) and two vanins involved in nitrogen compound metabolic processes (Fig. 4). Since none of these protein-coding genes have been functionally or genetically associated with insecticide resistance, we tested the possibility of a genetic sweep by rapid positive selection at vgsc rather than being selected due to their own activity as resistance loci. Fig. 4 shows the expected heterozygosity H exp of SNPs at vgsc and surrounding genes in knockdown-resistant and susceptible mosquitoes, where H exp = 1 − ∑ n i=1 p i 2 and n is the number of alternate nucleotides at a site. Heterozygosity at position L410 was low in knockdown-resistant mosquitoes and high in knockdown-susceptible mosquitoes. This is expected to occur when beneficial mutations (at vgsc ) are rare; but once a beneficial mutation has occurred it increases in frequency rapidly, thereby drastically reducing genetic variation in the population. To determine whether this pattern was genome wide, we tested for this pattern throughout all three chromosomes. Table 1 and Fig. S1 display the result of three pairwise t -tests, one for each chromosome testing whether the mean of the difference between heterozygosities in resistant and susceptible mosquitoes is different from zero. For each pairwise comparison of SNPs with −log 10 (prob) greater than the FDR we subtracted the expected heterozygosity of susceptible mosquitoes, Het (dead), from the heterozygosity of resistant mosquitoes, Het (alive). When Het (alive) > Het (dead) then the difference is positive, but when Het (dead) > Het (alive) the difference is negative. We would expect genes in surviving mosquitoes to have been swept and, therefore, have a lower heterozygosity. Conversely, we would expect genes in dead mosquitoes (homozygous susceptible) not to have been swept or swept at a lower rate (heterozygotes) and, therefore, have a higher heterozygosity. Along chromosomes 1 and 3 the mean difference was negative and significant in the paired t -test (Table 1), consistent with genetic sweeps.
However, along chromosome 2 the mean difference was positive and significant in the paired t -test (Table 1). Most importantly, in all three chromosomes there are distinct peaks where Het (alive) > Het (dead). This suggests that genes along all three chromosomes have not been swept and that the pattern in Fig. 4 is not genome wide.

Insecticide metabolism: detoxification genes
We identified 268 insecticide detoxification genes  in the AaeL5 assembly map and questioned which significant SNPs in our libraries were included in this categorical group. Genes with insecticide detoxification function surpassed the FDR cut-off calculated for each chromosome. Fig. 5 shows the −log 10 (prob) calculated for the SNPs comprised in 2, 57 and 46 detoxification genes dispersed across chromosomes 1, 2 and 3 respectively. Overall, 62% of the SNPs are classify as CYP genes, 20% as CCE, 16% as Red/ox and 4% as GSTs. In chromosome 1, only 13 SNPs were significant, with 11 SNPs in HPX8C and two SNPs in an oxidase peroxidase (HPX6) (Fig. 5A). Most SNPs were annotated as synonymous substitutions, except for an S766A in HPX8C [−log 10 (prob) = 3.76]. Both genes are currently annotated as chorion peroxidases. Several peroxidase isozymes have been described in Drosophila melanogaster as cyclooxygenase-like protein for both the actin filament bundle formation and the cortical actin strengthening during Drosophila oogenesis.
Because previous studies focused mostly on gene expression, we compared our list of genes with significant SNPs with 22 studies that had measured gene expression in mosquitoes exposed to sublethal doses of insecticides (Poupardin et al. , 2008;Riaz et al. , 2009;David et al. , 2010). We also included published articles that compared insecticide-resistant versus susceptible strains Marcombe et al. , 2009b;Saavedra-Rodriguez et al. , 2012;David et al. , 2014;Reid et al. , 2014;Estep et al. , 2017) and, most recently, allelic variants (David et al. , 2014;Dusfour et al. , 2015;Faucon et al. , 2015;. Fig. S3 shows the list of significant detoxification genes in our study and the expression that has been reported in previous studies. For example, in our study, CYP6N12 (LOC5571529) was significantly associated with deltamethrin resistance; this gene was significantly down/overexpressed in 12 studies, becoming a clear candidate for further study.
Furthermore, Campbell et al. (2019) analysed a compilation of high-throughput sequencing data from three pyrethroid treatment groups and interrogated them to identify functional categories that were overrepresented among all coding genes. Of eight functional categories, three were overrepresented, as calculated by a hypergeometric test, which tests the probability of a gene being represented more often than expected by chance. The three overrepresented groups were metabolism, transport and signalling. The detox gene subset would encompass both the metabolic category and the redox category, which was not overrepresented.

Discussion
With increased insecticide spraying to protect public health, the detection and characterization of resistance loci are a crucial requisite for adequate monitoring of the rise of insecticide resistance that could hamper vector control efforts. Currently, kdr is a major mechanism of pyrethroid resistance; around 11 amino acid replacements in vgsc have been associated with this trait in A. aegypti , and quantification of their frequencies in field populations is becoming a common practice to monitor pyrethroid resistance. Unfortunately, continued use of pyrethroids in the field is leading to the fixation of kdr mutations and possibly the selection of novel or compensatory mutations in vgsc and other genes associated with pyrethroid resistance. Identifying these novel or compensatory mutations is necessary to comprehensively monitor pyrethroid resistance in field populations. In this study, 52-80% of vgsc was sequenced across libraries (out of 7794 nucleotides covered by enrichment probes covering exons 5-30 and introns 4-29) and 19 SNPs were significantly associated with deltamethrin resistance. Notably, three nonsynonymous mutations in vgsc had significant association. Two are novel mutations, and their role in pyrethroid resistance is unknown. Contrastingly, V410L was highly associated with resistance in our study, bringing further evidence in its role in pyrethroid resistance (Haddi et al. , 2017). We recently tracked V410L in a population from Mexico collected in 2002, and frequencies have increased in the last decade in strong linkage disequilibrium to mutations I1016 and C1534 (Saavedra-Rodriguez et al. , 2018). Interestingly, I1016 and C1534 were not significantly associated with deltamethrin resistance in our study. To understand this observation, we individually genotyped mosquitoes included in our libraries and calculated allele frequencies of I1016 and C1534. I1016 had identical frequencies to those found in L410; however, coverage (>25 reads) was not reached at this site. The improvement of library preparation and DNA pool optimization will address the identification of this locus in association with pyrethroid resistance. Contrastingly, the frequency of C1534 was high (>0.9) in both susceptible and resistant libraries, so that −log 10 (prob) values were not high enough to surpass the chromosome FDR threshold. In Mexico, strong selection has pushed C1534 to fixation in many field populations (Vera-Maloof et al ., 2015). Therefore, the use of C1534 as a marker for pyrethroid resistance is no longer informative in Mexico, even though the role of C1534 in pyrethroid resistance has been validated by electrophysiology studies Haddi et al. , 2017). Furthermore, C1534 appears to be necessary for other mutations (I1016) to occur and increase in frequency in the population (Vera-Maloof et al ., 2015). Interestingly, we identified a possible genetic sweep driven by selection of vgsc nonsynonymous mutations associated with pyrethroid resistance. Similar genetic sweeps have been reported in D. melanogaster driven by selection of acetylcholinesterase mutations associated with carbamate and organophosphate resistance and by DDT resistance mediated by an insertion of a transposon in the 5′ regulatory region of CYP6g1 (Garud et al. , 2015). Similarly, strong signals of recent positive selection at several genes that are known to have a role in resistance have been identified in Anopheles , including vgsc , a cluster of GST genes including Gste2 (implicated in the metabolism of DDT and pyrethroids), and Cyp6p , a cluster of genes that include Cyp6p3 , which is upregulated in permethrin and bendiocarb-resistant mosquitoes (Anopheles gambiae 2017 Genomes Consortium, 2017). In A. aegypti , evidence of a genetic sweep (or hitchhiking effect) was demonstrated in a region of chromosome 1 that contains organophosphate-resistance-conferring genes (Yan et al. , 1998). The genetic sweep was reflected by low DNA polymorphisms and gene deletions for loci surrounding the EST-4 locus gene region. In contrast, a hitchhiking effect associated with DDT was not evident for the para (vgsc ) locus, in which heterozygosity was similar to other regions in the genome. This result was consistent with the hypothesis that, in the years since DDT was abandoned, the populations have had time to re-equilibrate. Luckily, with the A. aegypti physical map completed, we can test the gain/loss of genomic variability in populations exposed to different pressures of insecticides.
The role of metabolism in conferring pyrethroid resistance is unquestionable. Even if kdr mutations at vgsc protect mosquitoes from knockdown, pyrethroids are metabolized before reaching vgsc or after dissociation from vgsc by different enzymes to less toxic and more hydrophilic metabolites for excretion. The metabolism of pyrethroids occurs by ester hydrolysis (by CCE) and oxidation at methyl, methylene, alkenyl or aryl substituents (CYP), and ≥80 metabolites have been identified from cis and trans -permethrin (Casida et al. , 1983). To date, fewer than 10 genes have been functionally characterized and confirmed to metabolize pyrethroids in A. aegypti , including members of the CYP9 family (Stevenson et al. , 2012), CYP9M6 and CYP6BB2 , whereas other enzymes, such as CYP6Z8 and the Red/ ox (aldehyde dehydrogenases), catalyse the oxidation of permethrin intermediates (Chandor-Proust et al. , 2013;Lumjuan et al. , 2014). Noticeably, genes CYP6Z8, CYP9J26-28 and CYP9M6 belong to clusters of tandem genes. In our study, a large cluster of CYP325 genes at the p arm of chromosome 3 was strongly associated with deltamethrin. Although, the physiological function of the CYP325 in insects still remains unclear, some CYP325s (eg CYP325-E, K and G) are conserved among mosquito species, some only expanded in the Culex and Aedes genomes (CYP325X and Y) and some are species specific (Yan et al. , 2012). In our review, at least 13 CYP325 genes (Fig. S3) have been associated with changes of expression (two to threefold) in insecticide-resistant A. aegypti. In one study, CYP325N1 was 313-fold more expressed in a pyrethroid-resistant strain from Puerto Rico (Estep et al. , 2017). In other insects, only CYP325A3 was reported in a permethrin-resistant strain of Anopheles gambiae (David et al. , 2005); however, no clear orthologues were found in Aedes or Drosophila . Further research of the CYP325 cluster might help elucidate models for mechanisms of cluster expression and adaptation to specific insecticides, if any.
Similarly, identification and isolation of specific esterases that metabolize pyrethroids is ongoing, mainly because esterases have unspecific substrates that complicate protein isolation. In this study, several clustered CCE genes in chromosomes 2 (CCEaeA, CCEaeB and CCEaeD) and 3 (CCEunk and CCEglt) showed high association with deltamethrin resistance. CCEaeA family members are commonly reported in expression studies (Fig. S3). Genes from this cluster showed strong amplification in two deltamethrin-resistant strains from Thailand but were not present in two pyrethroid-  (Faucon et al. , 2015). Similarly, CCEae3A and 6A have been associated with temephos resistance in an A. aegypti strain for Thailand that was also resistant to pyrethroids . Although amplification of these CCE clusters are a possible mechanism to overexpress enzymes that metabolize insecticides, our study and others (Faucon et al. , 2015; have also identified allele variants within these genes. The roles of these nonsynonymous mutations in catalytic activity is of great interest. For example, amino acid substitutions increased diazinon hydrolysis by esterase E7 in Lucilia cuprina (Newcomb et al. , 1997) and M. domestica (Claudianos et al. , 1999), whereas a different substitution in the same gene was responsible for malathion resistance in L. cuprina (Campbell et al. , 1998). Characterization of CCE genes and their allelic variants will increase our knowledge of possible routes of metabolism.
We have mapped genomic regions associated with deltamethrin resistance in mosquitoes from Mexico. Our results validate the major role of vgsc in pyrethroid resistance, but fixation of replacement mutations in vgsc suggest that kdr is becoming well established in field populations. Additional mechanisms of deltamethrin resistance might play increasingly important roles in the future, following genomic changes in field populations that are routinely pressured with insecticides. It is therefore necessary to identify genes and possible mechanisms being selected by pyrethroids in general. So far, we have identified several detoxification genes that are highly variable in resistant mosquitoes. Research to further associate allelic variants with deltamethrin resistance is ongoing.

Mosquito samples, bioassays and DNA isolation
The A. aegypti population Viva Caucel from Yucatán State in southern Mexico (longitude −89.71827, latitude 20.99827) was collected in 2011 by Universidad Autónoma de Yucatán. This mosquito population was under field selection by a type-1 pyrethroid (permethrin) in vector control programmes from 1999 to 2010. We calculated the deltamethrin concentration to kill 50% of the mosquitoes (LC 50 ) by bottle bioassays to establish the level of resistance in the Viva Caucel population (F 3 ) relative to the susceptible New Orleans (NO). Around 20 non-blood-fed females (2-3 days old) of the Viva Caucel F 3 and NO were exposed to seven deltamethrin (Chem Service, West Chester, PA, USA) concentrations for 1 h, in three replicates. Mortality was scored following a recovery time of 24 h. The deltamethrin LC 50 at 24 h was 47.6-fold higher in the Viva Caucel population than in the NO susceptible strain (10.49 µg vs. 0.22 µg respectively) (Fig. S2). The NO strain was not used further in this study.
To discriminate the deltamethrin-resistant and susceptible individuals within the Viva Caucel population, 390 adult mosquitoes were exposed to 3 µg of deltamethrin for 1 h. This concentration is below the LC 50 calculated at 24 h post-exposure; however, 3 µg allowed us to discriminate the three phenotypes at 4 h of observation (instead of 24 h) while ensuring good DNA quality in 'susceptible' mosquitoes. Additionally, 4 h is a reasonable time for recovery, since longer periods increase the chances of predation and dehydration in knocked-down mosquitoes in the field.
The bioassay was performed by aspirating groups of 50 non-blood-fed female mosquitoes (3-4 days old) into deltamethrin-coated bottles for 1 h. After this time, active mosquitoes were transferred to 1 pint cardboard cups and placed into an incubator (28 °C and 70% humidity) for 4 h. These mosquitoes constituted the 'knockdown-resistant' group. Knocked-down mosquitoes were then transferred to a second cardboard cup. Four hours later, newly recovered mosquitoes were aspirated, frozen and labelled as 'recovered'. These were not used in the current study. The knocked-down mosquitoes that remained inactive at 4 h post-treatment were scored as 'susceptible'. DNA was isolated from individual mosquitoes that were categorized as either resistant or susceptible (Black and DuTeau, 1997) and resuspended in 150 µl of TE buffer (10 mM Tris hydrochloride, 1 mM ethylenediaminetetraacetic acid pH 8.0).

Sample pooling and library preparation
We constructed four genomic DNA libraries. A library was generated from the DNA of 25 knockdown-resistant females (Fig. 1), and this was replicated in a second set of 25 mosquitoes. The same was done with susceptible females (susceptible replicates 1 and 2) (Fig. 1). Before pooling, DNA from each individual mosquito was quantified using the Quant-IT Pico Green kit (Invitrogen, Eugene, OR, USA), and ~40 ng from 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 and 500 bp (Covaris Ltd, Brighton, UK). We prepared one library for each of the four DNA pools following the Low Sample (LS) protocol from the Illumina TrueSeqDNA PCR-Free Sample preparation guide (Illumina, San Diego, CA, USA).

Whole exome hybrid capture
Because of the large 1.38 Gbp size of the A. aegypti genome (Nene et al. , 2007), we used an exome capture method. To identify the proportion of L5 genes covered by the capture probes, the analogous L5 probe sequences were blasted (blastn) against the L5 reference transcriptome. The returned genes were cross-checked against the total number of L5 genes (18 081). Thus, 77% (13 941/18 081) of L5 genes were accounted for among the probe set. We performed an exome-capture hybridization to enrich for coding sequences using custom SeqCap EZ Developer probes (Roche NimbleGen Inc., Madison, WI, USA). Probes covered protein coding sequences (not including UTRs) in the AaegL1.3 genebuild using the exonic coordinates specified previously (Juneja et al. , 2015). In total, 26.7 Mb of the genome (2%) containing 15 735 genes was targeted for enrichment. TruSeq libraries were hybridized to the probes using the xGen®Lock®Down recommendations (Integrated DNA Technologies

Analysis pipeline
The raw sequence files (*.fastq) for each pair-ended gDNA library were aligned to a custom reference physical map generated from the assembly AaegL5 (Matthews et al. , 2018) using the package gsnap (Genomic Short-read Nucleotide Alignment; Wu and Nacu, 2010), which performed an SNP-tolerant alignment of the resulting (39-42) × 10 6 paired reads/library. Additional processing details for our in-house fortran pipeline are provided in Dickson et al. (2017). Briefly, we removed SNPs that were significantly different between replicates. Next, we grouped replicates that share the same sites together to build contingency tables and the heterogeneity χ 2 was calculated with (n − 1) nucleotides degrees of freedom. The probability derived from that analysis was −log 10 transformed to provide a −log 10 (prob) value. A Benjamini-Hochberg correction for FDR (Benjamini and Hochberg, 1995) was applied for chromosomes 1, 2 and 3 separately (α = 0.01) using sas.

Validation of nonsynonymous substitution V410L in vgsc
We determined the individual genotype of the mosquitoes used in our four libraries (n = 100) at two nonsynonymous substitutions in the vgsc . We performed melting curve analyses of allele-specific PCR products to genotype at locus V1016I (Saavedra-Rodriguez et al. , 2007) and a recently reported mutation at locus V410L (Haddi et al. , 2017) following an allele-specific PCR described in Saavedra-Rodriguez et al. (2018).

Supporting Information
Additional supporting information may be found in the online version of this article at the publisher's web site: Figure S1. Distribution of the differences in heterozygosity in knockdown resistant versus susceptible mosquitoes (H exp alive (HETA) H exp dead (HETD)) at each significant SNP in chromosomes 1, 2 and 3. Figure S2. Deltamethrin dose response curve for female Ae. aegypti from Viva Caucel and the reference strain from New Orleans. Exposure was performed using bottles coated with increasing concentrations of deltamethrin. Exposure time was 1 h and mortality was recorded at 24 h. Figure S3. List of significant detoxification genes in our study and the expression that has been reported in previous studies. Upregulated genes are highlighted in orange, downregulated are highlighted in blue, downregulated or upregulated in the same study are highlighted in grey. A number represents the number of polymorphisms found at a particular SNP in that gene. Table S1. Number of reads aligned AaeL5 across chromosome 1, 2 and 3 for each replicated libraries. The number of common reads between Table S2. List of Individual SNP, physical position of each SNP, gene description, mutation, frequencies of alternant nucleotide in resistant and susceptible groups, log 10 (prob), nucleotide, type of SNP, substitution, amino acid residue and codon position of each SNP for chromosome 1. Table S3. List of Individual SNP, physical position of each SNP, gene description, mutation, frequencies of alternant nucleotide in resistant and susceptible groups, log 10 (prob), nucleotide, type of SNP, substitution, amino acid residue and codon position of each SNP for chromosome 2. Table S4. List of Individual SNP, physical position of each SNP, gene description, mutation, frequencies of alternant nucleotide in resistant and susceptible groups, log 10 (prob), nucleotide, type of SNP, substitution, amino acid residue and codon position of each SNP for chromosome 3. Table S5. List of non synonymous mutations in deltamethrin resistant Ae. aegypti . The genes, description, SNPID, residue number, frequencies in resistant and susceptible and log 10 (prob) is shown for chromosomes 1, 2 and 3.