The impact of agrochemical pollutant mixtures on the selection of insecticide resistance in the malaria vector Anopheles gambiae: insights from experimental evolution and transcriptomics

Background There are several indications that pesticides used in agriculture contribute to the emergence and spread of resistance of mosquitoes to vector control insecticides. However, the impact of such an indirect selection pressure has rarely been quantified and the molecular mechanisms involved are still poorly characterized. In this context, experimental selection with different agrochemical mixtures was conducted in Anopheles gambiae. The multi-generational impact of agrochemicals on insecticide resistance was evaluated by phenotypic and molecular approaches. Methods Mosquito larvae were selected for 30 generations with three different agrochemical mixtures containing (i) insecticides, (ii) non-insecticides compounds, and (iii) both insecticide and non-insecticide compounds. Every five generations, the resistance of adults to deltamethrin and bendiocarb was monitored using bioassays. The frequencies of the kdr (L995F) and ace1 (G119S) target-site mutations were monitored every 10 generations. RNAseq was performed on all lines at generation 30 in order to identify gene transcription level variations and polymorphisms associated with each selection regime. Results Larval selection with agrochemical mixtures did not affect bendiocarb resistance and did not select for ace1 mutation. Contrastingly, an increased deltamethrin resistance was observed in the three selected lines. Such increased resistance was not majorly associated with the presence of kdr L995F mutation in selected lines. RNA-seq identified 63 candidate resistance genes over-transcribed in at least one selected line. These include genes coding for detoxification enzymes or cuticular proteins previously associated with insecticide resistance, and other genes potentially associated with chemical stress response. Combining an allele frequency filtering with a Bayesian FST-based genome scan allowed to identify genes under selection across multiple genomic loci, supporting a multigenic adaptive response to agrochemical mixtures. Conclusion This study supports the role of agrochemical contaminants as a significant larval selection pressure favouring insecticide resistance in malaria vectors. Such selection pressures likely impact kdr mutations and detoxification enzymes, but also more generalist mechanisms such as cuticle resistance, which could potentially lead to cross-tolerance to unrelated insecticide compounds. Such indirect effect of global landscape pollution on mosquito resistance to public health insecticides deserves further attention since it can affect the nature and dynamics of resistance alleles circulating in malaria vectors and impact the efficacy of control vector strategies. Supplementary Information The online version contains supplementary material available at 10.1186/s12936-023-04791-0.


Background
Malaria is one of the deadliest mosquito-borne diseases in Africa mainly transmitted by Anopheles spp.[1,2].The recent World Health Organization (WHO) report shows 619,000 deaths due to malaria in Africa in 2021 [3].Malaria control includes prevention through insecticide-based vector control in order to limit the transmission of pathogens [4].Five classes of insecticides are used in malaria vector control but pyrethroids are the most widely used [4][5][6][7].Their use relies essentially on two indoor methods: indoor residual spraying (IRS) and long-lasting insecticidal nets (ITNs).IRS and ITNs showed success between 2000 and 2015 by reducing the number of malaria cases in several countries [8].Over the years however, the use of chemical insecticides has become less effective against malaria vectors, which now show strong resistance to most public health insecticides [9].
Insecticide resistance is widespread in insects [10].Previous studies have shown that resistance in mosquitoes is mainly due to physiological adaptations, which include target-site mutations, increased insecticide metabolism or sequestration and altered insecticide penetration [11,12].Target site resistance is caused by non-synonymous mutations affecting the neuronal proteins targeted by insecticides [13,14].Such mutations are highly conserved in insects and well known in mosquitoes and their genotyping provides useful information to track resistance in the field [13][14][15][16][17]. Multiple target-site mutations conferring resistance to public health insecticides have been identified in African malaria vectors [18,19].The target protein of both pyrethroids and DDT is the voltage-gated sodium channel (VGSC) which is affected by knock-down resistance (kdr) mutations.A leucine-phenylalanine substitution at position 995 (L995F, orthologous to Musca domestica Vgsc codon 1014) [20], was first identified in West Africa and, therefore, named the 'kdr West' mutation [13].A second mutation at the same codon (L995S), was also identified in East Africa.Nowadays, there is increasing evidence of the spread of the L995S mutation in West Africa and vice-versa [10,[21][22][23][24].The targets of both organophosphate and carbamate insecticides are the acetylcholinesterases which include the ace1 gene affected by the G119S resistance mutation in Anopheles spp.[25][26][27][28].
Beside target-site mutations, increased insecticide metabolism (i.e.metabolic resistance) has been reported in several African countries [29].Such resistance phenotype is often caused by an increased activity of detoxification enzymes including cytochrome P450 monooxygenases, carboxylesterases, glutathione S-transferases, and UDP glycosyltransferases, though other protein families can also be involved [11].Such increased activity is often caused by gene over-expression though structural modifications contributing to enhanced insecticide metabolism populations have also been identified in malaria vectors [29].Given the high diversity and functional redundancy of insect detoxification enzymes, the identification of those conferring resistance to insecticides proved to be challenging [30].Although target-site mutations and metabolic resistance clearly play a major role in conferring resistance to public health insecticides, other physiological changes such as altered insecticide penetration caused by cuticle thickening or structural modification have also been shown to contribute to resistance [31].Indeed cuticular proteins like CPAP3E, CPLCG4 and CPLCG5 have been frequently associated with mosquito resistance through cuticle thickening [32][33][34].Genes coding for P450 cytochromes of the Cyp6M and Cyp6Z families, together with GSTE2, are involved in mosquito resistance to insecticides [35][36][37].
Though the mass distribution of LLINs and IRS acted as a major selection pressure leading to pyrethroid resistance in Africa [38,39], the intensive use of the same insecticide families for crop protection also represents a significant selection pressure undergone by Anopheles populations located in agricultural areas [40].Indeed, mosquito larvae found in these ecosystems are exposed to a wide range of pesticides used against crop pests while adults may also be impacted by agricultural spraying operations [41,42].In addition, non-insecticide molecules such as herbicides and fungicides can also have adverse effects on mosquitoes and may contribute to resistance selection through chemical stress response mechanisms [43].Overall, an increasing number of studies support the role of agricultural xenobiotics as a key selection pressure contributing to insecticide resistance in mosquitoes, though the underlying mechanisms have not been fully characterised [44,45].
In this context, the aim of this work was to explore the potential of agrochemical mixtures present in Anopheles gambiae breeding sites to select for inherited resistance to vector control insecticides at the adult stage.A field derived An. gambiae line was experimentally selected at the larval stage for 30 generations with three different mixtures containing agrochemicals commonly used in agriculture in Africa: (i) a mixture of insecticide-based formulations, (ii) a mixture of fungicide and herbicide formulations, and (iii) a mixture containing both insecticide and non-insecticide formulations.Comparative bioassays with the pyrethroid deltamethrin and the carbamate bendiocarb were used to monitor the resistance of adults to vector control insecticides across generations.The impact of each selection regime on resistance mechanisms was investigated by genotyping target-site mutations and whole transcriptome analysis.Results are discussed in regards to the impact of agriculture on the management of insecticide resistance in malaria vectors.

Mosquitoes
A field-derived colony of An. gambiae (form S) from Tiassalé in southern Côte d'Ivoire was used as the parental strain in the present study.This colony, hereafter Tiassale-S, has been maintained since 2015 at the Centre Suisse de Recherches Scientifiques (CSRS) in Côte d'Ivoire without selection and now exhibits a phenotype of low resistance to public health insecticides compared to contempory field mosquitoes from the locality of Tiassalé.The mosquitoes were reared under standard tropical rearing conditions (27 ± 3 °C and 75 ± 10% humidity under 12:12 photoperiod).Larvae were fed on cat food and adults on 5% honey solution.Adult females were blood fed membrane feeding (Hemotek).

Controlled selection
The parental line was divided into four distinct lines, each subjected to a different selection regime.The Control line (Cntrl line) was maintained without selection pressure and served as control in all experiments.The Insecticide line (Ins line) was produced by exposing L2 larvae to a mixture of five insecticides from the pyrethroid, neonicotinoid, carbamate and organophosphate chemical families (Table 1).The Non-insecticide line (Non-ins line) was produced by exposing L2 larvae to a mixture of five noninsecticide molecules, two fungicides from the organochlorine and carbamate families and three herbicides from the amide, pyridine and aminophosphonate families.The mixture line (mix line) was selected using the mixture of Insecticide and Non-insecticide formulations.For each line, selection for resistance in mosquitoes was performed by exposing stage II larvae (~ 300) in tanks containing 300 mL of water and agrochemical mixtures for 24 h to the different agrochemical mixtures.Surviving larvae were rinsed and transferred to clean tap water and reared.As soon as they emerged, the adults of both sexes, which were placed in wire cages, were allowed to mate freely and membrane feeding (Hemotek) was used to generate the eggs of the next generation.Selection was carried out over 30 successive generations using a dilution killing 20% of the larvae (LD 20 ) in the different lines (see Table 2).The selection pressure was maintained around the LD 20 throughout the selection process by adjusting the doses every five generations.The new LD 20s were determined using PoloPlus software.

Insecticide resistance monitoring
The resistance level of each line to bendiocarb and deltamethrin, two insecticides used in vector control, was monitored at the adult stage every five generations.Bioassays were carried out according using test tubes with filter papers impregnated with either 0.1% bendiocarb or 0.05% deltamethrin [46].At least four batches of 20 to 25 non-blood fed 2-5 days old females were used.Mortality was recorded after one hour of insecticide exposure and a 24 h recovery time during which the mosquitoes were provided a 5% honey solution.The Abbot formula correction was applied when the control mortality rate was between 5 and 20% and assays were discarded if mortality in control exceeded 20% [47].The mortality of each selected line to each insecticide was compared to that of the unselected line at the same generation using a Fisher test (N ≥ 4).

Target-site mutations
The frequencies of the kdr West L995F and ace1 G119S target-site mutations were monitored in each line by individual genotyping at generations G0, G10, G20 and G30.Genomic DNA was extracted using 2% CTAB as previously described [48].Kdr L995F and ace1 G119S mutations were genotyped using the allele-specific TaqMan qPCR methods as described.Each reaction mixture contained 5 μL of 2X sensimix (Bio Rad), 3.875 μL of nuclease-free water, 0.125 μL of TaqMan probes and 1 μL of gDNA.Quantitative PCR reactions were performed on a CFX 96 Real Time system (Bio-Rad technologies, California, USA) with the following amplification conditions: 95 °C for 10 min, followed by 40 cycles of 95 °C for 10 s and 60 °C for 45 s.For each mutation, individuals were scored as homozygous susceptible/resistant or heterozygous based on the intensity of the HEX/FAM channels at the end of the PCR reaction as compared to positive and negative samples of known genotypes.For each generation, the kdr and ace1 genotype frequencies were compared between selected and unselected lines using Genepop sofware 4.0.10 based on a chi 2 test.

RNA-seq library preparation and sequencing
The transcriptome of each selected line was compared to the control line using RNA-seq at generation G30.For each line, four pools of 30 three-day-old non-blood fed females (not exposed to insecticide) were collected and stored in RNA-later at − 20 °C.Total RNA was extracted from each pool using Trizol (Life Technologies) according to manufacturer's instructions, and treated with DNase to remove genomic DNA contaminants.RNAseq libraries were prepared from 150 ng total RNA using NEBNext ® Ultra ™ II Directional RNA library Prep Kit for Illumina (New England Biolabs) following manufacturer's instructions.Libraries were quantified using the Qubit DNA BR assay (Thermofisher Scientific) and quality checked using Bioanalyzer DNA 1000 assay (Agilent).
Libraries were sequenced in multiplex as single 75 bp reads on a NextSeq 500 sequencer (Illumina) by Helixio (Clermont-Ferrand, France).After demultiplexing and quality check using FastQC, reads were loaded into Strand NGS V 3.2 (Strand Life Sciences) and mapped against the AgamP4 assembly and AgamP4.12geneset using the following parameters: min identity 90%, max gaps 5%, min aligned length 35 bp, ignore reads with more than 5 matches, trim 3' ends of reads with average quality < 20, Kmer size 11, match score 1, mismatch score 4, gap opening penalty 6, gap extension penalty 1. Mapped reads were then filtered based on their sequence quality and mapping quality as follows: Mean read quality ≥ 20, number of N ≤ 5, alignment score ≥ 90, mapping quality ≥ 120, number of match = 1.The remaining reads (~ 90% of sequenced reads) were used for subsequent analyses.

Differential gene transcription
Differential transcription analysis was performed on Strand NGS V3.2.This analysis was performed on all protein coding genes with normalisation and quantification being performed according to the DESeq algorithm [49].Only the 10357 genes showing a coverage ≥ 4 reads/ kb in all replicates across all conditions were kept for further analysis.Transcription levels between each selected line and the control line were then compared across the four biological replicates using an ANOVA followed by a Tukey HSD test.P values were adjusted for multiple testing corrections using the false discovery rate method [50].Genes showing a transcription ratio ≥ 1.5 fold in either direction and a P value ≤ 0.005 in any selected line as compared to the parental line were considered differentially transcribed following selection.For each selected line, genes significantly over-and under-transcribed as compared to the control line were subjected to a Gene Ontology (GO) term enrichment analysis using the functional annotation tool DAVID [51].Reference gene list was constituted by the 10357 genes detected by RNA-seq and test lists were constituted from over-and under-transcribed genes in each line.GO term frequencies in selected vs unselected line were compared in a Fisher's Exact test, and terms with a P value < 0.05 upon FDR multiple testing correction were kept [50].A panel of 193 genes were selected from Agam P4.12 geneset as candidates possibly contributing to xenobiotic resistance.These genes included known insecticide targets, detoxification enzymes (cytochrome P450s, carboxylesterases and transferases), ABC-transporters, cuticle proteins, enzymes associated with redox stress, nervous receptors and putative insecticide binding proteins (see Additional file 4: Table S1).Heat maps reflecting transcription profiles of differentially expressed resistance candidate genes across all lines were generated using TM4 Multi-experiment Viewer (MeV) software [52].

Sequence polymorphism
Small Nucleotide Polymorphisms (SNPs) were called from transcriptomic sequence data using Strand NGS V 3.2 against all protein-coding genes of the AgamP4.12geneset using standard parameters (ignore homopolymer stretches greater than 4 bp and adjacent positions, coverage ≥ 30 and ≤ 5000, reads supporting the variant allele ≥ 2, base quality ≥ 20, variant confidence score ≥ 200 and strand bias ≤ 25).Among the variations called, only substitutions and indels were retained for further analyses.A principal component analysis [53] was then used to visualise the genetic divergence of each line to the AgamP4 reference genome.PCA was performed on the frequency of all bi-allelic variations identified in each replicate of all lines using the Ade4 R package [54].Genic effects were then computed based on the longest transcript for each gene according to the AgamP4.12geneset and sorted as affecting (non-synonymous) or not (synonymous) the protein sequence.Selection signatures were investigated using the bi-allelic SNPs that were polymorphic (i.e.showing a > 5% allele frequency variation between the control parental line and at least one selected line) in two ways.First, a selection was made of differential SNPs, i.e. those showing a clear difference in frequency between the selected and non-selected lines.For this purpose, the mean variant frequencies were compared in a Student's T test, and SNPs with a P-value < 0.0005 after FDR correction were retained.This stringent P-value threshold still retained 2.5 to 4.8% of the SNPs, depending on the selected line.A SNP score was then computed for each differential SNP based on its absolute frequency variation between the selected line and the control line.The score of non-differential SNPs was set to 0 while differential SNP scores were calculated as follows: Score = Abs[(%freq Selected ) − (%freq control )]/50, where '%freq' is the frequency in % of the variant allele.In this way, an allele showing a 50% frequency variation in a selected line scores 1, and an allele absent in the control line and fixed in a selected line scores 2. A gene Diff score was finally computed by summing SNP scores and dividing by the number of polymorphic SNPs in each gene.A second approach consisted in assessing F ST departure from neutrality using the Bayesian method implemented in BayeScan version 2.1 [55].A separated analysis was performed consisting in contrasting the selected line versus the control line across their four replicates.Default settings were used except that prior odd was set to 1000 in order to increase stringency.SNPs showing a Bayescan Q-value of zero were considered as 'Outliers' .Outliers represented 2.9 to 4.1% of all SNPs, depending on the selected line, and were counted per gene.The percentage of outlier SNPs in each gene was then plotted along chromosomes.

Insecticide resistance dynamics during the selection process
The resistance of adults to deltamethrin and bendiocarb was monitored in each line during the selection process (Fig. 1).High mortality rates to deltamethrin (94.1%) and to bendiocarb (92.7%) were obtained with the parental line at G0, supporting the low frequency of resistance alleles at the beginning of the selection process.Larval selection by the agrochemical mixtures had no impact on bendiocarb resistance, despite the presence of carbamates and organophosphates (both targeting the acetylcholinesterase) in the agricultural insecticide mixture.Conversely, an increased deltamethrin resistance was observed in response to larval selection with agrochemical mixtures.This increased resistance was significant in all selected lines from G5 onwards and continued to increase up to G30.The highest resistance level was reached in the ins line selected with the insecticide mixture (29.8% mortality at G30) while the two other nonins and mix lines were less resistant at G30 (non-ins line 54.74% and mix line 47.36% mortality at G30).

Target-site mutations
The evolution of kdr mutations affecting the voltagegated sodium channel and conferring resistance to pyrethroids and DDT was monitored through the selection process.The kdr East (L995S) mutation was not detected in the parental line and thus not further quantified during selection.The kdr mutation L995F was present in the parental line at a frequency of 60% at G0 (Fig. 2).Its frequency significantly decreased to reach 47% in the non-selected line at G30. Conversely, the frequency of the kdr L995F mutation remained stable in all selected lines (G30 frequencies of 63%, 65% and 62% in the Ins, Non-ins and Mix lines respectively).Genotypes carrying the kdr L995F mutation were significantly more represented in all selected lines as compared to the control line at G30.The ace1 mutation G119S affecting the acetycholinesterase and conferring resistance to carbamates and organophosphates was present in the parental line at a low initial frequency (17%).Its frequency gradually decreased through generations in all lines selected or not with agrochemical mixtures (Fig. 3).At G30, the final ace1 mutation frequencies were 0%, 0%, 5% and 0% in the ins, non-ins, mix and control lines respectively.However, a significant difference of genotype frequency was found between the control line and the non-ins line at G10 (chi 2 test P < 0.05).

Differential gene expression
Among the 10357 genes detected by RNA-seq, 775 were considered as significantly differentially transcribed in at least one selected line as compared to the control line (≥ 1.5 fold-change in either direction and corrected P value ≤ 0.005, Additional file 4: Table S1).A total of 472 genes were over-transcribed in at least one selected line with 294, 299 and 216 identified in the Ins, Non-ins and Mix lines respectively (Additional file 1: Figure S1).Among them, 111 genes were shared by two selected lines and 113 genes were shared by the three selected lines.Only 304 genes were under-transcribed in at least one selected line with 182, 156 and 177 identified in the Ins, Non-ins and Mix lines, respectively.Among them, 83 genes were shared by two selected lines and 64 genes were shared by the three selected lines.
Gene ontology enrichment analyses identified biological processes enriched from genes significantly over-and under-transcribed in each selected line (Additional file 2: Figure S2).Only a few GO terms were found significantly Fig. 1 Evolution of adult resistance to bendiocarb and deltamethrin in each selected line compared to the control line.For each line, insecticide resistance levels are shown as % mortality to 0.1% bendiocarb and to 0.05% deltamethrin.* Indicate significantly distinct mortalities using Fisher's exact test (P < 0.05), and error bars show SD in means 80 < n > 100 computed from each replica tube enriched from under-transcribed genes with a large overlap across the three selected lines.These included two terms associated with endopeptidase activity enriched in all lines and the GO term 'hydrolase activity' identified in both the Ins and the Mix lines.Multiple GO terms were found significantly enriched from over-transcribed genes, with again a good overlap across the three selected lines.These essentially included terms associated with P450 activity and detoxification ('oxidoreductase activity'; 'monooxygenase activity'; 'iron ion binding'; 'haem binding'; 'flavin adenine dinucleotide binding') together with terms associated with the insect cuticle ('structural constituent of cuticle';'chitinase activity'; 'chitin binding').
A total of 63 candidate genes potentially involved in insecticide resistance were over-transcribed in at least one selected line while only 23 candidate genes were found under-transcribed (Additional file 1: Figure S1).Among over-transcribed candidate genes, half were over-transcribed in at least two selected lines, including 14 over-transcribed in the three selected lines.Overtranscribed candidate genes include 17 P450s, 3 carboxylesterases, 10 transferases, 3 ABC transporters, 17 cuticle proteins and 13 other candidates (Fig. 4).Most P450s showed an over-transcription in the Ins line with three of them (CYP6P3, CYP6M2 and CYP6Z2) being known as able to metabolise insecticides [37,56,57].Four P450s (CYP9M1, CYP325D1, CYP12F4 and CYP4H25) were over-transcribed in the three selected lines, with CYP12F4 also showing a significant selection signature (see below).Other over-transcribed detoxification genes include four GSTs, five UDPGTs, one sulfotransferase, three carboxylesterases, three ABC transporters and other enzymes including an aldehyde oxidase and an epoxide hydrolase.GSTE2, known as able to metabolise DDT, was over-transcribed in the Ins line.Most of the 17 over-transcribed cuticle proteins were identified in multiple selected lines with CPLCX2, CPLCG5 and CPLPCP10 over-transcribed in all selected lines, and CPLCG5 previously shown to play a key role in cuticle resistance [34].Three cuticle proteins (CPR130, CPLCG4 and CPLCX3) were also associated with genomic selection signatures (see below).Among genes likely involved in response to oxidative stress, three haem peroxidase (HPX3, HPX5 and HPX12) and one thioredoxin peroxidase were over-transcribed in one or multiple selected lines.Finally, two cholesterollike transporters (Niemman-Pick type C2 proteins,

Polymorphism variations
More than 60 K SNPs were detected across all lines.When projected by Principal Component Analysis (PCA), the three first axes accounted for 96.3% of the total variance (Additional file 3: Figure S3).This included > 80% for the first axis which did not separate the selected and unselected lines and rather reflected the polymorphism between the parental line and the reference genome.The second axis (12.5%) opposed the control line to all selected lines in a balanced manner, supporting a common adaptive response to xenobiotics.Finally, the third axis (2.6%) opposed the Ins line and the Non-Ins line, implying specific components of adaptive response of lesser magnitude; the Mix line was located in-between, supporting its intermediate adaptive response.
The 50 K SNPs that were polymorphic between the control line and at least one selected line were used (Additional file 5: Table S2).As expected from transcriptomic data, most of these SNPs fell (98%) within gene boundaries, covering 4406 genes (i.e.42% of all RNA detected genes).Differential SNPs (Diff SNPs) were defined as those whose frequency varied significantly between any selected line and the control line (see "Methods").These include 3744, 4280 and 3862 Diff SNPs for the Ins, Non-ins and Mix lines, respectively.Summing up the Diff SNP scores supported a highest genetic divergence from the control line for the Ins line (total weight ~ 100), followed by the Mix line (total weight ~ 90) and the Non-ins line (total weight ~ 70).The Fst-based approach identified 2179, 1554 and 1994 Outlier SNPs, in the Ins, Non-ins and Mix lines, respectively, supporting the same divergence ranking between the selected lines as compared to the control line.Diff SNP and Outlier SNP densities often coincided between the two approaches, revealing multiple regions potentially under selection (Fig. 5).These regions often coincided across selected lines though no decreased genetic diversity was observed in the control line (53 K SNPs detected in the control line versus 45 K to 47 K in the selected lines), rather supporting a common multi-genic adaptive response to chemical stress than drift in the control line only during insectarium rearing.Candidates genes affected by Diff/Outlier SNPs in regions showing shared selection signatures included three transporters (ABCB7, ABCB4 and ABCF3), three UDPGTs (AGAP006775  2R).Among the 14 P450s of the CYP6M cluster on Chr 3R (which includes P450s known to metabolise insecticides and other xenobiotics) [56,58], three of them were affected in the Ins or the Mix lines (CYP6Y1, CYP6M4 and CYP6Z1) but not in the Non-ins line.Another interesting signal was observed on Chr 3R in a gene cluster containing 28 cuticle proteins (which includes CPLCG5 known to contribute to pyrethroid resistance) and from which the neighbouring genes CPLCG3 and CPLCG4 were identified in multiple selected lines.Finally, no selection signature was observed in the vicinity of the VGSC gene AGAP004707 (containing kdr mutations), but its low expression level prevented the detection of polymorphic SNPs in this region.

Discussion
The role of agriculture in the development of insecticide resistance in mosquitoes is increasingly being recognized.Indeed, laboratory work has shown the potential of agrochemicals to induce or select for an overexpression of resistance genes in both larvae and adults [44,59].In the field, susceptibility testing of adult mosquitoes in areas of high agricultural activity often revealed an increased resistance associated with the expression of detoxification enzymes [19,42].This is because some mosquito larvae developing in these ecosystems are subject to selection pressure as well as adult mosquitoes present during agricultural spraying campaigns [60].The long term impact of agriculture on the selection of resistance in mosquitoes is also supported by the fact that insecticides used in agriculture are often similar (same families and modes of action or same molecules) as those used in public health [40].Although non-insecticidal agrochemicals were also shown to affect mosquito tolerance to insecticides [61,62], the effect of complex agrochemical mixtures has less been studied.In this context, the present study aimed at combining controlled selection and molecular approaches to study the impact of larval selection by insecticide and non-insecticide agrochemical mixtures on the selection of insecticide resistance mechanisms in An. gambiae.

Larval selection with agrochemical mixtures select for increased resistance in adults
The present study confirmed that larval selection with agrochemicals can lead to an increased resistance to pyrethroids in adults, and also that agrochemical formulated products not sold as insecticides can indeed kill insects.Resistance to deltamethrin increased over the generations in all selected lines (ins line, non-ins line, mix line).Such increased resistance was not associated with a significant increase of the L995F kdr mutation affecting the voltage-gated sodium channel targeted by pyrethroids [63][64][65].However, one should note that this kdr mutation significantly decreased through generations in the non-selected line while its frequency remained stable in selected lines.Though drift effect may have occurred, this supports fitness costs associated with this mutation in absence of insecticide [66].Hence, the stable frequency observed in all selected lines may indicate that insecticide and non-insecticide mixtures still exerted a moderate selection pressure on the VGSC.In An. gambiae, the L995F 'kdr' mutation has been widely observed in association with pyrethroid and DDT resistance throughout Africa [39,63,67].In addition, a high frequency of this mutation is also frequently observed in intensive agricultural areas where crop protection strategies mainly rely on the use of chemical insecticides [18,19,41,68].The maintaining of the kdr L995F mutation in the noninsecticide line may be explained by the presence of an organochlorine (chlorothalonil, sold as a fungicide) in the non-insecticide agrochemical mixture which alike pyrethroids and DDT may exert a selection pressure on insect VGSC [69].Altogether, the early rise of deltamethrin resistance in selected lines can hardly be explained by the presence of the L995F kdr mutation suggesting that other deltamethrin resistance alleles were selected by pesticide mixtures.
In contrast to deltamethrin and despite the presence of carbamates in both the insecticide and non-insecticide mixtures (carbofuran and carbendazime), no resistance of adults to bendiocarb was observed.This absence of bendiocarb resistance was associated with a slow decrease of the G119S ace1 mutation commonly associated to carbamate and organophosphate resistance in An. gambiae [17,24,25,36,70].This trend might be explained by the low initial frequency of the ace1 mutation in the parental line and its significant fitness cost [17,71].This trend may also indicate a lower selection pressure exerted by carbamates and organophosphates present in the agrochemical mixtures as compared to other insecticides, such as pyrethroids and organochlorines.In accordance with this, the resistance of adult mosquitoes to bendiocarb and malathion is often less marked in agriculture intensive areas as compared to DDT and pyrethroids [18,19].In addition, it is likely that negative interactions occurred between the different chemicals present in the insecticide mixture, which might have decreased the selection pressure exerted by carbamates.In An. gambiae, negative interference between insecticides from different families has been shown with opposite effects on a key detoxification enzyme [56].

Larval selection with agrochemical mixtures selects for a broad chemical stress response in adults
GO terms enrichment analysis showed a marked enrichment in molecular functions associated with xenobiotic detoxification or cuticle exoskeleton in all selected lines.This trend was also evident from the large overlap of over-transcribed candidate genes between the three selected lines.Such a broad adaptive response to different chemical mixtures was also supported by polymorphism data showing common selection signatures between the three selected lines and low genetic distances between them as inferred by PCA.It is very likely that the common variations observed in the selected lines are adaptative responses rather than caused by genetic drift events affecting solely the unselected line, because polymorphism rate was still higher in the control line than in selected lines.Altogether, both gene expression and polymorphism data support the selection of a broad and generalist response to chemical stress affecting multiple loci and acting on various traits such as xenobiotic penetration and metabolism.
Among phase I detoxification enzymes, several P450s were over-transcribed in one or multiple selected lines.These include key resistance genes like CYP6P3 or CYP6M2 whose role in pyrethroid resistance has been functionally or genetically validated [58,72] together with other P450s (e.g.CYP12F4, CYP6Z2 and CYP9M1) previously associated with resistance using laboratory or field approaches [73,74].The response of P450s to larval selection with agrochemical mixtures is further supported by the selection signature observed at the CYP6M resistance locus on Chromosome 3R which contains 14 CYP6 genes from the CYP6Y, CYP6M and CYP6Z subfamilies.
Among phase II enzymes (transferases), multiple GSTs and UDPGTs were over-transcribed or affected by selection signatures following selection while agrochemical mixtures.An over-transcription of GSTe2 was detected upon selection with the insecticide mixture together with a positive selection signature at the GSTE locus.The role of GSTe2 in DDT resistance has been demonstrated in An. gambiae and other mosquito species [35,75,76].This gene and other epsilon GSTs have also been implicated in resistance to various insecticides including pyrethroids and organophosphates in mosquitoes [77] supporting their adaptive role toward various insecticides.Multiple UDPGTs were over-transcribed in one or multiple selected lines while others were associated with selection signatures.These phase II conjugating enzymes are thought to play a key role in xenobiotic detoxification pathways in most organisms [78][79][80].In insects including mosquitoes, the association of UDPGTs and P450s in pyrethroid metabolism pathways has frequently been observed [75,79,81] but their role in the detoxification of other agrochemicals is likely.Other proteins likely contributing to xenobiotic metabolism were over-transcribed and/or affected by selection signatures in selected lines.This includes various phase I enzymes such as aldehyde oxidase or epoxide hydrolase [12,82] but also multiple ABC transporters known to contribute to the excretion of xenobiotics and their conjugated metabolites [83,84].In mosquitoes, ABC transporters have been frequently associated with pyrethroid resistance though their complexity and late positioning in detoxification pathways makes their functional validation challenging [85,86].Among non-enzymatic binding proteins, sensory appendage proteins (SAP) appear as likely involved in the broad defence against xenobiotics, since a selection signature appears at the SAP locus in all selected lines.SAP proteins were shown to bind various xenobiotics, among which SAP2 was shown to confer pyrethroid resistance in An. gambiae [87].Also, two Niemann Pick type C2 (NPC2) genes were strongly upregulated in the three selected lines.NPC proteins, initially identified as cholesterol-like transporters, have been suggested to bind xenobiotics and might therefore contribute to their sequestration [88,89].Interestingly, the neurodegenerative Niemann-Pick type C (NPC) disease caused by mutations in NPC genes was also associated with a defective P450-mediated drug metabolism in mouse supporting a cross talk with detoxification pathways [90].Xenobiotic response has also been associated with a higher tolerance to oxidative stress in various insects including mosquitoes [91][92][93].Such response was also apparent in the transcriptomic dataset with multiple red/ox enzymes (haem peroxidases, thioredoxin peroxidase, superoxide dismutase) being differentially transcribed in selected lines.
Finally, several structural cuticle protein genes were over-transcribed in the selected lines.Expression patterns were relatively conserved between lines, supporting the hypothesis of a generalist adaptation to chemical stress.Among over-transcribed genes were two members of the CPLCG gene cluster located on chromosome 3R, including CPLCG5 known to play a key role in pyrethroid resistance [34].This gene cluster also shows a clear selection signature upon insecticide mixture selection.A cuticular component of xenobiotic resistance is further supported by the over-transcription of chitin synthase, an enzyme playing a key role in cuticle formation [94].Deciphering whether the over-transcription of these multiple cuticle proteins is associated with physiological cuticle alteration (cuticle thickening and/or altered insecticide penetration) in selected lines deserves further work.

Agrochemicals as a key selection pressure contributing to insecticide resistance in malaria vectors
Overall, the present work confirms that agrochemical mixtures contaminating mosquito breeding sites represent a significant selection pressure enhancing the ability of adult mosquitoes to resist vector control insecticides.The deltamethrin resistance phenotype observed upon selection with agrochemical mixtures was associated with a broad adaptative response to chemical stress involving detoxification-and cuticle-related pathways.Such multigenic adaptation to chemical stress was largely conserved among the different selection regimes suggesting that a large proportion of selected genes do not specifically respond to a particular agrochemical but were rather selected by multiple compounds.
One major difference between selection pressures represented by vector control versus agriculture stands on their specificity: chemical insecticides used for vector control can be considered as a specific selection pressure (limited number of active ingredients used at a time, most of the time as a single product) while chemical used in agriculture likely represent a broader selection pressure (higher diversity of active ingredients used sequentially or as mixtures).In addition, different agrochemicals (or their metabolites) may accumulate in mosquito breeding sites leading to the exposure of mosquito larvae to complex xenobiotic mixtures [40,42].In such situation, the higher complexity of agriculturebased selection pressures likely favours the selection of generalist resistance mechanisms (i.e.broad spectrum detoxification enzymes, sequestration proteins, cuticle modifications) as opposed to more specific resistance mechanisms (e.g.target-site mutations and a few detoxification enzymes) that are often selected by vector control interventions [38,39,44,59].Though no crossresistance between deltamethrin and bendiocarb was observed in the selected lines, the diversity and generalist nature of the resistance alleles selected by agrochemical mixtures agrees with the multi-resistance phenotypes frequently observed in intense agriculture areas [45,95].

Conclusion
Altogether, the present study confirms that mosquitoes can cope with the variety of anthropogenic xenobiotics encountered in their larval environment through multigenic adaptive trajectories which in turn may impact various adult traits including resistance to vector control insecticides.Given the limited number of active ingredients available for public health, this may have a significant impact on the management of resistance in malaria vectors and calls for an integrated management of resistance between agriculture and vector control.Whether other key vector biological functions (e.g.reproduction, development, aging, behaviour) are impacted by agrochemicals may deserve further attention as this may affects the global ecology of vectors and malaria transmission though Africa.
Additional file 5: Table S5.Gene-level polymorphism data obtained from RNA-seq.
• fast, convenient online submission • thorough peer review by experienced researchers in your field • rapid publication on acceptance • support for research data, including large and complex data types • gold Open Access which fosters wider collaboration and increased citations maximum visibility for your research: over 100M website views per year

•
At BMC, research is always in progress.

Learn more biomedcentral.com/submissions
Ready to submit your research Ready to submit your research ?Choose BMC and benefit from: ? Choose BMC and benefit from:

Fig. 4
Fig. 4 Expression profiles of candidate resistance genes in each selected line.Gene transcription levels were quantified by RNA-seq after 30 generations of selection.Only genes showing a significant differential transcription level between at least one selected line and the control line are shown (*: FC ≥ 1.5-fold in either direction and corrected P value ≤ 0.005).Red dots indicate genes known as contributing to insecticide resistance in malaria vectors.Black dots indicate genes affected by differential or outlier SNPs

Fig. 5
Fig. 5 Selection signatures observed in each selected line.SNPs diverging between the control line and each selected line were identified using a frequency-based approach (Diff SNPs) and a FST-based approach (outlier SNPs) and then averaged by gene (see methods).The upper Y axis shows the mean Diff SNP score per gene.The lower Y axis shows the proportion of outliers per gene.Symbol size increases with the number of polymorphic SNPs per gene.Triangles and circles denote candidate and non-candidate genes, respectively.Filled symbols indicate the presence of at least one differential or outlier SNP affecting the protein sequence.Blue and red symbols indicate candidate genes with a mean differential score > 0.4 or > 20% outliers, respectively; the corresponding gene names are indicated.Loci commonly associated with insecticide resistance in An. gambiae are indicated by dashed lines.The genomic scale shows chromosome arms with ticks every 10 Mb

Table 1
Composition of agrochemical mixtures used for