Insecticide-Driven Patterns of Genetic Variation in the Dengue Vector Aedes aegypti in Martinique Island

Effective vector control is currently challenged worldwide by the evolution of resistance to all classes of chemical insecticides in mosquitoes. In Martinique, populations of the dengue vector Aedes aegypti have been intensively treated with temephos and deltamethrin insecticides over the last fifty years, resulting in heterogeneous levels of resistance across the island. Resistance spreading depends on standing genetic variation, selection intensity and gene flow among populations. To determine gene flow intensity, we first investigated neutral patterns of genetic variability in sixteen populations representative of the many environments found in Martinique and experiencing various levels of insecticide pressure, using 6 microsatellites. Allelic richness was lower in populations resistant to deltamethrin, and consanguinity was higher in populations resistant to temephos, consistent with a negative effect of insecticide pressure on neutral genetic diversity. The global genetic differentiation was low, suggesting high gene flow among populations, but significant structure was found, with a pattern of isolation-by-distance at the global scale. Then, we investigated adaptive patterns of divergence in six out of the 16 populations using 319 single nucleotide polymorphisms (SNPs). Five SNP outliers displaying levels of genetic differentiation out of neutral expectations were detected, including the kdr-V1016I mutation in the voltage-gated sodium channel gene. Association tests revealed a total of seven SNPs associated with deltamethrin resistance. Six other SNPs were associated with temephos resistance, including two non-synonymous substitutions in an alkaline phosphatase and in a sulfotransferase respectively. Altogether, both neutral and adaptive patterns of genetic variation in mosquito populations appear to be largely driven by insecticide pressure in Martinique.


Introduction
During their life-span, mosquitoes face multiple selective pressures, including climatic constraints linked to an obligate aquatic larval stage, together with intense insecticide pressure. Insecticide pressure has led to the spreading of resistance alleles in all target species, and the managing of insecticide resistance in vectors of pathogens represents a major challenge for humanity. Originating from Africa, Aedes aegypti (Diptera, Culicidae) has spread throughout the tropics in historical times [1]. This mosquito species is the primary vector for dengue fever, yellow fever, and chikungunya virus in urban and sub-urban tropical areas across the world. In Martinique Island, a French overseas department located in the eastern Caribbean in the Lesser Antilles, the dengue virus (Flaviviridae) is transmitted by Ae. aegypti following an endemo-epidemic pattern and is now a major public health issue [2]. In the absence of vaccine or specific treatment against the disease, vector control by the use of chemical or biological agents against Ae. aegypti is currently the main method for reducing dengue transmission, in addition to environmental management and educational programs Over the last fifty years, larval and adult Ae. aegypti populations from Martinique have been intensively targeted with insecticide applications. Organochlorines such as DDT were first used since the 1950s and were replaced in the early 1990s by pyrethroid-based formulations (i.e. deltamethrin) for adult control by aerial spraying applications [3]. In addition to adult control, the organophosphate temephos (Abate®) was used for decades for larval control and there was a switch in 2009 to the bio-insecticide Bacillus thuringiensis var. israelensis (Bti). Because of the intense pesticide selection pressure in Martinique, vector control is now facing operational challenges with the emergence and development of insecticide resistance in Ae. aegypti populations [4] [5]. Resistance to organophosphates and pyrethroids has been reported since the 1980s and the 1990s respectively. A previous study has shown that deltamethrin and temephos resistances were widespread in Martinique. Molecular screening for common insecticide target-site mutations revealed the presence of the ''knock-down resistance'' V1016I kdr mutation at high frequency (>87%), and real time quantitative RT-qPCR showed the potential involvement of several candidate detoxification genes in insecticide resistance [3]. This study also showed that the insecticide resistance in Ae. aegypti populations was not homogeneously distributed across Martinique and that there were potential links existing between resistance distribution and vector control operations but also other factors such as climate and urbanization. Indeed, these factors are likely affecting the dynamics and the selection of insecticide resistance mechanisms in multiple ways [6]. The speed of resistance spreading depends on standing genetic variation in natural populations, selection intensity, local adaptation and gene flow among populations. The different climate and geological areas of Martinique with the northern part, mountainous, volcanic and rainy and the southern part slightly rugged, sunny and dry, together with the constant insecticide selection pressures and the well-developed road network that may favour passive gene flows, could lead to specific patterns of genetic diversity and population structure. The analyse of both neutral and adaptive genetic variation using genome-wide distributed molecular markers represents a powerful tool to better understand the genetic structure of these populations as well as the molecular mechanisms underlying insecticide resistance in Martinique Island.
First, using extensive field sampling (16 populations and 31 sub-populations) in different ecotypes of the island and seven microsatellite loci, we investigated neutral patterns of genetic variability and population structure at macro and microgeographic scales across Ae. aegypti populations from Martinique. Then, we investigated the association between insecticide resistance (temephos and deltamethrin) and 319 single nucleotide polymorphisms (SNPs) in 6 populations, in order to identify candidate genes for temephos and deltamethrin resistances in Ae. aegypti. We found that both neutral (microsatellites) and adaptive (SNPs) genetic variation is affected by insecticide pressure in Ae. aegypti populations in Martinique.

Sampling localities and bioassays
Two laboratory strains from French Polynesia (Bora-Bora) and Benin (SBE) and sixteen field-caught populations were used in the study. The two laboratory strains are susceptible to all insecticides and have been used as reference strains for resistance assays. Ae. aegypti was sampled from individual houses as larvae or pupae in sixteen localities of Martinique from January to March 2009 ( Figure 1, Table S1). Each population was constituted with about one thousands of larvae originated from 10 to 20 larval breeding sites (domestic breeding habitats) located in a 250 meters radius. Larvae and pupae were reared until the adult stage (F0 generation) and then morphologically identified. The adults were kept in a cage for copulation and females were blood-fed on rabbit (protocol CEEA-LR-13002 agreed by The French National Committee of ethics for animal experimentation) to obtain a F1 progeny. When a convenient number of F1 eggs were obtained, F0 were stored at -20°C until molecular analysis. Larvae and adults obtained from the F1 progeny were used for bioassays. The chosen populations covered most of the island ecotypes (coastal, mountainous, rural and urban) and two of them were sampled in two small islands (Ilet Long and Ilet Anonyme) located one kilometer away from the Atlantic coast ( Figure 1, Table S1). This sampling also covered a large range of insecticide pressures, from no treatment up to 81 interventions between 2006 and 2009 (Table S1). No specific permits were required for mosquito collection. We confirm that most of the sampling locations were not privately-owned. We asked the owners to sample in their back yards if needed. We also confirm that the locations were not protected in any way and that the field and laboratory studies did not involve endangered or protected species.
As described in [3], larval and adult bioassays were performed following WHO protocols [7] to measure the resistant ratios (RR 50 and RR 95 ) of twelve field-caught populations, the SBE strain (from Benin) and using the laboratory Bora-Bora strain as the susceptible reference strain. Larval and adult bioassays were carried out using temephos (Organophosphate) and deltamethrin (Pyrethroid), respectively.

Microsatellite genotyping and population genetic analyses
Microsatellite genotyping and population genetic analyses were implemented at the macro-geographic level using 16 F0 samples of 30 females and 10 males. In 3 localities, Fort-de France, St-Anne and Vauclin, individuals from 10, 4 and 4 distinct larval breeding sites respectively, were used for a micro-geographic level analysis. DNA was extracted from adult mosquitoes of both sexes using 2% CTAB as described in [8]. DNA pellets were resuspended in 20 µL sterile water and stored at -20 °C until analysis. Individual genotypes were scored at a set of seven microsatellite markers 34/72, AEDGA, AED19, AEDC, A10, 38/38 and H08. The microsatellites and PCR conditions have previously been detailed [9,10]. DNA was amplified in a GeneAmp 9600 thermal cycler (PerkinElmer, USA). PCR products were analyzed using ABI Prism™ 3100 (Applied Biosystems), and genotyped with GeneMapper (Applied Biosystems) (for details see 9). The microsatellite 38/38 was not polymorphic in our populations and was not kept for further analyses.
Genetic diversity by locus and sample and overall was characterized by estimating the number of alleles and allelic richness (Rs, [11]) using the software FSTAT2.9.3.2 [12]. Linkage disequilibrium, deviations from Hardy-Weinberg equilibrium (HWE) and genetic differentiation indices were assessed using GENEPOP4.7 software [13]. F IS and F ST estimates were calculated according to [14] and tested for statistical significance with exact tests available in GENEPOP4.7. The overall significance of multiple tests was estimated by Fisher's combined probability tests. Nominal significance levels for multiple testing were corrected using the sequential Bonferroni procedure [15]. Genetic isolation by geographic distance was assessed with Mantel's test in GENEPOP4.7, using the correlation between genetic and geographic distances by the regression of pairwise F ST /(1-F ST ) on the natural logarithm (Ln) of straight-line geographic distance [16]. Genetic structure was also studied at the macro and micro-geographic levels by examining the partition of neutral genetic variation using an AMOVA (Analysis of Molecular Variance, Arlequin 3.1 [17]). A bayesian clustering method implemented in STRUCTURE 2.3 [18] was used to assign individuals from all 16 populations to a pre-determined number of clusters (K=1-16) based on multi-locus microsatellite data. For each run, a burn-in period of 10,000 steps was followed by 30,000 iterations under the admixture model and the assumption of correlated allele frequencies among populations. To assist the clustering, sampling locations were used as prior information (LOCPRIOR model). For each K of 1-16, 20 runs were performed. Estimated log probabilities (Ln P(D)) were averaged across runs and compared to determine the posterior probability of each K. K with the highest likelihood was retained. Pearson's correlation between genetic diversity estimates (allelic richness and F IS ) and environmental variables (altitude, human density) and resistance ratio were tested using the "stat" package in the R software version 2.10 [19].

SNP genotyping
Six out of the 16 populations (GRMN, ILA, ILO, LAM, RSAL and SAN) representative of the range of insecticide resistance levels observed in Martinique were analysed for 319 SNPs. 318 of these SNPs were developed in the Ae. aegypti laboratory strain Bora-Bora using transcriptomic data in [20]. Briefly, short cDNA tags were obtained from pools of larvae and sequenced using a Genome Analyzer I (Illumina Corporation, USA). The reads were filtered, mapped to Ae. aegypti genome assembly AaegL 1.1 (see VectorBase, http:// aaegypti.vectorbase.org/) and 1520 SNPs with a coverage > 10x and a minimal allelic frequency of 5% were detected. The 318 SNPs used in this study were successfully transferred and genotyped in the natural Ae. aegypti populations of the Martinique Island [20]. They were distributed in 226 out of the 4,758 supercontigs in the reference genome. The knock down kdr-V1016I mutation in the voltage-gated sodium channel gene was also included in the SNP panel, as it was the only one found in a previous scan in Martinique populations among the four I1011M, I1011V, V1016I and V1016G kdr mutations [3] [21] . DNA was extracted from whole larvae or adults (30-40 individuals per population) from the F1 generation (also used for bioassays), using the DNeasy tissue Kit (Qiagen) according to the manufacturer's instructions. DNA concentration was determined for each sample with a NanoDrop TM ND-1000 (NanoDrop Technologies), normalized to 50ng/µl and genotyped using the DNA Illumina GoldenGate® array genotyping platform at the Genotyping Center GENTYANE in Clermont-Ferrand, France(http://www.ibisa.net/plateformes/ GENTYANE,101.html). The fluorescent signal intensities of the specific alleles were determined by the BeadExpress Reader (Illumina Inc), and quantified and matched to genotype classes using GenomeStudio v3.1.14 (Illumina Inc) with manual adjustment when necessary. To ensure a good reliability of the genotyping data we used the GenCall50 (GC50) which indicates the reliability of each genotype call, the GenTrain which measures the cluster quality of the SNPs, and the call rate (CR) which represents the fraction of the samples that has been successfully genotyped for a given SNP. We kept SNPs with GC50 higher than 0.35, GenTrain higher than 0.50 and CR higher than 0.50.

Detection of signature of positive selection (Fst outliers)
To detect SNPs displaying higher than expected differentiation under neutrality between the 6 populations (outliers) we used BayeScan v3 (freely available at http:// cmpg.unibe.ch/software/bayescan/). This Bayesian method is based on a logistic regression model in which each logit value of genetic differentiation Fstij at a given locus i in population j is decomposed as a linear combination of αi and βj, corresponding to a locus effect and to a population effect, respectively. To test for the significance of αi for each locus, the posterior probability that αi is different from zero is estimated using a Reversible Jump MCMC algorithm [22]. Priors for αi and βj were Gaussian distributions with means 0 and -2 respectively and standard errors 1 and 1.8 as in [23]. The estimation of model parameters was automatically tuned on the basis of short pilot runs (10 pilot runs, length 2000). Preliminary tests indicated that a burn-in of 10,000 iterations was enough to attain convergence. We used a sample size of 10,000 and a thinning interval of 50, as suggested by [22], resulting in a total length of the chain of 500,000 iterations. Four independent runs were performed to check the consistency between detected outliers. The loci were ranked according to their estimated posterior probability value and all loci with a Bayes factor >3 were retained as outliers.

Association with insecticide resistance
RR 95 and RR 50 were highly correlated (r 2 = 0.92 and 0.90, both P < 0.0001, for temephos and deltamethrin respectively) and we used RR 95 in association analysis in order to maximize the range of RR variation between the populations. Association analyses between RR 95 and SNPs were performed using the least squares fixed effects linear model implemented in the software TASSEL v3 [24]. This model tests for association between genotypes at each segregating marker and phenotype values, while accounting for population structure by using the 3 first principal components of a PCA derived from genetic markers as covariates. P values were corrected for multiple comparisons by permutation tests with 1000 permutations and 0.05 confidence level.

Markers location in the reference genome and candidate genes identification
All SNPs detected either as Fst outliers, or as being significantly associated with a resistance phenotype, were blasted on the reference genome of Ae. aegypti (http:// aaegypti.vectorbase.org/) to locate them inside or at proximity of potential candidate genes. The length of linkage disequilibrium can strongly vary from a genomic region to another and from a population to another, so that it is very difficult to assess an overall disequilibrium length in natural populations. In a closed population (a laboratory Ae. aegypti strain) under a standard pure-drift demographic model, overall linkage disequilibrium was estimated to be 20 Kb for neutral markers [25]. In a genomic region under selection, the amount of linkage disequilibrium can be considerably greater than that predicted for a neutral region. Furthermore, in natural population, gene flow can produce significant levels of disequilibrium. Therefore we have arbitrarily chosen a conservative threshold of 500 Kb to search for candidate genes potentially under selection in the vicinity of outlier markers, and all the genes present up to 500 Kb each direction from the detected SNP were considered as potentially involved in the resistance phenotype and discussed.

Bioassays
Bioassays showed that the twelve tested populations from Martinique were resistant to temephos and deltamethrin in comparison with the laboratory susceptible Bora-Bora strain and the SBE strain from Benin (Table 1). However, the resistance observed against temephos and deltamethrin was heterogeneous in the island. For temephos, RR 50 ranged from 13 fold (SJOS) to 36-fold (GMRN) and RR 95 ranged from 26fold (SJOS) to 153-fold (VCLN). The populations from the two islets showed a lower resistance level to temephos with RR 50 of 5-fold for both islets, and RR 95 of 17 and 21-fold for ILA and ILO, respectively. Mortality after 24 h deltamethrin exposure ranged from 19% for AJPB to 90% for the RSAL population compared to the 100% mortality in the susceptible strains Bora-Bora and SBE. For each insecticide, resistance ratios were not significantly explained by the number of insecticide treatments, or by other environmental variables such as the type of locality (rural, urban), human density or altitude (Table S1).

Microsatellite polymorphism analysis
Genetic variability and Hardy-Weinberg expectations. A total of 31 populations were successfully genotyped with the 7 microsatellites markers that were all polymorphic except 38/38 that was not kept for further analyses. Genetic diversity was relatively homogeneous among the different populations tested at both macro and micro-geographic levels ( Table 2 and Table  S2). Average allelic diversity was low (2.33 alleles/locus, range 1.73-2.38) as compared to the diversity found in African populations [10], but comparable to the diversity found in south American populations [26]. At the macro-geographic scale, when the populations were considered together, important heterozygote deficit were detected (34/72 and H08, and all loci; Table 2), suggesting that Martinique is not a single panmictic unit. However, within localities, departures from Hardy-Weinberg expectations and significant heterozygote deficits were detected only 3 times: in GRMN at locus 34/72, and in AJB and RS when all loci were considered together (Table 2). Overall, these results indicate that each sampled site can be considered as a panmictic population. At the micro-geographic level and considering the microsatellite markers individually, significant deviations were detected only in Fort-de-France samples (Table S2). Genetic diversity indices (Nall, Rs and F IS ) were not correlated (all P > 0.05) with the environmental variables measured (altitude, human density). However, allelic richness was negatively correlated with the level of resistance to deltamethrin (t = -2.66, df = 10, P = 0.02, r = -0.64) underlying a significant decrease of diversity in the most resistant populations. Furthermore, F IS was significantly correlated with the level of resistance to temephos (t = 2.43, df = 10, P = 0.035, r = 0.609) showing an increased consanguinity in populations more resistant to temephos.

Genetic differentiation and population structure
The global genetic differentiation index F ST was low (0.079) but the populations were significantly structured both at the macro-geographic and at the micro-geographic scales (Table  3). Pairwise F ST values were significant for 69% of the population pairs tested. Genetic and geographic distances were correlated at the Martinique scale (Mantel test; P = 0.009) and at the local scale in Vauclin (P = 0.042), but not in Fort-de-France and Sainte-Anne (P > 0.5). Four genetic clusters were identified in Martinique (Figure 1), one being mostly represented in the two islets and their closest eastern Vauclin population on the main island, and another corresponding to populations around the gulf of Fort-de-France (urban area). All the populations in Martinique were highly admixed, suggesting important gene flow across populations. The most homogenous populations were the two islets assigned >80% to one cluster and Fond-St-Denis assigned >90% to another cluster.

SNP diversity
Average allelic diversity per locus per population was 1.76 (range 1.6-1.9) and the 6 populations were at the HW equilibrium except for a slight deficit in heterozygotes in Rivière Salée. The global Fst = 0.031 was low, suggesting large gene flow among these populations; the structure analysis based on SNPs identified the same 4 main clusters as identified based on microsatellites (data not shown). Pairwise Fst were generally smaller for SNPs than for microsatellites (Table S3).

Detection of Fst outliers and association tests
Five out of the 319 SNPs were identified as Fst outliers across the 6 populations tested in Martinique (Table 4), including the kdr-V1016I mutation in the voltage-gated sodium channel. However, the kdr mutation was not detected as significantly associated with deltamethrin resistance. The resistance allele is found at high frequency in all populations on the main island (85% to 92%), confirming previous report [3], which may explain the failure of the GLM regression to detect significant association between deltamethrin resistance and the frequency of the resistant allele. From the five Fst outliers, two were significantly associated with deltamethrin resistance, and one with temephos resistance. Five other markers were associated with deltamethrin resistance and four other were associated with temephos resistance (Table 4). Taken separately, each marker explained from 9 to 17% of the resistance variability observed between populations.

Markers location in the reference genome and candidate genes identification
Most of the SNPs identified were either inside a gene annotated in Vectorbase, or in the flanking 3' extremity of a gene, which is not surprising because SNPs have been designed based on transcriptome sequences anchored on the 3' end of transcripts [27]. The seven SNPs associated with deltamethrin resistance corresponded to six genes encoding respectively a leucin-rich membrane protein, two proteins associated with ribosomes, and three hypothetical proteins of unknown function. All SNPs were synonymous. Among the five SNPs associated with temephos resistance, two were nonsynonymous mutations in an alkaline phosphatase and in a sulfotransferase, two were synonymous mutations in an aldehyde dehydrogenase and an unknown protein, and the last one was in a non-coding region. Genes located on the same supercontig and within ± 500 Kb from a marker associated with resistance are presented in Table S4. Among these, genes encoding enzymes potentially involved in detoxification processes such as transferases, phosphatases and dehydrogenases were identified. Genes encoding transcription factors and miRNA potentially regulating gene expression were also identified.

Population structure at local and global scales
Our study showed a low genetic diversity in Ae. aegypti populations from Martinique compared with levels observed in other parts of the world, especially in native African populations [9,10], but was consistent with a previous study based on isoenzymes in this Island [28]. The loss of genetic diversity observed here can be explained by the founder effect that occurred during the introduction of the mosquito in the Caribbean. Furthermore, this low genetic diversity may have been maintained since its introduction in the 17th century by several factors such as recurrent demographic crashes linked to the specific environmental variations in Martinique (i.e. rainfall and droughts) and more recently by anthropogenic actions such as vector control. Because of the recurrent Dengue epidemics occurring in Martinique, vector control (i.e. source reduction and insecticide applications) is intense and may induce local bottleneck effects and a reduction of the genetic diversity in Ae. aegypti populations. In support to this hypothesis, we found significant correlations between neutral genetic diversity estimates and insecticide resistance ratio: allelic richness was negatively correlated to deltamethrin Table 3. Genetic differentiation estimated with 6 microsatellite markers at two geographic levels.    Despite an overall low genetic differentiation among sites, most pairwise Fst were significant and there was a significant pattern of isolation-by-distance (IBD) at the whole Martinique scale. This suggests that gene flow mostly occurs at a local scale with high admixture occurring among neighbour sites. Furthermore, four genetic clusters were identified, with marked differences in the proportion of assignation of individuals to each cluster from a site to another. Northern populations were mostly assigned to one cluster (red on Figure 1) and Southern populations were more admixed with a predominance of another cluster (green on Figure 1). The genetic drift observed in Fond-Saint-Denis (mostly red cluster) and in the Islets (mostly blue cluster) could be explained by their geographic isolation, with restricted access in the mountainous area and only accessible by boat, respectively. At a more local scale, the IBD pattern was marginally significant in only one out of three sites analysed (Vauclin, P = 0.042), supporting the effect of genetic drift acting on mosquito populations. Genetic drift at local scale could be due to the recurrent insecticide applications leading to local bottlenecks. Furthermore, flight activity of Ae. aegypti is limited to few hundred meters [29] which is likely to increase the neutral genetic differentiation between the populations through genetic drift. No significant correlation was found between diversity estimates and environmental variables such as altitude or human density.

Outlier detection and association with insecticide resistance
The genome scan analysis was performed on 6 populations representative of the diverse mosquito habitats found across Martinique Island, and of the most extreme insecticide resistance ratio values observed. BayeScan detected five outliers out of 319 SNPs, all with positive αi values, indicative of directional selection. These included the knockdown resistance kdr mutation (V1016I), which is a well-characterized mechanism of resistance to pyrethroid insecticides in insects (including Ae. aegypti) and is caused by a point mutation of the para-type sodium channel [30], and three other genomic regions that were associated with resistance to either temephos or deltamethrin. Therefore, only one out of five BayeScan outliers was not related to insecticide pressure. This outlier was not found significantly associated to altitude or to human density, but it might be linked to other environmental or ecological pressures that were not measured in the present study, including climate and/or agricultural practices. Despite a highly diversified landscape in Martinique, and a complex life cycle involving an aquatic larval stage and an aerial adult stage, the mosquito genome appears to be more constrained by insecticide pressure than by any other environmental pressure. Although identified as an outlier by the BayeScan approach, the kdr mutation V1016I was not significantly associated with deltamethrin resistance; this could be due to the fact that this mutation was present at high frequency in all the main island populations (>85%, as also found in [3]), decreasing the power of the association study. Furthermore, other kdr mutations such as the kdr mutation F1534C found at a relative high frequency in Grand Cayman Island and shown to be linked with DDT and permethrin resistance [31] might also be involved in pyrethroid resistance in Martinique, and their occurrence should be searched for in Martinique. Association study revealed that nine additional SNPs not detected as outliers by BayeScan were significantly correlated with resistance to either insecticide. BayeScan is a conservative approach and is sensitive to the level of genetic differentiation across populations: if genetic differentiation is too small, the Reversible Jump MCMC algorithm used to estimate posterior probabilities can fail to reach convergence, and if it is too high, the power to detect outlier loci is low. In the association study, population structure is taken into account as a fixed effect in the general linear model fitting phenotypic to genotypic value at each locus under investigation. As it is a correlative approach, it is not sensitive to population processes such as drift, migration, and demographic fluctuations. On another hand, this correlative approach might detect false positives. We independently detected two markers significantly associated with deltamethrin resistance located at two different positions in the same gene (AAEL014562) suggesting that this association was not discovered by chance. A further way to test for the role of outliers in insecticide resistance would be to compare their relative frequencies in deads and alives after bioassays.

Candidate genes for resistance to deltamethrin and temephos
Among the SNPs significantly associated with resistance, only two were non-synonymous. Both were associated with resistance to temephos, and both involved genes encoding for enzymes. One induced a change from a threonine to a serine in an alkaline phosphatase while the other induced a change from a lysine to an asparagine in a sulfotransferase (Table S4). Alkaline phosphatases (ALP) catalyse the transfer of a phosphate group to water (hydrolysis) or alcohol (transphosphorylation) [32]. ALP-based assays have been used to detect pesticides in the environment, including organophosphates [33] and might thus be involved in their detoxification. Sulfotransferases are known to be involved in phase II detoxification processes by adding a sulphate group to various substrates to facilitate their excretion [34,35]. These two genes (and point mutations) were never described previously as being linked to temephos resistance, but they represent valuable candidates for further functional validation. All the other SNPs significantly associated with insecticide resistance were synonymous, and therefore not associated to protein variations potentially contributing to resistance, but rather in the regulation of gene expression (transcription factors, miRNA) or being in linkage disequilibrium with a functional mutation (hitch-hiking with the adaptive mutation). Because we genotyped individuals from natural populations, the length of linkage disequilibrium is likely to be low, and we restricted our search for candidates to 500 Kb upstream or downstream from each SNP (Table S4). Among candidates, we found several genes potentially involved in detoxification processes such as phosphatases, dehydrogenases and transferases.
In mammals, alcohol and aldehyde dehydrogenases (ADH and ALDH) can contribute to the metabolism of pyrethroids through further processing of the metabolites generated by esterases, phenoxybenzyl alcool and phenoxybenzyl aldehyde, to phenoxibenzoic acid [36,37]. However, no classical detoxification enzymes found overtranscribed in insecticide-resistant populations from Martinique were identified [3,4]. Indeed, recent data suggest that cytochrome P450 monoxygenases play a crucial role in pyrethroid detoxification in mosquitoes [38] while carboxylesterases are known to confer resistance to organophosphates [39]. It is possible that mutations causing the over-transcription of these detoxification genes in Martinique are located too far away from our markers, or that regulators of detoxification enzymes are located far away from the gene they regulate (trans-regulators); indeed, we identified potential regulators of transcription situated close to SNPs associated with resistance (Table S4). Although this genome scan provides valuable data to identify insecticide resistance genes in mosquitoes, the main limits are the partial annotation of the Ae. aegypti genome, with many genes annotated as 'hypothetical proteins' of unknown function, and the poor assembly of the genome with 4,758 supercontigs. As a result, we restricted our search to genes located 500 Kb in either direction in the same supercontig as the SNP associated to resistance, a relatively short distance, that was even reduced for markers located closer than 500 Kb to supercontig boundaries (see Table S4).
Finally, there was no SNP or genomic region associated with both resistance to deltamethrin and to temephos, which suggests that different resistance mechanisms are evolving in natural populations against these two commonly used insecticides, temephos being used against larvae and deltamethrin against adults. These two insecticides have different targets, the acetylcholinesterase and the sodium channel respectively, and the finding of different genes involved in the two types of resistance is not surprising. Functional validation of the precise role in each type of insecticide resistance for each candidate gene detected in this study is the next step to be undertaken.