Gene Amplification, ABC Transporters and Cytochrome P450s: Unraveling the Molecular Basis of Pyrethroid Resistance in the Dengue Vector, Aedes aegypti

Background Pyrethroid insecticides are widely utilized in dengue control. However, the major vector, Aedes aegypti, is becoming increasingly resistant to these insecticides and this is impacting on the efficacy of control measures. The near complete transcriptome of two pyrethroid resistant populations from the Caribbean was examined to explore the molecular basis of this resistance. Principal Findings Two previously described target site mutations, 1016I and 1534C were detected in pyrethroid resistant populations from Grand Cayman and Cuba. In addition between two and five per cent of the Ae. aegypti transcriptome was differentially expressed in the resistant populations compared to a laboratory susceptible population. Approximately 20 per cent of the genes over-expressed in resistant mosquitoes were up-regulated in both Caribbean populations (107 genes). Genes with putative monooxygenase activity were significantly over represented in the up-regulated subset, including five CYP9 P450 genes. Quantitative PCR was used to confirm the higher transcript levels of multiple cytochrome P450 genes from the CYP9J family and an ATP binding cassette transporter. Over expression of two genes, CYP9J26 and ABCB4, is due, at least in part, to gene amplification. Significance These results, and those from other studies, strongly suggest that increases in the amount of the CYP9J cytochrome P450s are an important mechanism of pyrethroid resistance in Ae. aegypti. The genetic redundancy resulting from the expansion of this gene family makes it unlikely that a single gene or mutation responsible for pyrethroid resistance will be identified in this mosquito species. However, the results from this study do pave the way for the development of new pyrethroid synergists and improved resistance diagnostics. The role of copy number polymorphisms in detoxification and transporter genes in providing protection against insecticide exposure requires further investigation.


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.

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 30fold 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 (26uC, 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 TM 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 280uC until microarray hybridizations.

Hybridizations, data acquisition and statistical analysis
Hybridizations were made to the 'Liverpool Aedes aegypti Agilent 8615K v1' microarray (A-MEXP-1966) designed by the Liverpool School of Tropical Medicine. Each array contains 60mer oligoprobes 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 65uC 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,0.01 after multiple testing correction were considered as significant. Descriptions and GO-terms of transcript-IDs were extracted from VectorBase [15] using BIOMART [16] and completed with Blast2GO software (BioBam Bioinformatics S.L. (Valencia, Spain)) [17]. GO term Enrichment analysis was performed on the significantly upregulated genes (72% of transcripts present on microarray have GO-terms) using Blast2GO software with Fisher's exact test and false discovery rate (FDR),0.05.

Quantitative Real-time PCR
Selected microarray data were validated using quantitative reverse transcription PCR (qRT-PCR). Primers were designed using the Oligo7 Primer Analysis Software (Molecular Biology Insights,Cascade, CO, USA) based on cDNA sequences retrieved from VectorBase. An aliquot of 4 mg total RNA from each of the three biological replicates, for each strain, served as a template for cDNA synthesis with Superscript III (Invitrogen, Carlsbad, CA, USA) using oligo-dT 20 , according to the manufacturer's instructions. The resulting cDNAs were diluted 20 times in ultra-high quality water for qRT-PCR reactions using a MiniOpticon System (Biorad, Hercules, CA, USA). PCR reactions of 25 ml contained Fast Start SYBR Green Master Mix (Roche, Penzberg, Germany), 0.3 mM of each primer (Table 1) and 5 ml of diluted cDNA. Melt curve analysis was performed to test the specificity of amplicons. A serial dilution of cDNA was used to generate standard curves for each gene in order to assess PCR efficiency and quantitative

Author Summary
Dengue is the most rapidly spreading arboviral infection of humans and each year there are 50-100 million cases of dengue fever. There is no vaccine or drug to prevent dengue infection so control of the mosquitoes that transmit this virus is the only option to reduce transmission. Removing mosquito habitats close to human homes can be effective but in reality most dengue control programmes rely on a small number of chemical insecticides. Therefore, when the mosquito vectors develop resistance to the available insecticides, dengue control is jeopardized. In this study we examined the causes of resistance to the insecticide class most commonly used in mosquito control, the pyrethroids. We found that a group of genes, which have been implicated in detoxifying these insecticides in other populations of dengue vectors, were highly over expressed in both these Caribbean populations and we investigated the molecular basis of this increased expression. The next steps, which will be a considerable challenge, are to utilize this information to develop effective means of restoring insecticide susceptibility in dengue vectors.  Table  S1. The fold-change of each target gene, normalized to the 60S ribosomal protein L8 (AAEL000987) and 40S ribosomal protein S7 (AAEL009496), and relative to NO, was calculated according to the 2 2DDCT method incorporating PCR efficiency [18]. In most cases, two independent primer sets were used for each gene (Table S1).

Gene copy number analysis
Genomic DNA (gDNA) from three batches of ten adult mosquitoes from each strain was extracted using DNAzol (Invitrogen) according to the manufacturer's instructions. DNA quality and quantity was assessed by Nanodrop ND-1000 spectrophotometry and by running an aliquot on a 1.5% agarose gel. Quantitative PCR reactions were performed as described above on the same genes chosen for transcript analysis. The PCR efficiency for target and control genes (RPS7 and RPL8) was calculated from standard curves generated from a pool of gDNA for all three strains. 150-200 ng of gDNA was used as template and primer concentrations were between 100-300 nM for all genes. Primer sequences are provided in Table S1. The relative copy number fold-change was calculated using the 2 2DDCt method [18].

kdr genotyping
Two target site mutations, both previously associated with pyrethroid resistance in Ae. aegypti, are present at high frequencies in the two resistant strains from the Caribbean. In the CAYMAN strain, kdr frequencies of 0.79 for the 1016I allele, and 0.68 for the 1534C allele have been reported previously [9]. In the current study, 38 individuals from the CUBA-DELTA strain were genotyped and the frequency of the resistant alleles were 0.51 (1016I) and 0.88 (1534C).    (Table S2). 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.

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,0.05) only the GO term designating moooxygenase activity was significantly differentially represented.
This initial screening via GO terms was followed up with a manual examination of the putative functions of each of the 107 transcripts up-regulated and 99-down regulated in each population. This gene set, ranked by fold change in the CAYMAN population, is listed in Table S2. Genes with putative detoxification functions are listed in Table 1. Detoxification genes comprised 15.8% of the commonly up-regulated subset but were not represented at all in the down -regulated subset of the commonly expressed genes. Several additional detoxification genes were found up-regulated in either the Cayman or Cuba strains (Table 1).

Detoxification genes
A total of 18 and 13 CYPs are over expressed .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 Pyrethroid Resistance in Aedes aegypti www.plosntds.org 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).
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 upregulated 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  [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).
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.    [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 6 SD (n = 3). (XLS)