Skip to main content
Download PDF

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

Malaria Journal volume 23, Article number: 69 (2024) Cite this article

Abstract

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.

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.

Methods

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

Agrochemical mixtures

Seven commercial formulations of products frequently used in agriculture in Côte d’Ivoire (each containing single or multiple insecticides, herbicides and fungicides and their adjuvants) were chosen according to surveys of farmers in Côte d’Ivoire to identify the pesticides commonly used by farmers [40]. These formulated products were combined to constitute the three agrochemical mixtures used for selection (Table 1). The stock insecticide mixture solution contained 3.3 ml/L Legumax® (deltamethrin); 2.7 ml/L K-Optimal® (lambda cyhalothrin, acetamiprid); 1.5 ml/L Verso 480® (chlopyriphos ethyl) and 66.7 g/L Furadan® (carbamate). The stock non-insecticide mixture contained 5.3 ml/L Banko plus® (fungicides, chlorothalonil and carbendazim); 5.3 ml/L Glyphader® (herbicide, glyphosate); 5.3 ml/L Garil® (herbicide, amide pyrimidine). A stock solution mixing all insecticide and non-insecticide compounds was also prepared (Table 1).

Table 1 Composition of agrochemical mixtures used for larval selection
Full size table

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 non-insecticide 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 (LD20) in the different lines (see Table 2). The selection pressure was maintained around the LD20 throughout the selection process by adjusting the doses every five generations. The new LD20s were determined using PoloPlus software.

Table 2 Concentrations of agrochemical mixtures used for larval selection
Full size table

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 chi2 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. RNA-seq 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.12 geneset 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 \(\ge\) 20, number of N ≤ 5, alignment score \(\ge\) 90, mapping quality \(\ge\) 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.12 geneset 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.12 geneset 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[(%freqSelected) − (%freqcontrol)]/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 FST 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.

Results

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 non-ins and mix lines were less resistant at G30 (non-ins line 54.74% and mix line 47.36% mortality at G30).

Fig. 1
figure 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

Full size image

Target-site mutations

The evolution of kdr mutations affecting the voltage-gated 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 (chi2 test P < 0.05).

Fig. 2
figure 2

Evolution of the kdr L995F mutation frequency during the selection process. Coloured bars show the genotype frequencies, as assayed from 30 individuals. Generation 0 corresponds to the parental line. Blue: 995FF kdr (resistant) homozygotes; Orange: 995FL heterozygotes; Grey: 995LL (wildtype) homozygotes. Black line: 995F allele frequency

Full size image
Fig. 3
figure 3

Evolution of the Ace1 mutation frequency during the selection process. Coloured bars show the genotype frequencies, as assayed from 30 individuals. Generation 0 corresponds to the parental line. Blue: 119SS (resistant) homozygotes; Orange: 119GS heterozygotes; Grey: 119GG (wildtype) homozygotes. Black line: 119S allele frequency

Full size image

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 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. Over-transcribed 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 cholesterol-like transporters (Niemman-Pick type C2 proteins, NPC2) potentially capable of binding xenobiotics were over-transcribed in all selected lines.

Fig. 4
figure 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

Full size image

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 AGAP007028 and AGAP012261), two cuticular proteins (CPLCX3 and CPLCG3) and the sensory appendage protein SAP1. Other candidate genes were found in regions showing more specific selection signals such as GSTE8 (Ins line only, Chr 3R), CPR130 (Ins and Mix lines, Chr X) and ace1 (Non-ins line only, Chr 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.

Fig. 5
figure 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

Full size image

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 non-insecticide 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 agriculture-based 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 cross-resistance 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 multi-genic 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.

Availability of data and materials

RNA-seq sequence data reported in this study have been deposited to the European Nucleotide Archive (ENA; http://www.ebi.ac.uk/ena) under the accession numbers E-MTAB-12756.

Abbreviations

WHO:

World Health Organization

IRS:

Indoor residual spraying

ITN:

Long-lasting insecticidal nets

VGSC:

Voltage-gated sodium channel

kdr:

Knock-down resistance

ace1:

Acetylcholinesterases1

References

  1. Sinka ME, Bangs MJ, Manguin S, Rubio-Palis Y, Chareonviriyaphap T, Coetzee M, et al. A global map of dominant malaria vectors. Parasit Vectors. 2012;5:69.

    Article  PubMed  PubMed Central  Google Scholar 

  2. Wolie RZ, Koffi AA, Ahoua Alou LP, Sternberg ED, N’Nan-Alla O, Dahounto A, et al. Evaluation of the interaction between insecticide resistance-associated genes and malaria transmission in Anopheles gambiae sensu lato in central Côte d’Ivoire. Parasit Vectors. 2021;14:581.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. WHO. World malaria report 2022. Geneva: World Health Organization; 2022.

    Google Scholar 

  4. Wilson AL, Courtenay O, Kelly-Hope LA, Scott TW, Takken W, Torr SJ, et al. The importance of vector control for the control and elimination of vector-borne diseases. PLoS Negl Trop Dis. 2020;14: e0007831.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. van den Berg H, Velayudhan R, Yadav RS. Management of insecticides for use in disease vector control: lessons from six countries in Asia and the Middle East. PLoS Negl Trop Dis. 2021;15: e0009358.

    Article  PubMed  PubMed Central  Google Scholar 

  6. WHO. World malaria report 2019. Geneva: World Health Organization; 2019. https://www.who.int/news-room/fact-sheets/detail/malaria

  7. Guillet P, Chandre F, Mouchet J. L’utilisation des insecticides en santé publique: état et perspectives. Med Mal Infect. 1997;27:552–9.

    Article  Google Scholar 

  8. WHO. Global malaria programme. Eliminating malaria. Geneva: World Health Organization; 2015.

    Google Scholar 

  9. WHO. World malaria report 2020. Geneva: World Health Organization; 2020.

    Google Scholar 

  10. Hemingway J, Ranson H. Insecticide resistance in insect vectors of human disease. Annu Rev Entomol. 2000;45:371–91.

    Article  CAS  PubMed  Google Scholar 

  11. Hemingway J, Hawkes NJ, McCarroll L, Ranson H. The molecular basis of insecticide resistance in mosquitoes. Insect Biochem Mol Biol. 2004;34:653–65.

    Article  CAS  PubMed  Google Scholar 

  12. Hemingway J, Coleman M, Paton M, McCarroll L, Vaughan A, DeSilva D. Aldehyde oxidase is coamplified with the World’s most common Culex mosquito insecticide resistance-associated esterases. Insect Mol Biol. 2000;9:93–9.

    Article  CAS  PubMed  Google Scholar 

  13. Martinez-Torres D, Chandre F, Williamson MS, Darriet F, Bergé JB, Devonshire AL, et al. Molecular characterization of pyrethroid knockdown resistance (kdr) in the major malaria vector Anopheles gambiae s.s. Insect Mol Biol. 1998;7:179–84.

    Article  CAS  PubMed  Google Scholar 

  14. Silva APB, Santos JMM, Martins AJ. Mutations in the voltage-gated sodium channel gene of anophelines and their association with resistance to pyrethroids—a review. Parasit Vectors. 2014;7:450.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Berthomieu A, Marquine M, Raymond M. The unique mutation in ace-1 giving high insecticide resistance is easily detectable in mosquito vectors. Insect Mol Biol. 2004;13:1–7.

    Article  PubMed  Google Scholar 

  16. Bourguet D, Pasteur N, Bisset J, Raymond M. Determination of ace.1 genotypes in single mosquitoes: toward an ecumenical biochemical test. Pestic Biochem Physiol. 1996;55:122–8.

    Article  CAS  PubMed  Google Scholar 

  17. Assogba BS, Djogbénou LS, Milesi P, Berthomieu A, Perez J, Ayala D, et al. An ace-1 gene duplication resorbs the fitness cost associated with resistance in Anopheles gambiae, the main malaria mosquito. Sci Rep. 2015;5:14529.

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  18. Camara S, Koffi AA, Ahoua Alou LP, Koffi K, Kabran JPK, Koné A, et al. Mapping insecticide resistance in Anopheles gambiae (s.l.) from Côte d’Ivoire. Parasit Vectors. 2018;11:19.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Fodjo BK, Koudou BG, Tia E, Saric J, N’dri PB, Zoh MG, et al. Insecticides resistance status of An. gambiae in areas of varying agrochemical use in Côte D’Ivoire. BioMed Res Int. 2018;2018:2874160.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Miles A, Harding NJ, Bottà G, Clarkson CS, Antão T, Kozak K, et al. Genetic diversity of the African malaria vector Anopheles gambiae. Nature. 2017;552:96–100.

    Article  ADS  Google Scholar 

  21. Chouaïbou M, Kouadio FB, Tia E, Djogbenou L. First report of the East African kdr mutation in an Anopheles gambiae mosquito in Côte d’Ivoire. Wellcome Open Research. 2017;2:5–9.

    Article  Google Scholar 

  22. Verhaeghen K, Van Bortel W, Roelants P, Backeljau T, Coosemans M. Detection of the East and West African kdr mutation in Anopheles gambiae and Anopheles arabiensis from Uganda using a new assay based on FRET/Melt Curve analysis. Malar J. 2006;5:16.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Dabiré RK, Namountougou M, Diabaté A, Soma DD, Bado J, Toé HK, et al. Distribution and frequency of kdr mutations within Anopheles gambiae s.l. populations and first report of the Ace1G119S mutation in Anopheles arabiensis from Burkina Faso (West Africa). PLoS ONE. 2014;9: e101484.

    Article  ADS  PubMed  PubMed Central  Google Scholar 

  24. Keïta M, Sogoba N, Kané F, Traoré B, Zeukeng F, Coulibaly B, et al. Multiple resistance mechanisms to pyrethroids insecticides in Anopheles gambiae sensu lato population from Mali, West Africa. Africa J Infect Dis. 2021;223:S81-90.

    Article  PubMed  Google Scholar 

  25. Elanga-Ndille E, Nouage L, Ndo C, Binyang A, Assatse T, Nguiffo-Nguete D, et al. The g119s acetylcholinesterase (Ace-1) target site mutation confers carbamate resistance in the major malaria vector Anopheles gambiae from Cameroon: a challenge for the coming irs implementation. Genes (Basel). 2019;10:790.

    Article  CAS  PubMed  Google Scholar 

  26. Essandoh J, Yawson AE, Weetman D. Acetylcholinesterase (Ace-1) target site mutation 119S is strongly diagnostic of carbamate and organophosphate resistance in Anopheles gambiae s.s. and Anopheles coluzzii across southern Ghana. Malar J. 2013;12:404.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Osta MA, Rizk ZJ, Labbé P, Weill M, Knio K. Insecticide resistance to organophosphates in Culex pipiens complex from Lebanon. Parasit Vectors. 2012;5:132.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Zoh DD, Ahoua Alou LP, Toure M, Pennetier C, Camara S, Traore DF, et al. The current insecticide resistance status of Anopheles gambiae (s.l.) (Culicidae) in rural and urban areas of Bouaké, Côte d’Ivoire. Parasit Vectors. 2018;11:118.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Vontas J, Katsavou E, Mavridis K. Cytochrome P450-based metabolic insecticide resistance in Anopheles and Aedes mosquito vectors: muddying the waters. Pestic Biochem Physiol. 2020;170: 104666.

    Article  CAS  PubMed  Google Scholar 

  30. Strode C, Wondji CS, David JP, Hawkes NJ, Lumjuan N, Nelson DR, et al. Genomic analysis of detoxification genes in the mosquito Aedes aegypti. Insect Biochem Mol Biol. 2008;38:113–23.

    Article  CAS  PubMed  Google Scholar 

  31. Ranson H, Guessan RN, Lines J, Moiroux N, Nkuni Z, Corbel V. Pyrethroid resistance in African anopheline mosquitoes: what are the implications for malaria control ? Trends Parasitol. 2010;27:91–8.

    Article  PubMed  Google Scholar 

  32. Vannini L, Reed TW, Willis JH. Temporal and spatial expression of cuticular proteins of Anopheles gambiae implicated in insecticide resistance or differentiation of M/S incipient species. Parasit Vectors. 2014;7:24.

    Article  PubMed  PubMed Central  Google Scholar 

  33. Yahouédo GA, Chandre F, Rossignol M, Ginibre C, Balabanidou V, Mendez NGA, et al. Contributions of cuticle permeability and enzyme detoxification to pyrethroid resistance in the major malaria vector Anopheles gambiae. Sci Rep. 2017;7:11091.

    Article  ADS  PubMed  PubMed Central  Google Scholar 

  34. Huang Y, Guo Q, Sun X, Zhang C, Xu N, Xu Y, et al. Culex pipiens pallens cuticular protein CPLCG5 participates in pyrethroid resistance by forming a rigid matrix. Parasit Vectors. 2018;11:1–10.

    Article  PubMed  PubMed Central  Google Scholar 

  35. Riveron JM, Yunta C, Ibrahim SS, Djouaka R, Irving H, Menze BD, et al. A single mutation in the GSTe2 gene allows tracking of metabolically based insecticide resistance in a major malaria vector. Genome Biol. 2014;15:R27.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Ranson H, Edi CVA, Koudou BG, Jones CM, Weetman D. Multiple-insecticide resistance in Anopheles gambiae mosquitoes, Southern Côte d’Ivoire. Emerg Infect Dis. 2012;18:1508–11.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Chandor-Proust A, Bibby J, Régent-Kloeckner M, Roux J, Guittard-Crilat E, Poupardin R, et al. The central role of mosquito cytochrome P450 CYP6Zs in insecticide detoxification revealed by functional expression and structural modelling. Biochem J. 2013;455:75–85.

    Article  CAS  PubMed  Google Scholar 

  38. Czeher C, Labbo R, Arzika I, Duchemin J. Evidence of increasing Leu-Phe knockdown resistance mutation in Anopheles gambiae from Niger following a nationwide long-lasting insecticide-treated nets implementation. Malar J. 2008;7:189.

    Article  PubMed  PubMed Central  Google Scholar 

  39. Aïzoun N, Aïkpon R, Akogbéto M. Evidence of increasing L1014F kdr mutation frequency in Anopheles gambiae s.l. pyrethroid resistant following a nationwide distribution of LLINs by the Beninese National Malaria Control Programme. Asian Pac J Trop Biomed. 2014;4:239–43.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Chouaïbou MS, Fodjo BK, Fokou G, Allassane OF, Koudou BG, David JP, et al. Influence of the agrochemicals used for rice and vegetable cultivation on insecticide resistance in malaria vectors in southern Côte d’Ivoire. Malar J. 2016;15:426.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Hien AS, Soma DD, Hema O, Bayili B, Namountougou M, Gnankiné O, et al. Evidence that agricultural use of pesticides selects pyrethroid resistance within Anopheles gambiae s.l. populations from cotton growing areas in Burkina Faso, West Africa. PLoS ONE. 2017;12: e0173098.

    Article  PubMed  PubMed Central  Google Scholar 

  42. Yadouleton A, Martin T, Padonou G, Chandre F, Asidi A, Djogbenou L, et al. Cotton pest management practices and the selection of pyrethroid resistance in Anopheles gambiae population in Northern Benin. Parasit Vectors. 2011;4:60.

    Article  PubMed  PubMed Central  Google Scholar 

  43. Hawkins NJ, Bass C, Dixon A, Neve P. The evolutionary origins of pesticide resistance. Biol Rev Camb Philos Soc. 2019;94:135–55.

    Article  PubMed  Google Scholar 

  44. Nkya TE, Poupardin R, Laporte F, Akhouayri I, Mosha F, Magesa S, et al. Impact of agriculture on the selection of insecticide resistance in the malaria vector Anopheles gambiae: a multigenerational study in controlled conditions. Parasit Vectors. 2014;7:480.

    PubMed  PubMed Central  Google Scholar 

  45. Oumbouke WA, Pignatelli P, Barreaux AMG, Tia IZ, Koffi AA, Ahoua Alou LP, et al. Fine scale spatial investigation of multiple insecticide resistance and underlying target-site and metabolic mechanisms in Anopheles gambiae in central Côte d’Ivoire. Sci Rep. 2020;10:15066.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. WHO. Test procedures for insecticide resistance monitoring in malaria vector mosquitoes, 2nd edn; 2016. https://fctc.who.int/publications/i/item/9789241511575

  47. Abbott WS. A method of computing the effectiveness of an insecticide. J Am Mosq Control Assoc. 1987;3:302–3.

    CAS  PubMed  Google Scholar 

  48. Khumallambam D, Kshetrimayum P, Nandeibam S, Huidrom S. An efficient protocol for total DNA extraction from the members of order Zingiberales- suitable for diverse PCR based downstream applications. Springerplus. 2013;2:669.

    Article  Google Scholar 

  49. Nadler SG, Tritschler D, Haffar OK, Blake J, Bruce AG, Cleaveland JS. Differential expression and sequence-specific interaction of karyopherin α with nuclear localization sequences. J Biol Chem. 1997;272:4310–5.

    Article  CAS  PubMed  Google Scholar 

  50. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B (Methodological). 1995;57:289–300.

    MathSciNet  Google Scholar 

  51. Huang DW, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: Paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009;37:1–13.

    Article  PubMed  Google Scholar 

  52. Saeed AI, Sharov V, White J, Li J, Liang W, Bhagabati N, et al. TM4: A free, open-source system for microarray data management and analysis. Biotechniques. 2003;34:374–8.

    Article  CAS  PubMed  Google Scholar 

  53. Pearson KLIII. On lines and planes of closest fit to systems of points in space. The London, Edinburgh, and Dublin Philos Magaz J Sci. 1901;2:559–72.

    Article  Google Scholar 

  54. Dray S, Dufour AB. The ade4 package: implementing the duality diagram for ecologists. J Stat Softw. 2007;22:1–20.

    Article  Google Scholar 

  55. Foll M, Gaggiotti O. A genome-scan method to identify selected loci appropriate for both dominant and codominant markers: a Bayesian perspective. Genetics. 2008;180:977–93.

    Article  PubMed  PubMed Central  Google Scholar 

  56. Adolfi A, Poulton B, Anthousi A, Macilwee S, Ranson H, Lycett GJ. Functional genetic validation of key genes conferring insecticide resistance in the major African malaria vector, Anopheles gambiae. Proc Natl Acad Sci USA. 2019;116:25764–72.

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  57. Müller P, Warr E, Stevenson BJ, Pignatelli PM, Morgan JC, Steven A, et al. Field-caught permethrin-resistant Anopheles gambiae overexpress CYP6P3, a P450 that metabolises pyrethroids. PLoS Genet. 2008;4: e1000286.

    Article  PubMed  PubMed Central  Google Scholar 

  58. Wagah MG, Korlević P, Clarkson C, Miles A, Lawniczak MKN, Makunin A. Genetic variation at the Cyp6m2 putative insecticide resistance locus in Anopheles gambiae and Anopheles coluzzii. Malar J. 2021;20:234.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Zoh MG, Tutagata J, Fodjo BK, Mouhamadou CS, Sadia CG, McBeath J, et al. Exposure of Anopheles gambiae larvae to a sub-lethal dose of an agrochemical mixture induces tolerance to adulticides used in vector control management. Aquat Toxicol. 2022;248: 106181.

    Article  CAS  PubMed  Google Scholar 

  60. Wang Y, Cheng P, Jiao B, Song X, Wang H, Wang H, et al. Investigation of mosquito larval habitats and insecticide resistance in an area with a high incidence of mosquito-borne diseases in Jining, Shandong Province. PLoS ONE. 2020;15:1–14.

    Google Scholar 

  61. Poupardin R, Riaz MA, Vontas J, David JP, Reynaud S. Transcription profiling of eleven cytochrome p450s potentially involved in xenobiotic metabolism in the mosquito Aedes aegypti. Insect Mol Biol. 2010;19:185–93.

    Article  CAS  PubMed  Google Scholar 

  62. Bara J, Montgomery A, Muturi J. Sublethal effects of atrazine and glyphosate on life history traits of Aedes aegypti and Aedes albopictus (Diptera: Culicidae). Parasitol Res. 2014;113:2879–86.

    Article  PubMed  Google Scholar 

  63. Diallo B, Lee Y, Reimer L, Fondjo E, Patchoke S, Ng A, et al. Relationship between kdr mutation and resistance to pyrethroid and DDT insecticides in natural populations of Anopheles gambiae. J Med Entomol. 2008;2:260–6.

    Google Scholar 

  64. Perera YC, Garcia GP, Segura KV, Monroy BL, Sanchez IPR, Lenhart A, et al. Impact of deltamethrin selection on kdr mutations and insecticide detoxifying enzymes in Aedes aegypti from Mexico. Parasit Vectors. 2020;13:224.

    Article  Google Scholar 

  65. Auteri M, La Russa F, Blanda V, Torina A. Insecticide resistance associated with kdr mutations in Aedes albopictus:. Biomed Res Int. 2018;2018:3098575. an update on worldwide evidences

  66. Amer K, Saavedra-Rodriguez K, Black WC, Gray EM. Effect of selection for pyrethroid resistance on abiotic stress tolerance in Aedes aegypti from Merida, Yucatan, Mexico. Insects. 2021;12:124.

    Article  PubMed  PubMed Central  Google Scholar 

  67. Lynd A, Oruni A, Van’T Hof AE, Morgan JC, Naego LB, Pipini D, et al. Insecticide resistance in Anopheles gambiae from the northern Democratic Republic of Congo, with extreme knockdown resistance (kdr) mutation frequencies revealed by a new diagnostic assay. Malar J. 2018;17:412.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Bamou R, Sonhafouo-Chiana N, Mavridis K, Tchuinkam T, Wondji CS, Vontas J, et al. Status of insecticide resistance and its mechanisms in Anopheles gambiae and Anopheles coluzzii populations from forest settings in South Cameroon. Genes (Basel). 2019;10:741.

    Article  CAS  PubMed  Google Scholar 

  69. Narahashi T. Neuronal ion channels as the target sites. Pharmacol Toxicol. 1996;79:1–14.

    Article  CAS  PubMed  Google Scholar 

  70. Edi CV, Djogbénou L, Jenkins AM, Regna K, Muskavitch MAT, Poupardin R, et al. CYP6 P450 enzymes and ACE-1 duplication produce extreme and multiple insecticide resistance in the malaria mosquito Anopheles gambiae. PLoS Genet. 2014;10: e1004236.

    Article  PubMed  PubMed Central  Google Scholar 

  71. Djogbénou L, Noel V, Agnew P. Costs of insensitive acetylcholinesterase insecticide resistance for the malaria vector Anopheles gambiae homozygous for the G119S mutation. Malar J. 2010;9:12.

    Article  PubMed  PubMed Central  Google Scholar 

  72. Djouaka RF, Bakare AA, Coulibaly ON, Akogbeto MC, Ranson H, Hemingway J, et al. Expression of the cytochrome P450s, CYP6P3 and CYP6M2 are significantly elevated in multiple pyrethroid resistant populations of Anopheles gambiae s.s. from Southern Benin and Nigeria. BMC Genomics. 2008;9:538.

    Article  PubMed  PubMed Central  Google Scholar 

  73. Nardini L, Christian RN, Coetzer N, Koekemoer LL. DDT and pyrethroid resistance in Anopheles arabiensis from South Africa. Parasit Vectors. 2013;6:229.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Toé KH, N’Falé S, Dabiré RK, Ranson H, Jones CM. The recent escalation in strength of pyrethroid resistance in Anopheles coluzzi in West Africa is linked to increased expression of multiple gene families. BMC Genomics. 2015;16:146.

    Article  PubMed  PubMed Central  Google Scholar 

  75. Lumjuan N, McCarroll L, Prapanthadara LA, Hemingway J, Ranson H. Elevated activity of an Epsilon class glutathione transferase confers DDT resistance in the dengue vector, Aedes aegypti White star. Insect Biochem Mol Biol. 2005;35:861–71.

    Article  CAS  PubMed  Google Scholar 

  76. Lumjuan N, Rajatileka S, Changsom D, Wicheer J, Leelapat P, La-aied Prapanthadara L, et al. The role of the Aedes aegypti Epsilon glutathione transferases in conferring resistance to DDT and pyrethroid insecticides. Insect Biochem Mol Biol. 2011;41:203–9.

    Article  CAS  PubMed  Google Scholar 

  77. Atoyebi SM, Tchigossou GM, Akoton R, Riveron JM, Irving H, Weedall G, et al. Investigating the molecular basis of multiple insecticide resistance in a major malaria vector Anopheles funestus (sensu stricto) from Akaka-Remo, Ogun State, Nigeria. Parasit Vectors. 2020;13:423.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  78. Zhou Y, Fu WB, Si FL, Yan ZT, Zhang YJ, He QY, et al. UDP-glycosyltransferase genes and their association and mutations associated with pyrethroid resistance in Anopheles sinensis (Diptera: Culicidae). Malar J. 2019;18:62.

    Article  PubMed  PubMed Central  Google Scholar 

  79. Antonio-Nkondjio C, Poupardin R, Tene BF, Kopya E, Costantini C, Awono-Ambene P, et al. Investigation of mechanisms of bendiocarb resistance in Anopheles gambiae populations from the city of Yaoundé, Cameroon. Malar J. 2016;15:424.

    Article  PubMed  PubMed Central  Google Scholar 

  80. Helvecio E, Romão TP, de Carvalho-Leandro D, de Oliveira IF, Cavalcanti AEHD, Reimer L, et al. Polymorphisms in GSTE2 is associated with temephos resistance in Aedes aegypti. Pestic Biochem Physiol. 2020;165: 104464.

    Article  CAS  PubMed  Google Scholar 

  81. Kaplanoglu E, Chapman P, Scott IM, Donly C. Overexpression of a cytochrome P450 and a UDP-glycosyltransferase is associated with imidacloprid resistance in the Colorado potato beetle, Leptinotarsa decemlineata. Sci Rep. 2017;7:1762.

    Article  ADS  PubMed  PubMed Central  Google Scholar 

  82. Coleman M, Vontas JG, Hemingway J. Molecular characterization of the amplified aldehyde oxidase from insecticide resistant Culex quinquefasciatus. Eur J Biochem. 2002;269:768–79.

    Article  CAS  PubMed  Google Scholar 

  83. Wu C, Chakrabarty S, Jin M, Liu K, Xiao Y. Insect ATP-binding cassette (ABC) transporters: roles in xenobiotic detoxification and Bt insecticidal activity. Int J Mol Sci. 2019;20:2829.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Chen X-Y, Yang Y, Wang J-Q, Wu Z-X, Li J, Chen Z-S. Overexpression of ABCC1 confers drug resistance to betulin. Front Oncol. 2021;11: 640656.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  85. Pignatelli P, Ingham VA, Balabanidou V, Vontas J, Lycett G, Ranson H. The Anopheles gambiae ATP-binding cassette transporter family: phylogenetic analysis and tissue localization provide clues on function and role in insecticide resistance. Insect Mol Biol. 2018;27:110–22.

    Article  CAS  PubMed  Google Scholar 

  86. Epis S, Porretta D, Mastrantonio V, Comandatore F, Sassera D, Rossi P, et al. ABC transporters are involved in defense against permethrin insecticide in the malaria vector Anopheles stephensi. Parasit Vectors. 2014;7:349.

    Article  PubMed  PubMed Central  Google Scholar 

  87. Ingham VA, Anthousi A, Douris V, Harding NJ, Lycett G, Morris M, et al. A sensory appendage protein protects malaria vectors from pyrethroids. Nature. 2020;577:376–80.

    Article  CAS  PubMed  Google Scholar 

  88. Gong Y, Duvvuri M, Duncan MB, Liu J, Krise JP. Niemann-Pick C1 protein facilitates the efflux of the anticancer drug daunorubicin from cells according to a novel vesicle-mediated pathway. J Pharmacol Exp Ther. 2006;316:242–7.

    Article  CAS  PubMed  Google Scholar 

  89. Naren D, Wu J, Gong Y, Yan T, Wang K, Xu W, et al. Niemann-Pick disease type C1(NPC1) is involved in resistance against imatinib in the imatinib-resistant Ph+ acute lymphoblastic leukemia cell line SUP-B15/RI. Leuk Res. 2016;42:59–67.

    Article  CAS  PubMed  Google Scholar 

  90. Nicoli ER, Eisa NA, Cluzeau CVM, Wassif CA, Gray J, Burkert KR, et al. Defective cytochrome p450-catalysed drug metabolism in Niemann-Pick type C disease. PLoS ONE. 2016;11: e0152007.

    Article  PubMed  PubMed Central  Google Scholar 

  91. Abdollahi M, Ranjbar A, Shadnia S, Nikfar S, Rezaie A. Pesticides and oxidative stress: a review. Med Sci Monit. 2004;10:RA141-147.

    CAS  PubMed  Google Scholar 

  92. Champion CJ, Xu J. Redox state affects fecundity and insecticide susceptibility in Anopheles gambiae. Sci Rep. 2018;8:13054.

    Article  ADS  PubMed  PubMed Central  Google Scholar 

  93. Oliver SV, Brooke BD. The role of oxidative stress in the longevity and insecticide resistance phenotype of the major malaria vectors Anopheles arabiensis and Anopheles funestus. PLoS ONE. 2016;11: e0151049.

    Article  PubMed  PubMed Central  Google Scholar 

  94. Yang X, Xu Y, Yin Q, Zhang H, Yin H, Sun Y, et al. Physiological characterization of chitin synthase A responsible for the biosynthesis of cuticle chitin in Culex pipiens pallens (Diptera: Culicidae). Parasit Vectors. 2021;14:234.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  95. Wipf NC, Duchemin W, Kouadio FPA, Fodjo BK, Sadia CG, Mouhamadou CS, et al. Multi-insecticide resistant malaria vectors in the field remain susceptible to malathion, despite the presence of Ace1 point mutations. PLoS Genet. 2022;18: e1009963.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This work was supported by the Wellcome Trust grants [103995] and [222542].

Author information

Authors and Affiliations

  1. University of Nangui Abrogoua, Abidjan, Côte d’Ivoire

    Christabelle G. Sadia, France-Paraudie A. Kouadio, Benjamin G. Koudou & Beatrice A. Adepo-Gourene

  2. Centre Suisse de Recherches Scientifiques (CSRS), Abidjan, Côte d’Ivoire

    Christabelle G. Sadia, Behi K. Fodjo, France-Paraudie A. Kouadio, Sebastien K. Oyou, Benjamin G. Koudou & Chouaibou S. Mouahamadou

  3. Laboratoire d’Ecologie Alpine (LECA) UMR 5553, Univ. Grenoble-Alpes, Univ. Savoie Mont Blanc, CNRS, 38000, Grenoble, France

    Jean-Marc Bonneville, Marius G. Zoh, Stephane Reynaud & Jean-Philippe David

  4. Vector Control Product Evaluation Centre (VCPEC)/Institut Pierre Richet, Bouaké, Côte d’Ivoire

    Marius G. Zoh

Authors
  1. Christabelle G. Sadia
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Jean-Marc Bonneville
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Marius G. Zoh
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Behi K. Fodjo
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. France-Paraudie A. Kouadio
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Sebastien K. Oyou
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Benjamin G. Koudou
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Beatrice A. Adepo-Gourene
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Stephane Reynaud
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Jean-Philippe David
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Chouaibou S. Mouahamadou
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

CGS contributed to study design, conducted experiments, contributed to data analysis and drafted the manuscript., JMB, JPD and SR contributed to RNA-seq data analysis and contributed to manuscript writing and revision. MGZ contributed to data analysis and contributed to manuscript writing. BKF and FPAK contributed to conduct experiments and contributed to manuscript writing. SKO provided technical help for mosquito rearing and bioassays. BGK and BAAG contributed to study design and to manuscript writing. CSM conceived the study, coordinated the study, contributed to manuscript writing and funding acquisition.

Corresponding author

Correspondence to Christabelle G. Sadia.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Figure S1.

Venn diagrams. Venn diagrams showing the numbers of genes significantly differentially transcribed in each selected line as compared to the control line (FC ≥ 1.5-fold in either direction and corrected P value ≤ 0.005). The number of resistance candidate genes are shown within brackets.

Additional file 2: Figure S2.

GO terms enrichment analysis. Enrichment analyses were performed for each line using genes significantly over- and under-transcribed separately (test lists) as compared to all genes detected by RNA-seq (reference list). Enrichment analyses were performed with the functional annotation tool DAVID (http://david.abcc.ncifcrf.gov) on all terms belonging to the ‘molecular function’ GO family. Only genes showing a corrected Fisher’s test P value ≤ 0.05 and represented by at least 4 genes are shown.

Additional file 3: Figure S3.

Principal Component Analysis of polymorphism data. PCA was performed using the frequency of all bi-allelic SNPs identified in each replicate of all lines as compared to the reference genome. Only the three first PCA axes are shown, accounting for 96.3% of the total variance.

Additional file 4: Table S4.

Transcription data obtained for all genes detected by RNA-seq.

Additional file 5: Table S5.

Gene-level polymorphism data obtained from RNA-seq.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Sadia, C.G., Bonneville, JM., Zoh, M.G. et al. The impact of agrochemical pollutant mixtures on the selection of insecticide resistance in the malaria vector Anopheles gambiae: insights from experimental evolution and transcriptomics. Malar J 23, 69 (2024). https://doi.org/10.1186/s12936-023-04791-0

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12936-023-04791-0

Keywords

Malaria Journal

ISSN: 1475-2875

Contact us