Introduction Aedes mosquitoes have shown a remarkable ability to develop resistance to insecticides [1]. Today, resistance to DDT, organophosphates and pyrethroids is widespread in the major dengue vector, Aedes aegypti [2] and this resistance is negatively impacting on control efforts. For example, in the Caribbean, resistance to pyrethroids is reducing the efficacy of pyrethroid space spraying in La Martinique and organophosphate resistance in Cuba is reducing the duration of control obtained by larviciding [3], [4]. Resistance to pyrethroids is of particular concern as this class of insecticides is increasingly replacing organophosphates in space spraying (WHO, 2011) and curtains, impregnated with pyrethroids, have also shown initial promise in reducing dengue transmission [5]. Most studies on the molecular basis of pyrethroid resistance focus on target site mutations [6], [7], [8]. Amino acid substitutions in the voltage gated sodium channel cause a resistance phenotype to pyrethroid insecticides known as knockdown resistance or kdr. At least four amino acid substitutions in the sodium channel (I1011M, V1016G, V1016I and F1534C) have been linked to resistance in Ae. aegypti [9], [10], [11], [12]. Two of these alleles, 1016I and 1534C, are widely distributed in the Caribbean [9]. The role of other resistance mechanisms is less clearly understood. Biochemical assays are frequently used to screen for metabolic resistance caused by elevated activities of cytochrome P450s, carboxylesterases and/or glutathione transferases. Although these assays lack sensitivity they have provided preliminary evidence that metabolic resistance is involved in conferring pyrethroid resistance in Caribbean populations of Ae. aegypti from Cuba, Grand Cayman, La Martinique and Trinidad . To determine the molecular basis of this metabolic resistance and, identify other pathways potentially involved in conferring the resistance phenotype, we utilized microarray and quantitative PCR to analyse the near complete transcriptome of pyrethroid resistant populations from Cuba and Grand Cayman. The results confirm that elevated cytochrome P450 activity is strongly associated with pyrethroid resistance in these populations. Comparative analysis of the data from this, and earlier studies on populations from Latin America and Southeast Asia, indicates that the CYP9J family of P450 enzymes is primarily responsible for metabolic resistance to pyrethroids in Ae. aegypti. Materials and Methods Mosquito strains and sample collection Three strains of Ae. aegypti were used in this study. The NEW ORLEANS (NO) strain is a laboratory strain that is susceptible to all known insecticides and was originally colonized by the Center for Disease Control and Prevention (CDC) Atlanta, USA. The pyrethroid resistant CAYMAN strain was colonized from larvae collected in routine field surveillance sites in Grand Cayman in 2008 . This strain has very high levels of resistance to DDT (>90% survival after 8 hours exposure to 4% DDT) and pyrethroids (resistance ratio of 109-fold to permethrin and 30-fold to deltamethrin compared with the susceptible New Orleans strain [9]). The CUBA-DELTA SAN 12 strain (CUBA-DELTA) was collected in 1997 in Santiago de Cuba. It was selected for 12 generations at the larval stage with deltamethrin at the Institute ‘Pedro Kouri’ in Havana, Cuba. CUBA-DELTA larvae were highly resistant to this insecticide (>1000-fold) and this resistance was also manifested at the adult stage [13]. Egg papers from the CAYMAN strain and the CUBA-DELTA strain were sent to the Liverpool School of Tropical Medicine, UK and the mosquitoes were reared under standard laboratory conditions (26°C, 80% RH) and a 12∶12 hours light∶dark cycle. Detection of target site mutations The prevalence of the 1016I and 1534C kdr mutations in the CAYMAN strain has been reported previously. For the CUBA-DELTA strain, 38 mosquitoes were genotyped for the 1534C mutation using the tetraplex assay described in [9] and for the 1016I mutation using the hot oligonucleotide ligation assay (HOLA) [11]. RNA extractions and labeling of cRNA For each strain, total RNA was extracted from three pools of 30, three day old, non blood-fed females using Pico Pure™ RNA Isolation Kit (Applied biosystems, Foster city, CA, USA). The strains were reared in parallel to minimize variation resulting from breeding conditions. Each biological replicate consisted of mosquitoes from distinct generations to control for stochastic variations. The quality and concentration of RNA was assessed using a 2100 Bioanalyzer (Agilent technologies, Santa Clara, CA, USA). Then, 100 ng of total RNA were used for RNA amplification and labeled with Cy-3 and Cy-5 fluorescent dyes using the Two Colors Low Input Quick Amp Labeling Kit (Agilent technologies) according to manufacturer's instructions. Labeled cRNAs were purified with the Qiagen RNeasy spin columns (Qiagen, Hilden, Germany). Quantification and quality assessment of labeled cRNA were performed with the Nanodrop ND-1000 (Thermo Scientific, DE, USA) and the Agilent 2100 Bioanalyser (Agilent Technologies). Purified labeled cRNAs were stored at −80°C until microarray hybridizations. Hybridizations, data acquisition and statistical analysis Hybridizations were made to the ‘Liverpool Aedes aegypti Agilent 8×15K v1’ microarray (A-MEXP-1966) designed by the Liverpool School of Tropical Medicine. Each array contains 60mer oligo-probes representing >14320 Aedes aegypti transcripts (93% of the putative gene count, 79% of putative transcripts –the lower coverage of transcripts is a consequence of the multiple putative transcripts for some genes). Labeled cRNA from CAYMAN and CUBA-DELTA were co-hybridized with age-matched NO samples, in direct pairwise comparisons. For two out of the three biological replicates, dye swaps were performed making a total of five hybridisations per comparison. Labeled targets were hybridized to the array for 17 h at 65°C and 10 rpm rotation and then washed according to Agilent protocol. Slides were scanned on Agilent G2565AA/G2565BA Microarray Scanner System using Agilent Feature extraction software (Agilent technologies). Genespring GX 11.1 software (Agilent technologies) was used for normalization and statistical analysis. To account for multiple testing , p-values were adjusted adopting the approach of Benjamini and Honchberg [14] to control for the false positives. Transcripts showing an absolute fold change >2-fold in either direction and a t-test P-value lower than P 2-fold in either direction and a t-test P-value lower than P 2 fold in either direction) in one or both resistant strains compared to the susceptible New Orleans strain. Upward arrows indicate over- transcribed in resistant strains, downward represent under-transcribed. The predicted functions of the genes differentially expressed in both populations were identified by BLAST2GO. More than 43% of the differentially expressed genes in the Cuba vs New Orleans comparison are annotated as ‘conserved hypothetical proteins’ in Vectorbase and 49% for Cayman vs New Orleans. These are listed in tables S2 to S4 but are not discussed further in this manuscript. For further analysis we focused primarily on the subset of genes that were differentially expressed in both populations (Table 2), although other genes of interest are also discussed. 10.1371/journal.pntd.0001692.t002 Table 2 Accession numbers and putative functions for the 20 genes showing the highest elevation in expression in the CAYMAN vs NO comparison and the 20 genes showing the highest decrease in expression in the same comparison. Cayman vs NO Cuba vs NO Transcript ID Description Fold change Corrected p-value Fold change Corrected p-value AAEL009076-RA nadh dehydrogenase subunit 4 72.55 0.00113 43.79 0.00890 AAEL012836-RA cytochrome b561 25.03 0.00147 20.43 0.00237 AAEL015136-RA niemann-pick type c- 21.31 0.00252 8.30 0.00635 AAEL002813-RA coupling factor, putative 14.71 0.00069 9.74 0.00237 AAEL014617-RA cytochrome p450 CYP9J28 14.21 0.00091 12.79 0.00346 AAEL000385-RA developmentally regulated rna-binding protein 9.64 0.00455 4.62 0.00364 AAEL007083-RA protein zer-1 homolog 9.16 0.00132 2.57 0.00464 AAEL014893-RA cytochrome p450 CYP6BB2 8.45 0.00040 7.07 0.00133 AAEL014616-RA cytochrome p450 CYP9J27 8.41 0.00058 6.80 0.00042 AAEL001668-RA enolase 8.39 0.00098 10.78 0.00138 AAEL010227-RA dolichol-phosphate mannosyltransferase 8.20 0.00143 10.93 0.00586 AAEL003099-RA glucosyl glucuronosyl transferases 6.99 0.00750 7.52 0.00292 AAEL009798-RA zinc finger protein 6.67 0.00160 4.03 0.00416 AAEL007849-RA hypothetical protein AaeL_AAEL007849 [Aedes aegypti] 6.62 0.00164 5.28 0.00242 AAEL008390-RA guanylate cyclase 6.50 0.00132 4.98 0.00496 AAEL004943-RA riken cdna isoform cra_a 6.28 0.00149 8.51 0.00493 AAEL010761-RA grip and coiled-coil domain-containing protein 1 6.09 0.00082 4.55 0.00042 AAEL001390-RA hypothetical protein AaeL_AAEL001390 [Aedes aegypti] 6.08 0.00113 6.87 0.00153 AAEL014609-RA cytochrome p450 CYP9J26 5.69 0.00088 4.86 0.00586 AAEL014645-RA hypothetical conserved protein 5.59 0.00304 5.36 0.00237 AAEL012440-RA sodium-bile acid cotransporter −5.56 0.00393 −2.79 0.00320 AAEL003785-RA uncharacterized protein kiaa0090 −5.57 0.00200 −2.23 0.00518 AAEL008025-RA cg16787 cg16787-pa −5.70 0.00029 −2.93 0.00320 AAEL011597-RA pyridoxal phosphate phosphatase phospho2 −5.79 0.00108 −3.96 0.00990 AAEL009949-RA homeotic antennapedia −5.86 0.00051 −12.09 0.00252 AAEL006515-RA bcdin3 domain containing −6.76 0.00235 −8.02 0.00242 AAEL014128-RA hypothetical protein AaeL_AAEL014128 [Aedes aegypti] −8.04 0.00070 −3.91 0.00392 AAEL009462-RA hydroxyacylglutathione hydrolase −8.42 0.00049 −6.81 0.00042 AAEL009335-RA adhesion regulating molecule 1 (110 kda cell membrane glycoprotein) −8.60 0.00031 −3.18 0.00078 AAEL003485-RA adhesion regulating molecule 1 (110 kda cell membrane glycoprotein) −8.77 0.00031 −3.17 0.00129 AAEL013403-RA hypothetical protein AaeL_AAEL013403 [Aedes aegypti] −8.81 0.00070 −11.83 0.00042 AAEL005204-RA isoform a −10.21 0.00132 −7.68 0.00735 AAEL000219-RA lactoylglutathione lyase −11.40 0.00031 −5.47 0.00236 AAEL002812-RA lyr motif-containing protein 5 −13.57 0.00164 −2.60 0.00501 AAEL000721-RA deoxynucleotidyltransferase terminal-interacting −13.99 0.00031 −12.24 0.00352 AAEL007244-RA zinc finger protein −16.33 0.00266 −2.46 0.00462 AAEL008079-RA trypsin- −17.70 0.00082 −6.14 0.00320 AAEL001336-RA charged multivesicular body protein 2a −36.40 0.00528 −3.42 0.00680 AAEL002811-RA kda midgut protein −77.13 0.00031 −3.39 0.00523 Detoxification genes are shown in bold. Genes differentially expressed in both populations Enrichment analysis was used to identify particular GO terms that were over represented in the subset of transcripts up regulated in both resistant populations. Thirteen GO terms were significantly over represented in the up-regulated subset (Figure 2). However, after Benjamini and Hochberg multiple testing correction (Pval 2-fold in the CUBA-DELTA strain and CAYMAN strain relative to the susceptible NO population respectively (Table 1). Seven of these CYPs were up-regulated in both strains: CYP6BB2, CYP9J9, CYP9J10, CYP9J26, CYP9J27, CYP9J28 and CYP329B1. Three P450s were down regulated in the CAYMAN strain and two in the CUBA-DELTA strain but none of these are common to both strains. Twelve of the 24 up-regulated CYPs (and five of the seven up-regulated in both strains) belong to the CYP9J family. Further genes with roles in oxidative metabolism of xenobiotics were amongst the subset of genes over expressed in both resistant populations. This included a dimethylanaline monooxygenase (AAEL00834), a member of the cytochrome b561 family (AAEL012836) and subunit 4 of NADH dehydrogenase (AAEL009076). Other detoxification genes included the glutathione transferase, GSTe4 (AAEL007962), and two glucosyl glucornosyl transferases (AAEL003099 and AAEL014246) (Table 1). Seven P450s (plus an ABC gene described below) were selected to validate the microarray results by qPCR. These included six out of the seven P450s over expressed in both strains and an additional P450 that is only over expressed in the CAYMAN strain, CYP9J19. In general there is good agreement between the qPCR and microarray data (Table S5) with the exception of the CYP6BB2 gene, (AAEL014893). The high level of over expression of this gene observed in the microarray could not be confirmed by qPCR using two alternative primer sets. However, the qPCR confirmation of over expression of the CYP9J genes adds further support for these enzymes playing a role in resistance to pyrethroids in these Caribbean populations. Other genes of interest An ABC transporter gene, AAEL006717, was expressed at approximately 5-fold and 2-fold higher levels in the pyrethroid resistant populations from Cayman and Cuba respectively, relative to the susceptible NO strain. This gene is potentially of interest because elevated ABC transporters have been linked to insecticide resistance in several species [20], [21], [22], [23] although the physiological mechanism by which these transporter proteins act to reduce insecticide susceptibility is unknown. The over expression of AAEL006717, which is an orthologue of the An. gambiae ABCB4 gene, [24] was confirmed by qPCR for the CAYMAN strain but not attempted in the Cuban population (Figure 3 Table S5). 10.1371/journal.pntd.0001692.g003 Figure 3 Quantitative PCR analysis of selected genes from the microarray experiments. Relative-fold change in transcript and gene copy number normalised to two ribosomal genes was compared between the resistant CAYMAN (A) and CUBA-DELTA (B) strains against the NEW ORLEANS susceptible strain. Transcript levels are shown by the white columns and gene copy numbers by the grey columns. Error bars represent 95% confidence intervals. In the CAYMAN strain, 8 transcripts for odorant binding proteins (OBPs) were up-regulated. Furthermore, when the CAYMAN population was analysed alone, the GO term ‘odorant binding’ was the most differentially represented term in the up-regulated set of transcripts (Figure 2). OBPs facilitate the passage of semio-chemicals across the antennae, and other sensory appendages, to the olfactory neurones. To date, no specific role for OBPs in insecticide resistance has been demonstrated but this is not the first time that OBPs have been identified as being over expressed in insecticide resistant populations. A study of bendiocarb resistance in Anopheles identified an OBP gene that was overexpressed in Ghanaian resistant populations (S Mitchell, unpublished data). Gene amplification associated with elevated levels of ABC transporter and P450 transcripts Quantitative PCR was used to compare gene copy number between the two resistant and the susceptible Ae. aegypti strains. Using the same cut off of >2 fold change in expression, gene amplification was observed for two genes, CYP9J26 and the ABC transporter, ABCB4 (AAEL006717). The copy number of CYP9J26, measured using two different primer pairs, was between 6.5 and 8.1-fold higher in the resistant CUBA and CAYMAN strains respectively compared with New Orleans (Figure 3). Similarly, the ABCB4 gene was amplified approximately 7-fold in the Cayman strain relative to the New Orleans (Figure 3). Discussion Pyrethroid resistance is widely distributed in Ae. aegypti throughout its range but relatively little is known about the mechanisms responsible for this resistance. Target site resistance is present in both the Cuban and Cayman populations although neither of the two mutations, 1016I or 1534C, were fixed in either population, despite several rounds of laboratory selection with deltamethrin in the Cuban strain. Interestingly the 1534C mutation has recently been shown to confer a selective advantage against type I pyrethroids but not affect the sensitivity to type II pyrethroids such as deltamethrin [25]. As far as we are aware, electrophysiological experiments have not been performed to examine the impact of the V1016I substitution. Thus, target site resistance may be partially responsible for the high levels of deltamethrin resistance in both these populations but it is likely that other mechanisms are involved. In this study, a microarray containing probes for the vast majority of annotated genes in the Ae. aegypti genome was used to compare gene expression in the two Caribbean populations with a standard lab susceptible strain. A potential limitation of this approach is the use of a single laboratory susceptible strain that originated from the United States. Ideally, a range of susceptible strains including wild populations from similar genetic regions would be included in the study. Unfortunately, such strains are becoming increasingly difficult to find. The differential gene expression observed may be partially attributed to the different genetic background of the strains although, encouragingly, earlier experiments have shown no significant difference in expression of detoxification genes between New Orleans and another well established laboratory susceptible strain, Rockefeller [19]. In light of this potential criticism, the analysis focused primarily on genes that were up-regulated in both resistant populations. Interestingly, a smaller number of genes were found differentially transcribed in the Cuban strain, which had been subject to extensive laboratory selection, than the Cayman strain, which was resistant upon colonization. By using GO term enrichment analysis, the functions or processes that were enriched in the subsets of genes up or down regulated in the pyrethroid resistant populations were identified. Only one GO term was significantly enriched in this analysis. Eleven transcripts with the GO term GO:0004497, monooxygenase activity, were found amongst the subset over expressed in both insecticide resistant populations. This supports the well documented role of cytochrome P450s in conferring pyrethroid resistance [26]. Aedes aegypti has an extensive repertoire of between 160 and 180 CYP genes [19]. The uncertainty over the exact gene count is partly due to the fragmented nature of the Ae. aegypti genome assembly; several supercontigs containing clusters of P450s most likely represent alternative haplotypes. P450s contained within these putative duplicate clusters have been assigned independent accession numbers in VectorBase but named as allelic variants of the same P450 by the P450 nomenclature committee (designated v1 or v2). This issue is discussed further in the supplementary material of Strode et al (2008) [19]. For ease of discussion, in the current manuscript, official P450 nomenclature has been used to discuss the P450 family and the v1/v2 nomenclature omitted. However, the finding that two variants of the same gene frequently show similar fold changes in expression, adds confidence to the analysis as the probes were not designed to be able to distinguish allelic variants. Seven cytochrome P450 genes were up-regulated in both resistant populations. This included two genes in the CYP6 clade, CYP6BB2 and CYP329B1, but neither was confirmed by qPCR. The five remaining P450 genes belonged to the CYP9J family. Four of these have been found to be up-regulated in resistant strains from other geographical localities. CYP9J9 and CYP9J10 were found over expressed in Thai and Latin American populations of Ae. aegypti. CYP9J28 is over expressed in pyrethroid resistant populations from Peru and Mexico [19], [27], and has also been shown to be over expressed in pyrethroid resistant Ae. aegypti from Vietnam (Warr and Ranson, unpublished data). CYP9J27 is over expressed in Thailand [19] and is also one of the candidates emerging from the Vietnam study. Other CYP9J genes have also been implicated in resistance (Figure 4). In fact, in total, ten CYP9J genes have been found over expressed in at least two pyrethroid resistant populations. Four of these, CYP9J24, 26, 28 and 32, have now been biochemically characterized [28] and have all been shown to metabolize pyrethroids (a single CYP6, CYP6CB1 was also expressed but had no activity against this insecticide class). 10.1371/journal.pntd.0001692.g004 Figure 4 Unrooted distance neighbour joining tree showing phylogeny of Aedes aegypti CYP9 genes. Nodes with >70% bootstrap support (500 pseudoreplicates) are indicated. Sequences in blue are up-regulated in both CUBA-DELTA and CAYMAN strains from the current study. Sequences underlined are over expressed in ≥2 populations versus laboratory susceptible strains from previously published studies (see text for further details). Sequences marked with * have proven ability to metabolise pyrethroids (Stevenson et al, submitted). Aside from the CYP9s, the only additional clade of Ae. aegypti P450s that is found repeatedly over expressed in resistant strains is the CYP6Z subfamily. Genes CYP6Z6, Z8 and Z9 have been found over expressed in multiple populations from southeast Asia, Latin America and Caribbean [19], [27], [29] and CYP6Z8 and Z9 were over expressed in Cuba and Cayman populations respectively in the current study. The CYP9 family in Ae. aegypti is greatly expanded compared to other insect species with over three times as many members as found in An gambiae and nearly six times as many as in D. melanogaster [19], [30]. The degree of genetic redundancy in the P450 family of Ae. aegypti makes it unlikely that a single gene responsible for pyrethroid resistance in all strains will be detected, particularly if resistance is emerging independently in different populations. However, the identification of a small subset of genes, consistently over expressed in resistant populations does suggest it should be possible to develop specific inhibitors of these metabolic pathways that could be used as insecticide synergists. Gene amplification was associated with the over expression of one of the P450 genes, CYP9J26, in both strains, with approximately 7 –fold increase in copy number compared to the susceptible strain. Increased gene copy numbers have been associated with P450 mediated resistance in An. funestus, D. melanogaster and Myzus persicae [31], [32], [33]. The increased transcript levels of the other CYP9 genes were not associated with an increase in gene copy number although copy number polymorphisms appear to be common in the P450 family in Ae. aegypti (Strode et al, 2008). CYP9J26, 27 and 28 are arranged sequentially in the Ae. aegypti genome within a large cluster of CYP9 genes on supercontig 1.1188. It is not yet known if the CYP926 duplications are found in tandem. Multiple copies of the ABC transporter were also present in the Cayman resistant population. Gene amplification is being increasingly recognized as an important mechanism conferring metabolic resistance to insecticides with examples reported in all the major families of detoxification enzymes from several insect species [34]. Resistance to pyrethroid insecticides is now widely established in Ae aegypti populations throughout its distribution [2]. Elucidating the mechanisms responsible for this resistance will facilitate resistance monitoring and pave the way for the development of effective resistance reversal approaches. The microarray approach used in this study is not itself a field applicable screening approach. However, using this tool to analyse additional pyrethroid resistant populations will help define a subset of genes that are responsible for pyrethroid resistance. As discussed above, although it is unlikely that a single diagnostic mutation will be detected, the identification of a panel of candidate resistance associated genes is an important prerequisite for developing simple, molecular diagnostics that are urgently needed by dengue control programmes. Given the key role that pyrethroids play in controlling this disease vector, and the lack of affordable, acceptable alternative insecticides, it is imperative that efforts are made to monitor for resistance and reduce the impact that this resistance may impose on vector control interventions. Supporting Information Table S1 Primer sequences used for quantitative PCR. (DOCX) Click here for additional data file. Table S2 List of genes differentially expressed in both the CAYMAN and CUBA-DELTA strain compared to New Orleans. (XLS) Click here for additional data file. Table S3 List of genes differentially expressed in the CAYMAN strain compared to New Orleans. (XLS) Click here for additional data file. Table S4 List of genes differentially expressed in the CUBA-DELTA strain compared to New Orleans. (XLS) Click here for additional data file. Table S5 Differential expression of genes determined in quantitative PCR experiments. Data analysis was performed according to the ΔΔCT method taking into account PCR efficiency [35] and normalized with two control genes 60S ribosomal protein L8 (AAEL000987) and the 40S ribosomal protein S7 (AAEL009496). Each experiment was analyzed separately. Results were expressed as mean transcription ratios ± SD (n = 3). (XLS) Click here for additional data file.
