Insecticide resistance mechanisms associated with different environments in the malaria vector Anopheles gambiae: a case study in Tanzania

Background Resistance of mosquitoes to insecticides is a growing concern in Africa. Since only a few insecticides are used for public health and limited development of new molecules is expected in the next decade, maintaining the efficacy of control programmes mostly relies on resistance management strategies. Developing such strategies requires a deep understanding of factors influencing resistance together with characterizing the mechanisms involved. Among factors likely to influence insecticide resistance in mosquitoes, agriculture and urbanization have been implicated but rarely studied in detail. The present study aimed at comparing insecticide resistance levels and associated mechanisms across multiple Anopheles gambiae sensu lato populations from different environments. Methods Nine populations were sampled in three areas of Tanzania showing contrasting agriculture activity, urbanization and usage of insecticides for vector control. Insecticide resistance levels were measured in larvae and adults through bioassays with deltamethrin, DDT and bendiocarb. The distribution of An. gambiae sub-species and pyrethroid target-site mutations (kdr) were investigated using molecular assays. A microarray approach was used for identifying transcription level variations associated to different environments and insecticide resistance. Results Elevated resistance levels to deltamethrin and DDT were identified in agriculture and urban areas as compared to the susceptible strain Kisumu. A significant correlation was found between adult deltamethrin resistance and agriculture activity. The subspecies Anopheles arabiensis was predominant with only few An. gambiae sensu stricto identified in the urban area of Dar es Salaam. The L1014S kdr mutation was detected at elevated frequency in An gambiae s.s. in the urban area but remains sporadic in An. arabiensis specimens. Microarrays identified 416 transcripts differentially expressed in any area versus the susceptible reference strain and supported the impact of agriculture on resistance mechanisms with multiple genes encoding pesticide targets, detoxification enzymes and proteins linked to neurotransmitter activity affected. In contrast, resistance mechanisms found in the urban area appeared more specific and more related to the use of insecticides for vector control. Conclusions Overall, this study confirmed the role of the environment in shaping insecticide resistance in mosquitoes with a major impact of agriculture activities. Results are discussed in relation to resistance mechanisms and the optimization of resistance management strategies.


Background
Recent years have shown a decline in malaria prevalence attributed to efficient vector control strategies implemented in endemic areas [1]. In Tanzania, malaria is still a major public health problem, leading to high mortality and morbidity [2] and claiming more than one third of all deaths in children under five years old [3]. Members of the Anopheles gambiae complex are the main malaria vectors in Tanzania, with An. gambiae sensu stricto (referred to as An. gambiae s.s.) being endophagic and Anopheles gambiae arabiensis (referred to as An. arabiensis) being more exophagic [4]. Anopheles funestus is also present in southern Tanzania mainly during dry seasons [5].
As in other African countries, control programmes are based on the application of chemical insecticides by the use of insecticide-treated nets (ITNs) or indoor residual spraying (IRS). Current programmes are largely dependent on synthetic pyrethroids, which are the only WHO-recommended insecticides for ITNs [6]. DDT and bendiocarb are also used for IRS in several African countries, including Tanzania [7][8][9]. Both DDT and pyrethroids target the voltage-gated sodium channel (VGSC) in the mosquito central nervous system [10], while carbamates block the degradation of the neuromediator acetylcholine by inhibiting the acetylcholinesterase [11]. Following years of intensive usage, insecticide efficacy is now threatened by the rise of resistance in target populations. Such a phenomenon is occurring throughout the whole African continent and spreads at a rapid rate [12][13][14][15]. In Tanzania, reports of insecticide resistance in malaria vectors were rare but the situation is now changing with the rise of pyrethroid and DDT resistance in several regions [9,16,17]. Indeed, in 2011, a national wide survey indicated that most An. gambiae populations were still susceptible to carbamates and organophosphates while resistance to DDT and pyrethroids was rising in some areas [16].
In malaria vectors, resistance can be the consequence of mutations in the target proteins (target-site insensitivity), a lower penetration or sequestration of the insecticide, or an increased biodegradation of the insecticide due to enhanced detoxification activities (metabolic resistance). Resistance to pyrethroids appears to rely mainly on the 'knock down resistance' target-site mutations (kdr) and metabolic resistance mechanisms, although other mechanisms, such as cuticle alteration have also been evoked [18][19][20]. Kdr mutations and elevated levels of detoxification enzymes also confer resistance to DDT while carbamate resistance can be conferred by ace1 mutation and detoxification [21][22][23].
The recurrent use of insecticides for vector control is believed to be the main cause of resistance in mosquito populations [15,[24][25][26]. Concomitantly, the use of insecticides for the control of crop pests in agriculture and the presence of pollutants in urban and industrial areas have been suggested to play a significant role in selecting for insecticide resistance in mosquitoes [18]. Indeed, most pesticides used in agriculture have the same targets as those used for vector control, and can therefore select for resistance mechanisms in mosquitoes breeding in areas of intense agriculture activities [27][28][29][30]. Such cross-selections may also occur in urban environments where rapid urbanization promotes small-scale vegetable farming as source of livelihood and income and has led to an uncontrolled use of pesticides [31,32]. In addition, several studies showed that urban pollutants can affect mosquito detoxification system leading to an enhanced tolerance to insecticides [33][34][35][36][37][38]. In Tanzania, urbanization is spreading dramatically around large cities, and several areas show an intense agriculture activity [39] potentially affecting the resistance status of malaria vectors toward insecticides. Populations of An. gambiae from an area of intense agriculture in lower Moshi (Kilimanjaro region) were found highly resistant to pyrethroids [17]. As in other African countries, elevated insecticide resistance levels of malaria vectors have also been identified in Tanzanian urban areas [16,27,32,35].
In this context, the present study aimed at investigating insecticide resistance levels and associated mechanisms in malaria vectors from different environments in East Tanzania. To achieve this, nine Anopheles populations were sampled across three large areas characterized by different environments: i) an agriculture area that has a long history of pesticide use in coffee and rice farms and low ITN coverage, ii) an urban area from a large city with high anthropogenic pollutants and use of insecticides for vector control and iii) a low pollution area characterized by a limited use of pesticides for agriculture and vector control and a low urbanization. Resistance levels of each population at both larval and adult stages were determined by bioassays with deltamethrin, DDT and bendiocarb. The frequency of An. gambiae s.s. and An. arabiensis were assessed together with the presence of kdr mutations using PCR assays. A genome-wide microarray approach was used to compare transcriptome variations between each area and identify genes whose expression profiles were linked to insecticide resistance. These results are discussed in relation to known resistance mechanisms and their relationship with different environments.

Study sites and mosquito collection
This study took place in three large geographical areas in the east of Tanzania (Figure 1). These three areas were characterized by different environments as follows: i) urban area ii) agriculture area and iii) low pesticide area depending on urbanization, agriculture and vector control activities. Three populations separated by distances from 10 to 25 Km were studied in each area. Each population consisted of several breeding sites.
Urban area was characterized by a high degree of human activity, the presence of urban and industrial pollutants and a high ITN coverage. Dar es Salaam was chosen as the urban area because of its density (population~2.5 million) and its intense economic activity. Three populations were sampled in different city areas; Ilala (ILA), Temeke (TEM) and Kinondoni (KIN). Agriculture area was characterized by an intense agriculture activity with a recurrent use of pesticides on crops and a low ITN coverage. The lower Hai district, below Kilimanjaro Mount and near Moshi was chosen because of its low population density and intense agriculture activity requiring the use of pesticides (rice fields, sugar cane and coffee). Three populations were sampled across the area; Rundugai (RUN), Kikafu (KIK) and Kawaya (KAW). Low-pesticide area was characterized by low urbanization, small-scale farming with limited use of pesticides and a medium to high ITN coverage. The Muheza district, in the Tanga region was chosen as low pesticide area. Three populations were sampled across the area; Zeneti (ZEN), Kilulu (KIL) and Muheza estate (MUH).
Various ecological parameters were estimated from each population (Additional file 1) including breeding site type, surrounding area type (5 km range), surrounding area ITN coverage (5 km range data based on community interview), use of pesticide in agriculture (based on interview of farmers and pesticide shop owners), urbanization, malaria prevalence [40] and vegetation coverage. The five last parameters were quantified as follows for pairwise correlation studies (Additional file 2): absent (1), very low (2), low (3), medium (4), high (5).
Mosquitoes were collected from each population in 2011 during the long and short rain seasons. Anopheles F 0 larvae were collected and brought back to the laboratory for rearing in standardized conditions (natural photoperiod, 30°C, larvae feed on Tetramin® fish food). F 0 larvae were reared to adulthood, allowed to mate and females were blood fed on rabbits to obtain F 1 eggs. In order to avoid transient effects due to environmental variations, only F 1 larvae and F 1 adults reared in standardized conditions were used for bioassays and molecular analyses. Mosquitoes used for molecular analyses were stored individually in RNAlater (Ambion) at -20°C.

Bioassays with insecticides
Larvae and adult bioassays were performed following WHO protocols. Larvae bioassays were carried out using the pyrethroid deltamethrin, the organochlorine DDT and the carbamate bendiocarb and third stage larvae of each population. Four replicates per concentration and 6 concentrations in the activity range of each insecticide were used (N = 50 per replicate per concentration). Larval mortality was recorded after 24 h exposure to each insecticide and LC 50 (lethal concentration required for killing 50% of larvae) calculated using a probit approach.
Adult bioassays were performed on 3-5 days-old females following WHO guidelines with plastic test tubes and insecticide impregnated papers at the following concentrations: 4% DDT, 0.05% deltamethrin, and 0.1% bendiocarb. Four replicates of 10-25 mosquitoes were exposed to each insecticide for 60 minutes. For deltamethrin and DDT, the number of knockdown (KD) mosquitoes was recorded after 5, 10, 20, 30, 45 and 60 minutes exposure and KDT 50 (time required for knocking down 50% of individuals) estimated. After 1 h exposure mosquitoes were provided with 10% glucose solution and allowed recovering for 24 h before recording mortality. For bendiocarb, mosquitoes were exposed for varying durations (5, 10, 20, 30, 45 or 60 min). For each exposure time, adults were allowed to recover as described above before measuring mortality. Controls consisted of larvae or adults from each population not exposed to any insecticide. Mortality was corrected using Abbot Formula when the mortality rate of controls was between 5-20%. Two fully susceptible reference strains were also tested in parallel: the Kisumu strain (An. gambiae s.s. S form isolated from west Kenya) and the Ifakara strain (An. arabiensis isolated from central Tanzania). For each insecticide, larval resistance ratios (RR 50 ) consisted of LC 50 obtained for each population divided by LC 50 obtained with the Kisumu strain. For adults, RR 50 consisted of KDT 50 of each population divided by the KDT 50 of the Kisumu strain. As bendiocarb KDT 50s could not be estimated because of high mortality, % mortality after 10 min and 24 h recovery was used for comparing populations.

Correlation between resistance levels and ecological factors
Correlation between environmental data (ITN coverage, agriculture intensity, urbanization intensity, malaria prevalence and vegetation coverage), insecticide resistance levels (larval LC 50 , larval LC 95 , adult KDT 50 , adult KDT 95 and adult % mortality to each insecticide), kdr mutation frequency and the frequency of An. gambiae s. s. and An. arabiensis specimens were investigated across the nine field populations. Pearson's correlation coefficients were calculated for all possible pairwise comparisons and then subjected to a t-test using the formula t = r × √(n-1) / √(1-r 2 ). P values were then adjusted using Holm's multiple testing correction. Correlations showing a r ≥ 0.8 or ≤ -0.8 and an adjusted P value ≤ 0.05 were considered significant.

Species identification and kdr genotyping
For each population, genomic DNA samples (gDNA) from 50 larvae and 50 adults were extracted individually. gDNA was extracted by grinding and heating the last larval segment (larvae) or mosquito legs (adults) at 65°C for 30 minutes in 100 μl Bender buffer (0.1 NaCl, 0.2 M sucrose, 0.1 M tris-HCL pH 7.5, 0.05 M EDTA pH 9.1, 0.5% SDS) according to the method described by Collins et al. [41]. These individual gDNA samples was used for discriminating between An. gambiae s.s. and An. arabiensis following the PCR-based method described by Scott et al. [42]. The same samples were used for estimating the frequency of the L1014F and L1014S kdr alleles using the TaqMan® PCR diagnostic assays described in Bass et al. [43] and a MX3005P qPCR system (Agilent).

Gene expression study with microarrays
Because An. arabiensis represented the majority of individuals and in order to avoid potential bias due to different species ratios across populations, only individuals identified as An. arabiensis were used for RNA extraction and subsequent gene expression analyses. Total RNA was extracted from pools of 25 three days-old adult females from each population and three replicates of 25 individuals from the Ifakara susceptible strain using the RNAqueous-4PCR kit (Ambion). Total RNAs were treated with DNaseI (Ambion) to remove any genomic DNA contaminant. The quantity and integrity of RNA was analysed using a NanoDrop ND1000 Spectrophotometer (NanoDrop Technologies) and a 2100 Bioanalyzer (Agilent). Hundred ng total RNA were used for cDNA synthesis and T7-RNA amplification with Cy3-CTP or Cy5-CTP labelling using the 'Two-colors Low Input Quick Amp' labelling kit (Agilent Technologies) following manufacturer's instructions. Labelled cRNA were purified and their Cy3/Cy5 specific activity and quantity measured using a Nanodrop ND1000 (NanoDrop Technologies) and a 2100 Bioanalyzer (Agilent). Cy3-or Cy5-labelled cRNAs obtained from 3 biological replicates of the susceptible Ifakara strain were pooled together in equal quantity to build a common reference for all hybridizations. cRNAs from each population were hybridized versus cRNA of the reference (300 ng/sample) using the microarray AGAM 15 K (Agilent) representing more than 15 K An. gambiae transcripts as described in Mitchell [44]. After hybridization, non-specific probes were washed off according to manufacturer's instructions and slides were scanned with an Agilent G2205B microarray scanner. Spot finding, signal quantification and Lowess normalization were performed automatically using the Agilent Feature Extraction software (Agilent). Two hybridizations where Cy3 and Cy5 dyes were swapped were performed for each population for a total of 18 hybridizations. Normalized data were loaded into Genespring GX version 12.5 (Agilent). Spots showing a signal-to-noise ratio > 2 for both colours were flagged as detected and only probes detected in 50% of hybridizations were considered for further analysis.
Gene expression analysis was first performed at the area level (3 populations per area: 'Urban' , ' Agriculture' and 'Low pesticide'). For each area, mean expression ratios compared to the reference strain Ifakara were calculated, log transformed and submitted to a one sample student's t-test against 0 (equal transcription compared to reference strain) followed by Benjamini and Hochberg multiple testing correction. Transcripts showing a mean transcription ratio > 3 fold in either direction and a corrected P value < 0.001 were considered significantly differentially transcribed. These 416 transcripts were then submitted to hierarchical clustering analysis based on their expression profile across all areas. Clustering was based on Euclidean distance between log 2 ratios versus reference in each area and complete linkage algorithm. Transcripts belonging to each main cluster were then submitted to a Gene Ontology (GO) term enrichment analysis using Bingo version 2.4 [45]. Enrichment was assessed by a hypergeometric test followed by Benjamini and Hochberg multiple testing correction and restricted to the GO category 'Molecular function'. In each cluster, GO terms showing a corrected P value < 0.05 were considered significantly enriched compared to those represented in all detected transcripts. Significantly enriched GO terms were then graphically represented as networks using Cytoscape version 2.8.2 (cytoscape.org).
The 27 differentially expressed transcripts encoding genes potentially involved in insecticide resistance (insecticide targets, cuticle proteins, detoxification enzymes, ABC transporters, red/ox enzyme…) were then submitted to a focused clustering analysis as described above.
As differential tolerance to deltamethrin was observed between the nine studied populations, correlations between the expression profile of particular transcripts and deltamethrin RR 50 or kdr mutation frequency were investigated within all detected genes across all populations. Transcripts showing transcription profile with a correlation coefficient > 0.85 or < -0.85 to deltamethrin RR 50 or kdr mutation frequency were retained and submitted to a focused clustering analysis as described above.

RT-qPCR validation
Reverse-transcription quantitative PCR (RT-qPCR) was used to confirm the expression profile of the cytochrome P450 monooxygenase CYP6P3 (AGAP002865) and the cuticle protein AGAP000987 across populations. Specific primers targeting each transcript were designed using Beacon Designer 7.02 and their specificity checked using vectorbase Blast function. Total RNA samples used for RT-qPCR were identical to those used for microarrays. Experimental procedures used for reverse transcription and real-time PCR are described in Poupardin [46]. Data analysis was performed according to the ΔΔ Ct method taking into account PCR efficiency [47] and using the housekeeping genes encoding L8 (AGAP005802) and S7 (AGAP010592) ribosomal proteins. Four replicates were performed per population and results were expressed as mean transcription ratio relative to the reference strain Ifakara.

Mosquito species distribution
A total of 547 Anopheles larvae and adults from the nine sampled populations were used for species identification by PCR [40]. Most individuals belong to An. gambiae complex with only a few samples identified as other Anopheles species in various populations (Table 1). An. arabiensis was predominant in most populations and An. gambiae s.s. was mostly found in the urban area of Dar Es Salaam representing up to 50% in the ILA population.

Insecticide resistance levels
Larval bioassays revealed low resistance levels to deltamethrin, DDT and bendiocarb compared to the susceptible Kisumu strain (Table 2). Deltamethrin resistance ranged from 3.8 to 15.4 fold. One could note that a resistance ratio (RR 50 ) of 3.8 was observed for the susceptible Ifakara strain suggesting that populations showing a RR 50 below 4 fold can be considered as susceptible. Populations from the urban and agricultural areas showed elevated resistance to deltamethrin with a mean RR 50 of 10.9 fold and 8.2 fold respectively. Regarding DDT, resistance was slightly elevated in the urban area (mean RR 50 3.1 fold) but remain low in other areas. Overall, bendiocarb resistance was low although the KIK population from the agriculture area showed a 3.5 fold increased resistance.
Adult bioassays confirmed the moderate resistance levels of populations from urban and agriculture areas to deltamethrin with mean RR 50 of 3.1 fold and 5.6 fold respectively and a mortality rate after 1 h exposure to WHO diagnostic dose between 84 and 100% ( Table 3). As for larvae, DDT resistance was slightly elevated in urban area (2.2 fold in TEM). All tested populations showed 100% mortality following 1 h exposure to bendiocarb WHO diagnostic dose. However, when considering shorter exposure times, most field populations showed a better survival than the two susceptible populations suggesting a slight increased resistance to this insecticide (Table 3).

Kdr mutation frequency
Although both L1014F (kdr west) and L1014S (kdr east) mutations were screened, only the L1014S kdr 'east' mutation was detected (Table 4). This kdr mutation was present at high frequency in An. gambiae s.s. while only one heterozygous kdr individual was detected in An. arabiensis in the KIL population. In consequence, kdr mutation frequency was moderate in urban populations containing An. gambiae s.s. (ILA and TEM) and almost absent in populations from other areas.

Correlation between resistance levels and ecological factors
The only significant correlation found between ecological factors and resistance data involved deltamethrin resistance and agriculture (Additional file 2) where adult mortality after 1 h exposure to deltamethrin was negatively  correlated to agriculture intensity (higher tolerance in agriculture area, r = -0.95 adjusted P value = 0.0018). Significant correlations were also found between deltamethrin resistance and DDT resistance (Additional file 2). As expected, a strong positive correlation was found between the frequency of the L1014S kdr mutation and the presence of An. gambiae s.s.
Transcription profiling using the Anopheles gambiae 15 k microarray By using the An. gambiae 15 K microarray and the susceptible strain Ifakara as reference, the expression profiles of An. arabiensis females from each population were compared. This analysis allowed detecting 8280 probes corresponding to 7904 distinct transcripts showing a high and reproducible signal-to-noise ratio across all populations. When grouping populations by area, a total of 416 transcripts were found significantly differentially transcribed in any area compared to the susceptible reference strain (FC > 3 and adjusted P value < 0.001). Differential expression was mainly found in the agriculture area where 375 transcripts were found differentially expressed while low pesticide and urban areas showed less variation (65 and 5 transcripts respectively). Within each area, transcription ratios were well balanced between over-and undertranscribed genes with most variations being within 20 fold. Clustering analysis based on expression profile across  the three areas confirmed the strong variations observed in the agriculture area compared to other environments ( Figure 2 and Additional file 3). Gene ontology term enrichment analysis performed on each main cluster revealed an enrichment of multiple GO terms. Cluster 2 representing genes over-transcribed in agriculture areas was enriched in GO terms linked to receptors, oxidoreductases including cytochrome P450 monooxygenases and transcription factors. Cluster 5, representing genes specifically under-transcribed in agriculture area was enriched in the GO term 'nucleic acid binding' linked to transcription factors. Clusters 6 to 9 representing genes under-transcribed in all areas were enriched in genes and GO terms related to sugar and protein catabolism, energy storage, and odorant binding proteins ('odorant binding'). Focusing on transcripts encoding proteins putatively involved in insecticide resistance (i.e. detoxification enzymes sensu lato, ABC transporters, insecticide targets, cuticle proteins) confirmed the strong response of populations from agriculture area with 19 transcripts being significantly overtranscribed compared to the susceptible strain (Figure 3). These include three P450s (CYP6P1, CYP9J4, CYP4C28), two cuticle proteins (AAGP001329 and AGAP000987), one esterase (AGAP006227), one UDPGT (AGAP006775), one sulfotransferase (AGAP005721) and one aldehyde oxidase (AGAP006775). The GABA and nicotinic acetylcholine receptors were also over-transcribed in mosquitoes from agriculture areas.
At the population level, looking for genes whose expression profiles were positively or negatively correlated to insecticide resistance levels in adults identified multiple transcripts correlated to pyrethroid selection pressure ( Figure 4). Several genes were correlated to deltamethrin resistance including the GABA receptor (AGAP006028) and the nicotinic acetylcholine receptor (AGAP008588), the 'still life' gene (AGAP006590) and two genes encoding cuticle proteins (AGAP006829 and AGAP008449). Although kdr mutations were absent in An. arabiensis specimens used for microarray study, correlation between gene expression profiles and kdr mutation frequencies recorded at the population level were investigated based on the assumption that such association may still identify genes associated with pyrethroid selection pressure. Indeed, height transcripts were correlated to kdr mutation frequency obtained from whole populations, including two known P450 candidates (CYP6P3 and CYP9J5), the cuticle gene AGAP000987 and the gene AGAP002667 coding for a protein homologous to a tumour-related protein. Finally, only two genes had their expression profile negatively correlated with deltamethrin resistance or kdr mutation, including the transcript AGAP004203 encoding for vitellogenin C. Cross-validation of expression profiles of the cuticle protein coding gene AGAP000987 and the cytochrome P450 CYP6P3 revealed a good correlation between microarray and RT-qPCR although microarray transcription ratios appear slightly over-estimated (Additional file 4). These data confirmed the marked over-transcription of these two genes in urban and agriculture areas compared to the reference strain.

Discussion
The present study aimed at investigating insecticide resistance levels and identifying associated mechanisms across multiple An. gambiae populations from different environments. The eastern Tanzania region appeared suitable such study since insecticide resistance levels are currently increasing in this region with a significant heterogeneity observed in populations from contrasted environments. Although the experimental design (nine populations across three areas of contrasted environments) was constrained by field limitations, results support the role of environment in shaping the distribution Figure 3 Candidate genes differentially transcribed in any area. Clustering analysis was performed on the 27 transcripts showing a significant differential expression in any area and potentially involved in insecticide resistance. Clustering of transcripts was based on differences between mean Log 2 ratios in each area versus reference using Euclidean distance and complete linkage algorithm. Colour scale (green to red) indicates fold change versus reference strain. Stars indicate a significant differential transcription compared to the susceptible reference strain Ifakara (P ≤ 0.001 and FC ≥ 3). and mechanisms of insecticide resistance of malaria vectors in Tanzania.

Species distribution and insecticide resistance
Most sampled mosquitoes were identified as An. arabiensis with An. gambiae s.s. only present in some urban populations of Dar es Salaam. Even though a decline of the traditional malaria vector An. gambiae s.s. is observed in Tanzania, residual malaria transmission remains due to its replacement by the more exophagic An. arabiensis throughout the country [48]. The present study confirmed the preference of An. gambiae s.s. for permanent urban breeding sites [49], such as the Jangwani Valley in Dar es salaam and its potency to breed in organically polluted sites while An. arabiensis was mainly represented in temporary rural breeding sites [50] including those from intensive agriculture areas [17].
Overall, moderate DDT resistance was found in larvae and adults in populations from the urban area of Dar Es Salaam while deltamethrin resistance was highest in the agriculture area of the Kilimanjaro region. Deltamethrin resistance ratios were higher in larvae than adults although resistance profile across areas was conserved between life stages. The link between agriculture and pyrethroid resistance was confirmed by a significant correlation between deltamethrin resistance levels and agriculture intensity across all populations. Urban populations containing An. gambiae s.s. individuals carrying the L1014S kdr showed the highest resistance levels to DDT. Indeed, several studies confirmed the strong genetic association between DDT resistance and this kdr mutation [15,22,51]. In contrast, An. arabiensis populations found in the agriculture area showed higher resistance levels to deltamethrin together with an absence of kdr mutation frequencies, suggesting that other resistance mechanisms are predominant. Although all populations showed 100% mortality after 1 h exposure of adults to bendiocarb, 10 min exposure mortality data and larval bioassays suggest that some individuals from urban and agriculture areas may carry resistance alleles. This was expected as carbamates are commonly used in agriculture and are also available for domestic use as revealed by a survey through local dealers (Theresia Nkya unpublished). In addition, it is likely that some metabolic resistance mechanisms selected by Figure 4 Genes with expression profiles associated with deltamethrin resistance. Clustering analysis was performed on the 46 transcripts showing a correlation between their expression profile and deltamethrin resistance (RR 50 ) and kdr mutation frequency. Clustering of transcripts was based on differences between mean Log 2 ratios in each population versus reference using Euclidean distance and complete linkage algorithm. Colour scale (green to red) indicates fold change versus reference strain. + and -indicate a positive (r > 0.85) or negative (r < -0.85) correlation with deltamethrin resistance or kdr mutation frequency.
insecticides used for vector control are also conferring resistance to carbamates [52]. Such findings highlight the need for investigating background resistance and crossresistance mechanisms in mosquito populations before switching to carbamate insecticides in high resistance areas.
Transcription variations associated to different environments A 15 k An. gambiae microarray was used to investigate transcription level variations associated to different environments and insecticide resistance across sampled populations. Because of the high prevalence of An. arabiensis in our samples and potential gene expression variations between the two sibling species, only individuals identified as arabiensis were used for comparison with the susceptible An. arabiensis strain Ifakara. In order to focus on inherited variations and minimize transient variations linked to different environments or sampling bias, bioassays and gene expression studies were performed on F 1 individuals (i.e. progeny of those collected from the field) bred in standard laboratory conditions.
Most transcription level variations compared to the susceptible strain were found in populations from the agriculture area suggesting that the intense use of agrochemicals represent a significant selection pressure for mosquito populations. Such hypothesis was confirmed by an enrichment of GO terms related to detoxification and receptors in genes over-transcribed in this area. Of particular interest is the strong over transcription of three known insecticide targets: the acetylcholinesterase (targeted by organophosphates and carbamates), the GABA receptor (targeted by cyclodienes) and the nicotinic acetylcholine receptor (targeted by neonicotinoids). Increased gene copy number of insecticide receptors has frequently been associated to resistance and often lead to their over-expression in resistant populations [53][54][55][56][57]. The fact that organophosphates and cyclodienes are only used for crop protection in this area supports previous studies highlighting the potential of agricultural pesticides to select for resistance in mosquitoes [18]. Multiple detoxification enzymes were over-transcribed in the agriculture area including the alpha esterase AGAP006227, the UDPGT AGAP006775 and the cytochrome P450s CYP9J4 and CYP6P1. The Aedes aegypti paralog of this esterase (AAEL010389) has been found over-transcribed in pyrethroid resistant strains [58] and in response to xenobiotic exposure [34]. Similarly, two Ae. aegypti paralogs of this UDPGT (AAEL008560 and AAEL014371) were found over-transcribed in insecticide resistant strains [38,59]. Regarding CYP genes, although the role of CYP9Js remains unclear in Anopheles, this subfamily has been frequently involved in insecticide resistance in other insects with several of them validated as pyrethroid metabolizers in Ae. aegypti and Drosophila melanogaster [60][61][62]. CYP6Ps have also been frequently associated with insecticide resistance with An. gambiae CYP6P3 and An. funestus CYP6P9 being validated as pyrethroid metabolizers [36,60,63,64]. The over transcription of multiple detoxification enzymes in the agricultural area supports previous findings suggesting that in absence of kdr mutation, resistance of An. gambiae to pyrethroids in this region mainly relies on metabolic mechanisms favored by agricultural practices [17]. The gene AGAP000987 encoding the cuticle protein CPAP3-A1b was also strongly over-transcribed in the agriculture area. Cuticle modifications are likely involved in the resistance of mosquitoes to insecticides and genes encoding cuticle proteins have been frequently found over-transcribed in pyrethroid resistant populations [18][19][20]. An over-transcription of this gene has recently been reported in An. gambiae populations from West Africa strongly resistant to deltamethrin supporting its role in resistance (Hyacinthe Toe, unpublished data).
Although resistance to DDT and to a lesser extent to deltamethrin occurred in the urban area, very few genes were found differentially transcribed in this area compared to the reference strain. This can be explained by a lower resistance of An. arabiensis compared to An. gambiae s.s. (carrying the L1014S kdr mutation) in this area which were not used for gene expression study. In addition, although several studies revealed that mosquitoes exposed to urban pollutants display an increased tolerance to insecticides through the induction of their detoxification system [18,34,37], such transient effect may not fully express in F 1 individuals used for the microarray study. Nevertheless, resistance to insecticides in the urban area may not be the sole consequence of the presence of the kdr mutation as few detoxification genes, such as the known pyrethroid metabolizer CYP6P3 showed their highest transcription level in urban populations and were correlated to deltamethrin resistance (see below).
Finally, a fraction of candidate genes identified from the agriculture area were also found over-transcribed in the low pesticide area despite a low agriculture activity and moderate usage of insecticides for vector control in this region. Such situation may be the consequence of the migration of resistant alleles to this region due to the relative proximity of the Kilimandjaro agricultural area. Whether these resistant alleles are stable in this area and maintained by the usage of insecticides for vector control or other factors requires further investigations.

Transcription variations associated to pyrethroid resistance
As bioassay data revealed significant variations in deltamethrin resistance between populations, correlation between gene expression and deltamethrin resistance profiles were investigated across all populations. This analysis identified 46 transcripts showing an expression profile highly correlated to variations of deltamethrin RR 50 or kdr mutation frequency (r > 0.85 or r < -0.85). As expected, most transcripts correlated to deltamethrin resistance or kdr mutation frequency were over-transcribed in the agriculture and urban areas respectively. The transcript encoding vitellogenin C protein (AGAP004203) was found negatively correlated to both factors. Vitellogenin is the major yolk protein of most insects including mosquitoes and its negative correlation to deltamethrin resistance may reflect resistance costs linked to energy reallocation and general stress response [65,66]. The transcription profiles of the cytochrome P450s CYP6P3 and CYP9J5 were positively correlated to kdr mutation frequency supporting their role in deltamethrin resistance in Tanzania. Indeed, CYP6P3 has been showed to metabolize both permethrin and deltamethrin [67]. Anopheles gambiae CYP9J5 has been previously found over-transcribed in pyrethroid resistant populations [35,68] and various Ae. aegypti CYP9Js were shown to metabolize pyrethroids [61,62]. Unexpectedly, CYP6M2 often presented as the Anopheles P450 mainly responsible for metabolic resistance to pyrethroids and also metabolizing DDT [44,60,69] was not found significantly upregulated in any area compared to the susceptible reference strain nor correlated to deltamethrin resistance or kdr mutation frequency. This supports the importance of other detoxification enzymes in pyrethroid resistance in East Africa. Multiple genes encoding proteins involved in synapse functioning also had their transcription profiles positively correlated to deltamethrin resistance. This includes the Still life gene (AGAP006590) encoding a protein that converts GTPase from inactive to active form and modulates synapse differentiation in Drosophila [70] and the other GTPase activating protein AGAP00394. Genes encoding the protein Sidestep (AGAP001674), the neurexin AGAP004066, the JNK interacting protein AGAP006021, the cyclic-nucleotide-gated calcium channel AGAP007008, the Munc-13 protein AGAP005816, the nicotinic acetylcholine receptor AGAP008588 and the GABA receptor AGAP006028, all linked to synapse functioning and neurotransmitter activity were also found correlated to deltamethrin resistance [71][72][73]. Taken together, these results suggest that elevated neurotransmitter activity can contribute to insecticide resistance in mosquitoes as previously suggested [74]. Finally, three genes encoding cuticle proteins had their expression profiles correlated to deltamethrin RR 50 variations (CPR59 AGAP006829, CPLCG5 AGAP008449 and CPAP3-A1b AGAP000987), supporting their possible role in pyrethroid resistance through altered insecticide penetration.

Conclusions
Insecticide resistance is well established in malaria vectors throughout Africa and represents a threat to malaria control. Previously circumscribed to West and Central Africa, resistance has spread all over Africa dramatically. In absence of new insecticides approved for vector control, developing efficient resistance management strategies is critical for the next decade. This study provides ecotoxicological and molecular data supporting the role of the environment in which the vectors are found in shaping insecticide resistance mechanisms. Both agriculture and urban areas are likely favouring the emergence of resistance to insecticides used for vector control as a consequence of different factors. In agriculture areas, the massive usage of pesticides from various chemical families appears to select for metabolic and cuticle resistance mechanisms to a wide range of insecticides including pyrethroids and their potential alternative carbamates. As insecticide classes used for vector control are frequently used earlier for crop protection, this may explain why resistance of malaria vectors to all insecticides develops so rapidly in Africa. In this concern, implementing molecules/formulations showing limited cross-resistance with pesticides used for crop protection is of high interest for resistance management. In urban areas, the high pyrethroid-impregnated nets coverage completed by uncontrolled indoor spraying with insecticides may strongly select for kdr mutations and metabolic resistance mechanisms with a potential role of pollutants in favouring the selection of particular detoxification enzymes. Further understanding of how environment modulates the selection and spread of insecticide resistance mechanisms will help reducing vector control failures and improving resistance management strategies.

Additional files
Additional file 1: Ecological characteristics of sampled population.
Additional file 2: Correlation matrix between environmental factors and insecticide resistance data. For each pairwise comparison, Pearson correlation (r) coefficient (above diagonal) and its associated adjusted P value ≤ are shown (below diagonal). Two factors were considered significantly correlated if r ≥ 0.8 or ≤ -0.8 and if adjusted P value ≤ 0.05. Source data are shown as a separate sheet. Additional file 4: Cross-validation of microarray data with by RT-qPCR. Gene transcription data are expressed as fold change versus Ifakara strain. CYP6P3 AGAP002865 data are shown as plain marks. Cuticle protein AGAP000987 data are shown as empty marks. Squares, triangles and circles represent populations from urban, agricultural and low pesticide usage areas respectively. Population names and correlation coefficient between microarray and RT-qPCR fold changes are indicated.
Abbreviations ITN: Insecticide-treated net; IRS: Indoor residual spraying; DDT: 4,4'-(2,2,2trichloroethane-1,1-diyl)bis(chlorobenzene); VGSC: Voltagegated sodium channel; WHO: World Health Organization; KDT: Knock down time; KD: Knock down; LC: Lethal concentration; RR: Resistance ratio; The present study and T.N. research fellowship were funded by the EU FP7 project 'AvecNet' African Vector Control New Tools (EU grant agreement number 265660). We also acknowledge additionnal funding from the federative structure Environnemental and Systems Biology (BEeSy) of Grenoble -Alpes University. We thank Prof. H. Ranson for valuable comments on the manuscript. We are also grateful to Dr. P. Müller for help with statistical analyses and microarray design. We also thank Mr. Omary Mpili for mosquito sampling, rearing and field work. For project management and administration we thank Dr. Eve Worrall and Dr. Denis Massue.

Data deposition
Microarray data have been deposited to ArrayExpress under accession number E-MEXP-3987.