Additional Selection for Insecticide Resistance in Urban Malaria Vectors: DDT Resistance in Anopheles arabiensis from Bobo-Dioulasso, Burkina Faso

In the city of Bobo-Dioulasso in Burkina Faso, Anopheles arabiensis has superseded Anopheles gambiae s.s. as the major malaria vector and the larvae are found in highly polluted habitats normally considered unsuitable for Anopheles mosquitoes. Here we show that An. gambiae s.l. adults emerging from a highly polluted site in the city centre (Dioulassoba) have a high prevalence of DDT resistance (percentage mortality after exposure to diagnostic dose = 65.8% in the dry season and 70.4% in the rainy season, respectively). An investigation into the mechanisms responsible found an unexpectedly high frequency of the 1014S kdr mutation (allele frequency = 0.4), which is found at very low frequencies in An. arabiensis in the surrounding rural areas, and an increase in transcript levels of several detoxification genes, notably from the glutathione transferase and cytochrome P450 gene families. A number of ABC transporter genes were also expressed at elevated levels in the DDT resistant An. arabiensis. Unplanned urbanisation provides numerous breeding grounds for mosquitoes. The finding that Anopheles mosquitoes adapted to these urban breeding sites have a high prevalence of insecticide resistance has important implications for our understanding of the selective forces responsible for the rapid spread of insecticide resistant populations of malaria vectors in Africa.


Introduction
Almost half of the human population in Sub-Saharan Africa are projected to live in urban areas by 2030 [1]. Although urbanisation is typically thought to reduce malaria transmission [2], the rapid and unplanned growth of towns and cities is generally associated with inferior housing, poor sanitation, and increased pollution, all of which could impact on the distribution and abundance of mosquito vectors [3,4]. An evidence-based approach to understanding the ecology of urban malaria vectors for better control is therefore required.
Anopheles arabiensis is one of five mosquito species known to transmit malaria in African urban settings [3]. This member of the Anopheles gambiae species complex has a wide geographical range extending across the continent and is associated with drier climates and extensive land clearance [5]. Deforestation and urbanisation have created an arid environment into which An. arabiensis has invaded and adapted [6] and this species is the dominant malaria vector in several African cities [7][8][9][10]. The invasion of An. arabiensis into these 'urban islands' may be attributable to local adaptation to atypical breeding sites such as polluted water pools or ditches [8,[11][12][13].
High coverage with insecticide-based interventions such as longlasting insecticide treated nets (LLINs) and indoor residual spraying (IRS) has made a major contribution to reducing the malaria burden in Africa [14]. Insecticide resistance represents a serious threat to this progress particularly in light of ambitious plans to further scale-up coverage with these interventions [15]. Pyrethroids remain the only compound licensed for use on LLINs while the limited arsenal of available insecticides has led to a resurgence in the use of DDT for IRS [16]. Reports of pyrethroid and DDT resistance in Anopheles are widespread throughout Africa [17][18][19]. Two mechanisms have been conclusively shown to cause insecticide resistance in mosquitoes: i) mutations in the target-site protein of the insecticide and ii) increased detoxification of the insecticide by metabolic enzymes (e.g. P450 monooxygenases, glutathione-S-transerferases (GSTs) and carboxyesterases). DDT and the pyrethroid insecticides target the voltage-gated sodium channels (VGSC) in the nerve membranes of insects. Cross-resistance to these compounds in Anopheles is associated with two variant knockdown resistance (otherwise known as kdr) mutations at position 1014 in the VGSC; a leucine to phenylalanine (L1014F) and a leucine to serine (L1014S) substitution [20,21]. These kdr alleles are widely distributed in An. gambiae s.l. throughout Africa and have been subject to strong positive selective sweeps [22,23]. Despite the presence of both 1014F and 1014S in West African An. arabiensis [24,25], allele frequencies are much lower than those reported in M-and S-molecular forms of An. gambiae s.s. [18,26].
Gene expression profiling using targeted microarrays has proven a powerful experimental tool for identifying candidate genes associated with insecticide resistance [27][28][29]. In An. gambiae, several detoxification genes have been identified that are overexpressed in multiple insecticide resistant populations and the ability of several candidate cytochrome P450s and glutathione transferases (GSTs) to metabolise DDT and/or pyrethroids confirmed [29][30][31][32]. Less attention has been paid to other putative resistance mechanisms but the development of robust wholegenome microarrays [33] permits a much broader approach to studying this area.
The city of Bobo-Dioulasso lies within a large cotton belt in the south-west of Burkina Faso. An. gambiae s.s. collected from central Bobo-Dioulasso previously displayed reduced susceptibility to pyrethroids and DDT [34]. The vector dynamics in this area have, however, changed over the past decade with An. arabiensis superseding An. gambiae s.s. as the main member of the An. gambiae complex [35]. In this paper, we describe a DDT resistant population of An. arabiensis from a polluted breeding site in Bobo-Dioulasso and explore the genetic basis of this resistance. The implications for understanding the selective pressures on insecticide resistance in urban malaria vectors are discussed.

Study Site and Mosquito Collections
Third and fourth stage Anopheles larvae were collected from Dioulassoba (11u1094299N; 4u1793599W), an urban district in the centre of Bobo-Dioulasso, the second largest city of Burkina Faso. This district is crossed by a permanent stream which is intermittently surrounded by vegetable cropping. The stream is polluted with domestic waste and organic refuge (e.g. empty containers and daily water-waste). Small puddles along the stream and the surrounding borders constitute permanent breeding sites for both Anopheles and Culex larvae. Anopheles larvae were collected using a standard dipping method in August 2009 and January 2010, during the rainy and dry season respectively. Larvae were reared to adults at the insectaries of IRSS/Centre Muraz and blood fed to generate the test population which was then used in bioassays for resistance to DDT.

Ethics Statement
Dioulassoba has been a sentinel site of IRSS/Centre Muraz research team since the 1980s and, therefore, no specific work permits were required to work in this area. The stream is not privately owned and the field study did not involve protected or endangered species in any way.

DDT Susceptibility Bioassays
Bioassays were conducted on adult mosquitoes emerging from larvae collected in the dry and wet seasons. Three to five day old adult females An. gambiae s.l. were exposed for one-hour to papers impregnated with DDT (4%) in World Health Organisation (WHO) susceptibility tests [36]. Mosquito mortality was recorded 24 hours later. Control assays were performed throughout the experiment with a minimum of 25 mosquitoes exposed to noninsecticide treated papers. Mosquitoes were left in holding tubes for a further 24 hours (48 hours post-exposure) to minimise the effect of induction on gene expression from insecticide exposure.
Any mosquitoes which died during this period were not considered for molecular analysis. Between one and two legs were removed from all survivors (exposed and non-exposed) and placed on silica gel. The remainder of the insect was preserved in RNAlater (Sigma) for microarray experiments.

Species Identification and kdr Genotyping
Genomic DNA (gDNA) was extracted from the legs of adult female An. gambiae s.l. by heating to 95uC for 30 minutes in 20 ml of 16 PCR buffer (Promega). This was subsequently used to distinguish between An. arabiensis and An. gambiae s.s. following the method described by Scott et al. [37].
The frequency of the 1014F and 1014S kdr alleles in An. arabiensis exposed and non-exposed to DDT was determined using the TaqManH PCR assays described in Bass et al. [38]. TaqManH PCR reactions were performed on the Agilent MX3005P qPCR system. Hardy-Weinberg equilibrium analysis was performed following Rodriguez et al. [39].

Microarray Experiments
A whole-genome microarray approach was used to determine whether DDT resistance in An. arabiensis is associated with differences in gene expression. The transcriptional profiles of three experimental treatments were compared against each other ( Figure 1): (i) DSSBA RES -An. arabiensis from Dioulassoba surviving one-hour exposure to DDT (4%), (ii) DSSBA CON -An. arabiensis from Dioulassoba exposed to insecticide-free control papers and (iii) MOZ SUS -a laboratory colony of An. arabiensis originating from Mozambique. MOZ SUS is maintained at the Liverpool School of Tropical Medicine, is fully susceptible to all WHO-approved insecticides and was not exposed to DDT as part of this study.
All mosquitoes used in the microarray were adult females, five to seven day-old (i.e preserved 48 hours post bioassay) identified as An. arabiensis by PCR. Total RNA was extracted from six replicate pools of eight individuals from the three adult An. arabiensis samples (DSSBA RES, DSSBA CON, MOZ SUS) using TRIzolH reagent (Invitrogen) following the manufacturer's instructions. RNA was treated with two units of TURBO DNase TM (Ambion) to remove any remaining genomic DNA. The quantity and integrity of RNA was analysed using a NanoDrop (ND1000) Spectrophotometer (NanoDrop Technologies) and 2100 Bioanalyzer (Agilent technologies, Santa Clara, CA, USA). Two independent extractions were pooled to generate single biological replicates (a total of 16 mosquitoes per replicate). In total, three biological replicates per strain were used in the microarrays. cRNA was amplified from a starting template of 100 ng of total RNA and labeled with Cy3 and Cy5 fluorescent dyes using the two-color Low Input Quick Amp Labeling Kit (Agilent Technologies) according to the manufacturer's instructions. Samples were column purified (Qiagen) and the Cy3/Cy5 specific activity and cRNA quantity measured using a spectrophotometer (NanoDrop Technologies) and Bioanalyzer (Agilent Technologies).
Labeled cRNAs were hybridized to the Agilent 8615 k 'Anopheles gambiae' array (AGAM_15K; ArrayExpress website. Available: http://www.ebi.ac.uk/arrayexpress. Accessed 2012 Sep 4) designed at the Liverpool School of Tropical Medicine [33]. The An. gambiae chip contains eight replicated arrays spotted with 14,071 probes for the 12,604 An. gambiae s.s. genes identified in the Ensembl P3.5 annotation (September 2009). On each array, 281 detoxification genes which were previously included on the Anopheles gambiae detox chip [40], are represented by four separate probes. Microarray comparisons were made between the three biological replicates of each strain with dye-swaps. Hybridizations were conducted for 17 hours at 65uC at 10 rpm rotation and washed following the manufacturer's protocol (Agilent Technol-ogies). Scanning of each microarray slide was performed with the Agilent G2565AA/G2565BA Microarray Scanner System using the Agilent Feature Extraction Software (Agilent Technologies).
All microarray analyses were performed in Genespring GX 11.5 software (Agilent technologies). Transcripts flagged as 'detected' or 'non-detected' in at least five out of six hybridizations were filtered from the probe set and from these, only transcripts with the background subtracted signal greater than the background standard deviation in at least five out of six hybridizations were taken forward for statistical analysis. The transcription ratios were subjected to a one sample Student's t-test against the baseline value of 1 (equal gene expression in both strains) with a Benjamini-Hochberg FDR multiple testing correction. Transcripts with over a two-fold difference in expression and a t-test p-value below 0.05 after multiple testing were considered significantly different. The fold-changes of replicated probes were averaged for each representative gene. The microarray data have been deposited in ArrayExpress (ArrayExpress website. Available: http://www. ebi.ac.uk/arrayexpress/. Accessed 2012 Sep 4).
The transcript IDs of significantly up-or down-regulated genes were entered into the Functional Annotation tool available in the Database for Annotation Visualization and Integrated Discovery (DAVID) 6.7. DAVID systematically maps large lists of genes according to the associated biological annotation and calculates statistically significant enriched biological annotations from the target data set [41]. The most enriched biological terms were identified from each individual microarray comparison and presented according to their enrichment p-value (p,0.05) following Benjamini-Hochberg multiple testing.

RT-qPCR Validation of Candidate Genes
Reverse-transcription quantitative PCR (RT-qPCR) was used to confirm the expression of two target genes selected from the microarray experiments; a glutathione-S-transferase, GSTD3 (AGAP004382-RA) and an ABC transporter, ABCG4 (AGAP001333-RA). Primers spanning exon-exon boundaries (to prevent gDNA amplification) were designed using PrimerBLAST (NCBI) to minimise off-target amplification. First-strand cDNA was reverse-transcribed from approximately 5 mg of total RNA using 200 units/ml of SuperScript TM III Reverse Transcriptase (Invitrogen) and 1 ml of oligo(dT) 20 (50 mM), in the presence of 1 ml of RNaseOUT TM (Invitrogen). Samples were run on the Mx3005P QPCR system (Agilent Technologies) using a thermal profile of 95uC for 3 minutes followed by 40 cycles of 95uC for 10 seconds and 60uC for 10 seconds. qPCR reactions (20 ml) composed of 5 ml of cDNA (1 ng/ml), 16 Brilliant III Ultra-Fast SYBRH QPCR Master Mix (Agilent Technologies) and 300 nM of forward and reverse primers. Disassociation curves were run after each qPCR reaction to check for non-specific PCR products. The PCR efficiency and dynamic range of each primer set was determined using calibration curves generated from running  qPCR reactions over five serial dilutions (five-fold) of input cDNA. Information on primers and qPCR efficiencies are given in Table  S1. The same three biological replicates used for microarrays were analysed per treatment group and all qPCR reactions were performed in triplicate. No-template controls (NTCs) were run throughout all experiments. The expression of these genes was quantified using the comparative C T method (otherwise known as the 2 2DDCt method) [42]. Each target gene was normalised against the mean of two reference genes; ribosomal protein S7 (AGAP010592) and actin-5C (AGAP000651).

DDT Resistance in An. arabiensis from Dioulassoba
To determine DDT resistance levels from the polluted site of Dioulassoba, a total of 313 and 262 female An. gambiae s.l. were exposed for one hour to DDT (4%) during the rainy (August-September, 2009) and dry season (January 2010), respectively. Mosquito mortality in the rainy season was 65.8% (N = 313; 95% CI = 60.3270.1%) and in the dry season was 70.4% (N = 262 95% CI = 64.4275.9%) with no significant difference observed (Fisher's

Genotyping of kdr
The frequencies of the kdr alleles 1014F and 1014S in An. arabiensis from Dioulassoba were determined from a sub-sample of mosquitoes surviving (N = 63) and non-exposed (N = 76) to DDT using TaqManH PCR ( Table 1). The 1014F mutation was not found in An. arabiensis. The frequency of 1014S was high and remains in HW-equilibrium (p.0.05) but did not differ between non-exposed (0.40) and DDT survivors (0.43) (x 2 p = 0.610). 1014F was at a high frequency (0.88) from a small number of An. gambiae s.s. tested (exposed, N = 16; unexposed, N = 1), while 1014S was not found.

Transcription Profiles in An. arabiensis from Dioulassoba
All mosquito samples included in the microarrays were preidentified as An. arabiensis. Gene expression profiles were compared between three groups of An. arabiensis; (i) mosquitoes from Dioulassoba surviving DDT exposure (DSSBA RES), (ii) mosquitoes from Dioulassoba unexposed to insecticide (DSSBA CON), (iii) a laboratory insecticide-susceptible strain of An. arabiensis (MOZ SUS) (Figure 1). This three-way experimental design enabled identification of genes differentially expressed in the DDT resistant proportion of the DSSBA population, by comparing to the global population with the same genetic and environmental background and, in addition, provided information on expression relative to a standard reference strain, important for cross population/cross country comparisons. As mosquitoes were stored 48 hours post DDT exposure, the induction effect of insecticide exposure was expected to be minimal and hence the primary objective was to compare constitutive changes in gene expression.
The number of significantly differentially expressed probes for each comparison is given in Figure 2. As expected, the greatest number of differentially expressed genes was observed between DSSBA RES/DSSBA CON and MOZ SUS. A total of 360 probes (2.6% of the total probes on the array) were differentially expressed between DSSBA RES and DSSBA CON. The full list of genes differentially expressed in each comparison is given in the Supporting Information (Tables S2, S3, S4).
Enrichment analysis was initially performed on probes commonly expressed in both DSSBA RES and DSSBA CON compared with MOZ SUS (N = 299). This revealed gene-terms consistent with increased monooxygenase activity ( Table 2), suggesting constitutive differences in P450 activity between the resistant field mosquitoes and the susceptible laboratory strain. Thirty probes representing 10 individual P450s were significantly differentially expressed in DSSBA RES and DSSBA CON compared to MOZ SUS. Of these, six were commonly upregulated (CYP9J5, CYP6Z2, CYP6Z3, CYP6M3, CYP4H24, CYP12F1) whereas only members of the CYP325 class were down-regulated (Table 3). Moreover, a minimum of three probes were consistently up-regulated, representing CYP9M1 (2.62-fold) and CYP4C36 (3.66-fold), in DSSBA RES compared with DSSBA CON.
When the enrichment analysis was performed on probe sets differentially expressed between DSSBA RES and DSSBA CON, only gene-terms associated with glutathione-S-transferase (GST) activity were overrepresented and these possessed a high foldenrichment score (.10) ( Table 2). Thirteen probes representing seven GSTs were differentially expressed between DSSBA RES  and DSSBA CON. Six of these GSTs were up-regulated in DSSBA RES: GSTE6, GSTU4, GSTD1, GSTD3, GSTD11 and GSTE5 (Table 4). Of these, GSTD3 (AGAP004382) was significantly over-expressed in all microarray comparisons with a foldchange (average of three over-expressed probes) in the following descending order: DSSBA RES vs. MOZ SUS (foldchange = 5.46), DSSBA CON vs. MOZ SUS (fold-change = 2.57) and DSSBA RES vs. DSSBA CON (fold-change = 2.28) ( Table 5). Three genes from the ABC transporter family are expressed at higher levels in the DSSBA mosquitoes compared to MOZ SUS. All of these are members of the G family of half transporters which include the eye pigment precursors. One ABCG gene, ABCG4 (AGAP001333) is found at 2-fold higher levels in DSSBA RES than DSSBA CON. The closest ortholog to this gene in humans is ABCG2, which has been associated with resistance to cancer drugs [43]. Expression of ABCB1 is also increased in the DDT resistant subset relative to the non-exposed group. The ABCB family in vertebrates is also associated with multidrug resistance and contains the MDRR1 transporter [44]. A list of all ABC transporters expressed in microarray experiments are presented in Table S5.
Real-time qPCR confirmed the over-expression of GSTD3 in each comparison (Figure 3), although, as reported previously for the Agilent array platform, expression ratios obtained from RT-qPCR were on the whole higher than those generated by the microarrays [28,45]. RT-qPCR also confirmed the over expression of ABCG4 (AGAP001333) in DSSBA RES compared with DSSBA CON (2 2DDCt fold change = 4.97).

Discussion
Anopheles breeding sites are not restricted to clearly defined habitats in urban ecosystems [8]. Members of the An. gambiae species complex are increasingly found ovipositing and breeding in atypical polluted or domestic water bodies [12]. Our study site (Dioulassoba) represents one such area where mosquitoes breed within and alongside a permanent stream heavily polluted with organic and domestic waste. An. arabiensis comprised 0.79 (rainy season) and 0.99 (dry season) of An. gambiae s.l. collected from Dioulassoba, confirming previous findings which show that this  species has become the predominant member of the An. gambiae complex in this site over the past decade [35]. A similar pattern is emerging in Burkina Faso's capital Ouagadougou [9] and other areas of West Africa [46]. Insecticide resistance is an emerging threat for LLIN and IRS efficacy in rural areas but much less attention has focussed on the urban setting. Although control of urban malaria vectors may involve targeting the immature stages via measures such as larviciding and environmental management, prevention of transmission via the use of LLINs is usually the primary control measure. Unplanned development of towns and cities often leads to an increase in the number of available breeding sites, many of which will contain higher levels of pollutants or other xenobiotics than the typically more pristine rural sites. Mosquito larvae exposed to anthropogenic chemicals have been shown to be better able to survive insecticide exposure in laboratory settings [47]. The use of pesticides in urban agriculture, and the exposure to petroleum products have both been identified as a potentially strong driving force behind the selection of resistance [48,49].
Resistance to DDT and permethrin was reported in Bobo-Dioulasso over 10 years ago but An. arabiensis only comprised 0-8.3% of the An. gambiae s.l. population at this time [34]. In bioassays conducted in 2009-2010 a high prevalence of DDT resistance (mortality of 65.8-70.4%) was detected in An. arabiensis in the city. At the same time, resistance to DDT was observed in An. arabiensis from four rural sentinel sites (mortality after diagnostic dose exposure ranging from 78.9-84.2% [50]). Thus resistance is slightly more prevalent in the urban setting. More striking, however, is the difference in the frequency of the 1014S kdr mutation between the urban and rural sites.
The 1014S allele frequency is 40-43% in An. arabiensis from Dioulassoba. Frequencies of only 6-11% for 1014S have been observed in An. arabiensis from the four rural sites in Burkina Faso in the same time period [50]. The 1014S mutation has been detected in An. arabiensis from Sudan [51], Uganda [52], Kenya [53] and Benin [54] but in each case, the reported frequencies are low. Although no dead mosquitoes were screened from the bioassays to establish a genotype-phenotype link, previous field data suggests that 1014S is associated with higher levels of DDT and lower levels of permethrin resistance [52,55,56]. In vitro expression work using modified Drosophila VGSCs in Xenopus oocytes has shown that the 1014S allele confers higher resistance to DDT than the 1014F variant [57]. Clearly, the high frequency of 1014S is indicative of a selective benefit for this mutation though whether this benefit is a consequence of improved survival to insecticide exposure remains unknown. Since the frequency of the 1014S allele was only 40% in the resistant population, this suggests that other resistant mechanisms are present.
We applied whole-genome microarrays to identify additional putative DDT resistance mechanisms in An. arabiensis. Although the microarray used here is designed for studying transcription in An. gambiae s.s., the high degree of homology between the sibling species permitted cross-species hybridisation [58,59]. In this study, we used a three-way comparison between field mosquitoes surviving (DSSBA RES) or non-exposed (DSSBA CON) to DDT with a laboratory susceptible strain (MOZ-SUS) to assess constitutive gene expression. Several probes were commonly transcribed in both DSSBA RES and DSSBA CON compared with MOZ-SUS (N = 299). The most over-represented biological terms from this gene list were consistent with P450 activity. It is possible that fast-evolving genes such as P450s have diverged sufficiently between An. gambiae s.s. and An. arabiensis to augment issues of cross hybridisation that are common to all microarray studies on large gene families and further qPCR analysis is needed to confirm the expression patterns observed. It is of interest that two of the six P450s up-regulated in An. arabiensis from Dioulassoba (CYP6Z2 and CYP12F1) were over-transcribed in the DDTresistant An. gambiae s.s. strain ZAN/U [40]. The over-transcription of P450 genes and their association with insecticide resistance has been well documented in insect pests of agriculture and human disease vectors [27,28,60].
In DSSBA RES, genes with putative glutathione transferase (GST) activity were enriched amongst the differentially expressed subset. Increased GST activity is known to confer DDT resistance in mosquitoes [61,62] by catalysing the removal of a chlorine group from the pesticide. Interestingly, GSTD3 was over-expressed in all three microarray comparisons (and confirmed by RT-qPCR). Delta class GSTs have been implicated in insecticide resistance [63] but their role has previously thought to be relatively minor compared with those from the epsilon class (for example GSTE2). Functional validation of the interaction between GSTD3 and DDT is required to understand the relevance of the overexpression of this gene in all microarray comparisons.
Finally, members of the ABC-transporter family (otherwise known as P-glycoproteins) were over-expressed in An. arabiensis from Dioulassoba. ABC-transporters pump foreign molecules out of insect cells using an ATP-dependent mechanism. ABCtransporters are increasingly implicated in pesticide (herbicide, fungicide and insecticide) resistance although the precise mechanism is not understood [64,65]. The over-expression of an ABCtransporter (ABCB4) through gene amplification has recently been reported from pyrethroid-resistant Aedes aegypti [66].
Insecticide resistance is well established in Anopheles vectors throughout Africa and presents a genuine threat to malaria control. Until recently, An. arabiensis was generally regarded as the more susceptible sibling species compared to sympatric populations of An. gambiae s.s. in Burkina Faso and nearby regions of West Africa -this is not the case in Bobo-Dioulasso. The precise factors driving the selection for DDT resistance in An. arabiensis from Dioulassoba remain to be elucidated but an increase in insecticidebased interventions (including personal sprays and coils) and exposure to organic and anthropogenic pollutants may be important. Understanding the selection pressures in urban environments is important given the rapid rate of urbanisation, and concomitant increasing role of urban malaria in the epidemiology of this disease.

Supporting Information
Table S1 Quantitative PCR (qPCR) primers for target genes selected from microarrays. (XLSX)