Introduction Malaria mortality has decreased substantially in sub-Saharan Africa over the last decade, attributed in part to a massive scale-up in insecticide-based vector control interventions [1]. As the only insecticide class approved for treatment of bednets (ITNs) and the most widely used for indoor residual spraying (IRS), pyrethroids are by far the most important class of insecticides for control of malaria vectors [2]. Unfortunately pyrethroid resistance is now widespread and increasing in the most important malaria-transmitting Anopheles species [3]–[5] and catastrophic consequences are predicted for disease control if major pyrethroid failure occurs [6]. With no entirely new insecticide classes for public health anticipated for several years [5], [6] preservation of pyrethroid efficacy is critically dependent upon strategies such as rotation or combination of pyrethroids with just three other insecticide classes, organochlorines, carbamates and organophosphates [6], [7]. In addition to logistical and financial issues, insecticide resistance management suffers from knowledge-gaps concerning mechanisms causing cross-resistance between available alternative insecticides, and more, generally how high-level resistance arises [8]. With strongly- and multiply-resistant phenotypes documented increasingly in populations of the major malaria vector Anopheles gambiae in West Africa [9]–[13] such information is urgently required. Of the four classes of conventional insecticide licensed by the World Health Organisation (WHO), pyrethroids and DDT (the only organochlorine) both target the same para-type voltage-gated sodium channel (VGSC). This creates an inherent vulnerability to cross-resistance via mutations in the VGSC target site gene [14]–[16], which are now widespread in An. gambiae [5]. In contrast, carbamates and organophosphates cause insect death by blocking synaptic neurotransmission via inhibition of acetylcholinesterase (AChE), encoded by the ACE-1 gene in An. gambiae. Consequently, target site mutations in the VGSC gene producing resistance to pyrethroids and DDT will not cause cross-resistance to carbamates and organophosphates. The carbamate bendiocarb is being used increasingly for IRS [17], [18], and has proved effective in malaria control programs across Africa targeting pyrethroid- or DDT-resistant An. gambiae [18]–[20]. A single nucleotide substitution of glycine to serine at codon position 119 (Torpedo nomenclature; G119S) in the ACE-1 gene, which causes a major conformational change in AChE, has arisen multiple times in culicid mosquitoes [21], [22], and is found in An. gambiae throughout West Africa [23]–[25]. The G119S mutation can produce carbamate or organophosphate resistance [26] but typically entails considerable fitness costs [27]–[30]. This is beneficial for resistance management because in the absence of carbamates or organophosphates, serine frequencies should fall rapidly [29], [31]. In Culex pipiens, duplications of ACE-1 create linked serine and glycine alleles, which, when combined with an unduplicated serine allele, creates highly insecticide resistant genotypes with near-full wild-type functionality, thus providing a mechanism that can compensate for fitness costs [28], [31]. Worryingly, duplication has also been found in An. gambiae [23] though the consequences of copy number variation for fitness in the presence or absence of insecticide are not yet known in Anopheles. Though far from complete, information is available for metabolic resistance mechanisms to pyrethroids and DDT in wild populations of An. gambiae [5], [6], [32]–[34]. Indeed, a specific P450 enzyme, CYP6M2, has been demonstrated to metabolize both of these insecticide classes, suggesting the potential to cause cross-resistance in An. gambiae [32], [35]. By contrast little is known about metabolic mechanisms of carbamate resistance in mosquitoes and, as a consequence, potential for mechanisms of cross-resistance are unknown. A particularly striking and potentially problematic example of insecticide resistance has been found in one of the two morphologically identical, but ecologically and genetically divergent molecular forms comprising the An. gambiae s.s. species pair (M molecular form, recently renamed as An. coluzzii [36]) in Tiassalé, southern Côte d'Ivoire. The Tiassalé population is resistant to all available insecticide classes, and displays extreme levels of resistance to pyrethroids and carbamates [11]. The VGSC 1014F (‘kdr’) and ACE-1 G119S mutations are both found in Tiassalé [11], [25]. Yet kdr shows little association with pyrethroid resistance in adult females in this population [11]. ACE-1 G119S is associated with both carbamate and organophosphate survivorship [11], but this mutation alone cannot fully explain the range of resistant phenotypes, suggesting that additional mechanisms must be involved. Here we apply whole genome microarrays, transgenic functional validation of candidates, insecticide synergist bioassays, target-site genotyping and copy number variant analysis to investigate the genetic basis of (1) extreme bendiocarb resistance and (2) cross-insecticide resistance in An. gambiae from Tiassalé. Our results indicate that bendiocarb resistance in Tiassalé is caused by a combination of target site gene mutation and duplication, and by specific P450 enzymes which produce resistance across other insecticide classes. Results Whole genome transcription analysis Our study involved two microarray experiments (hereafter referred to as Exp1 and Exp2), involving solely M molecular form An. gambiae (Table S1), to identify candidate genes involved in bendiocarb resistance (full microarray results for Exp1 and Exp2 are given in Table S2A). In Exp1 gene expression profiles of female mosquitoes from bendiocarb-susceptible laboratory strains (NGousso and Mali-NIH) and a bendiocarb-susceptible field population (Okyereko, Ghana), none of which were exposed to insecticide, were compared to those of Tiassalé females. Two Tiassalé groups were used: either without insecticide exposure (Figure 1A), or the survivors of bendiocarb exposure selecting for the 20% most resistant females in the population [11] (Figure 1B). We used a stringent filtering process to determine significant differential expression (detailed in the legend to Figure 1), which included criteria on both the probability and consistency of direction of differential expression, and also required a more extreme level of differential expression in the Tiassalé-selected than Tiassalé (unexposed) vs. susceptible comparisons. Inclusion of this third criterion enhanced the likelihood that genes exhibiting differential expression are associated with bendiocarb resistance, rather than implicated via indirect association with another insecticide. Moreover, the requirement for significance in comparisons involving both bendiocarb-exposed and unexposed Tiassalé samples (Figure 1A, B) negates the possibility that any differential expression identified was a result solely of induction of gene expression by insecticide exposure. 10.1371/journal.pgen.1004236.g001 Figure 1 Microarray experimental design. Arrows indicate pairwise comparisons with direction indicating an increasing level of bendiocarb resistance, which was used to predict the expected direction of differential gene expression (only solid arrows were used to determine significance). Coloured boxes indicate samples resistant to bendiocarb; the red box indicates the only bendiocarb-selected sample. In Exp2 (C) microarray probes were considered significantly differentially expressed in resistant samples if: (i) each sus vs. res comparisons showed a consistent direction of expression as predicted by arrow direction; and (ii) each sus vs. res comparison yielded corrected P Tiassalé unexposed>Kovié) was met qualitatively for all genes (Figure 3). 10.1371/journal.pgen.1004236.g003 Figure 3 qRT-PCR expression analysis of candidate genes. Bars show mean fold changes relative to the bendiocarb and organophosphate susceptible Okyereko population. Asterisks indicate significant over-expression. Expression differences between pairs of populations are significant where error bars do not overlap. N = 5 biological replicates except for Tia_sel (N = 3). Insecticide resistance phenotypes of CYP6 genes in Drosophila For functional validation via transgenic expression in D. melanogaster, we chose CYP6P3 and CYP6M2; both of which have been shown to metabolize pyrethroids [34], [35], and CYP6M2 also DDT [32]. The capacity of each gene to confer resistance to bendiocarb, to the class I and II pyrethroids permethrin and deltamethrin, respectively, and to DDT and was assessed by comparing survival of transgenic D. melanogaster, exhibiting ubiquitous expression of CYP6M2 or CYP6P3 (e.g. UAS-CYP6M2/ACT5C-GAL4 experimental class flies), to that of flies carrying the UAS-CYP6M2 or CYP6P3 responder, but lacking the ACT5C-GAL4 driver (e.g. UAS-CYP6M2/CyO control class flies). For CYP6M2 the relative expression level of the experimental flies was 4.0 and for CYP6P3 4.3 (Table S3). As indicated by elevated LC50 values (Figure S4), expression of either CYP6M2 or CYP6P3 produced pyrethroid resistant phenotypes, and CYP6M2 expression also induced significant DDT resistance (Table 1). Assays for CYP6P3 with DDT did not produce reproducible results (data not shown). Flies expressing the candidate genes exhibited greater survival across a narrow range of bendiocarb concentrations (Figure S4). However, at a discriminating dosage of 0.1 µg/vial [37] a resistance ratio of approximately seven was exhibited for CYP6M2/ACT5C: CYP6M2/CyO flies (Mann-Whitney, P = 0.0002; Figure 4) with a much smaller, but still significant, ratio of approximately 1.4 (Mann-Whitney, P = 0.019) for CYP6P3/ACT5C: CYP6P3/CyO flies. Caution is required in quantitative interpretation of the resistance levels generated, both because of the non-native genetic background and also ubiquitous expression of genes that may be expressed in a tissue-specific manner [38]. Nevertheless, the bioassays on transgenic Drosophila show that each P450s can confer resistance to more than one insecticide class. 10.1371/journal.pgen.1004236.g004 Figure 4 Survival of transgenic Drosophila expressing An. gambiae Cyp6M2 or CYP6P3 in the presence of bendiocarb. Boxes show interquartile ranges with median lines and whiskers (error bars) show 95th percentiles for test (Act5C driver) or control (CyO) lines following exposure to 0.1 µg bendiocarb. Note that whiskers and median lines coincident with interquartile limits are not visible. Individual points falling outside percentiles are marked as dots. Mann-Whitney tests: ***P 1 replicate probes. (PPTX) Click here for additional data file. Figure S3 Relationship between expression measured by qRT-PCR and microarrays for candidate genes. The overall correlation is r = 0.50 (P = 0.056). (PPTX) Click here for additional data file. Figure S4 Survival of transgenic D. melanogaster that express CYP6M2 or CYP6P3 in the presence of varied amounts of insecticides. Log-linear plots of insecticide concentration vs. survival are shown. Blue points show survival of transformed flies with the Act5C driver which exhibit ubiquitous expression; red points show CyO control class flies. Bars show SEM of percent survival. Owing to the sharp inflection for both bendiocarb plots the regression model could not be applied to either Act5C or CyO data. N = 5 for all insecticides and concentrations other than bendiocarb at 0.1 µg, for which N = 8 (see Fig. 4). The gap in the x-axis results from use of a log scale on which control vials (zero insecticide) have no value. (PPTX) Click here for additional data file. Figure S5 Distributions of Ace-1 copy number ratios. Estimates are calculated relative to the susceptible Kisumu strain in Tiassalé samples genotyped as G119S heterozygotes that survived or died in a bendiocarb bioassay following PBO pre-exposure. (PPTX) Click here for additional data file. Figure S6 Interwoven microarray experimental loop design used in Exp2 comparing field samples from Kovie (KOV) with Malanville (MAL) and Okyereko (OKY) Each pool, indicated by a circle, represents mRNA extracted from 10 female An. gambiae s.s. M form mosquitoes. Arrows indicate individual microarrays (N = 18 in total), with direction representing microarray cy dye labelling. (PPTX) Click here for additional data file. Table S1 Details of An. gambiae samples used in experiments 1 and 2 (.xlsx). (XLSX) Click here for additional data file. Table S2 Microarray results (.xlsx). Table S2A . Full microarray results from both experiments. Table S2B . Significant microarray probes from Exp1. Table S2C . DAVID functional annotation clustering for Exp1. Table S2D . Significantly microarray probes from Exp2. Table S2E . Microarray probes significant in both Exp1 & Exp2. (XLSX) Click here for additional data file. Table S3 qRT-PCR expression results for transformed Drosophila melanogaster. (DOCX) Click here for additional data file. Table S4 GLM testing factors effecting bioassay mortality. (DOCX) Click here for additional data file. Table S5 Resistance association of the G119S target site mutation. (DOCX) Click here for additional data file. Table S6 qRT-PCR primer details for copy number variant analysis. (XLSX) Click here for additional data file. Table S7 qRT-PCR primer details for gene expression analysis. (XLSX) Click here for additional data file.
