Gene enrichment in individual body part versus whole organism

Transcripts showing enriched transcription in each of the four body parts were determined using a multiple test correction significance cut-off of p ≤0.05 for both the resistant and susceptible strains. As expected, a clear positive correlation can be observed for local transcription between the resistant and susceptible strains in each body part. This is illustrated in Figure 1, where the vast majority of probes follow a y = x trend.

In all body parts, similar numbers of transcripts showed enriched or depleted transcription compared to the whole organism, with a range of 1.4% to 4.4% of the total probes (Additional file 2: Table S1). The magnitude of change in transcription of individual transcripts in the midgut and remaining body parts is relatively low with no transcripts exceeding log₂ fold change of 2, compared to the abdomen integument and malpighian tubules where 7 and 42 genes respectively are above log₂ 3-fold enriched (Figure 1). The full lists of transcripts enriched in each body parts are listed in Additional file 3: Table S2.

The An. gambiae genome contains 211 genes encoding four major detoxification families (111 cytochrome P450s, 31 GSTs, 43 CCES (only putatively catalytically active enzymes included) and 26 UGTs), together comprising 1.5% of the probes on the Agilent array [6]. Each of the dissected body parts have an over representation of detoxification transcripts with more members of these gene families overtranscribed relative to the whole organism than found in the ‘remaining’ (undissected body part). Additional file 2: Table S1 shows that 8.5% of the total number of detoxification transcripts are enriched in the malpighian tubules of the resistant strain, 13.7% in the midgut and 7.6% in the abdomen integument. In contrast, the ‘remaining’ (undissected body part) show a depletion of detoxification transcripts in the overtranscribed subset when compared to the whole organism. These data reinforce the importance of the selected body parts in xenobiotic detoxification.

The cytochrome P450 family has been most strongly linked with insecticide resistance in Anopheles mosquitoes, with several enzymes capable of detoxifying insecticides from more than one class [2, 8]. Identifying the primary sites of transcription of this enzyme family will aid prediction of function [18] and help identify the key organs largely responsible for insecticide detoxification in resistant mosquitoes. Within the P450 gene family, body part enrichment shows some relationship with the gene tree clustering (Additional file 4: Figure S2) with the CYP9J and 6P families largely enriched in the midgut, the CYP4Gs, 6Ys and 325Cs enriched in the abdomen integument and the CYP6Z family enriched in the malpighian tubules. The diversity in enrichment patterns within this gene family led us to look specifically at thirteen cytochrome P450s that had been implicated in insecticide resistance to see if these were enriched in a particular body part (Figure 2). The criteria for inclusion of these P450s as candidate insecticide resistant genes was that they had been found to be significantly over expressed in pyrethroid and/or DDT resistant An. gambiae populations in more than one independent study. Four of these transcripts were significantly (p < 0.05) and highly (Log₂ Fold Change >1.5) enriched in the malpighian tubules (CYP6M3, CYP6Z1, CYP6Z2 and CYP6Z3) compared to whole organism and one was enriched in the midgut and two in the abdomen integument (CYP4H24 and CYP4G16/CYP4G17 respectively). The remaining 5 showed no significant tissue enrichment.

The enrichment of a subset of seven cytochrome P450s in particular body parts was confirmed by qPCR, although for two of these midgut enriched transcripts, a much greater over expression was observed by qPCR than for microarray (Additional file 5: Figure S3) (A level of discordance between qPCR and microarray is frequently observed [19] (but it in the current study, the direction of change was consistent between the two methodologies for all seven genes). Antibodies were available for three cytochrome P450s, which were used to confirm the major sites of transcription within the abdomen integument of the Tiassalé resistant strain. In agreement with microarray transcription data, CYP6Z1 and CYP6Z2 were detected in the malpighian tubules of resistant mosquitoes, whilst CYP4G17, identified as enriched in the abdomen integument by both microarray and qPCR, was found only in the oenocytes, a cellular layer located under the abdomen cuticle (Figure 3). A CYP4G from Drosophila melanogaster, CYP4G1, has also been shown to be highly enriched in oenocyte cells; this Drosophila enzyme catalyses a key step in the formation of cuticular hydrocarbons [20]. The potential role of over expression of 4G17, and its paralogue 4G16, in altering the cuticular structure in insecticide resistant mosquitoes is currently under investigation.

A new list of candidate insecticide resistance genes from the body part specific arrays

The two strains used in this study originate from sites separated by approximately 2,500 km, so substantial variation in gene transcription between them is to be expected, regardless of the difference in their insecticide resistance profile. Thus caution must be applied when correlating gene transcription levels with the resistance phenotype. However, with a goal of identifying further candidate insecticide resistance genes from the direct comparisons between dissected body parts from the resistant and susceptible strain for further functional validation, we applied a stringent selection process to derive a new gene list. This involved selecting only genes detected as enriched in body parts from Tiassalé versus N’Gousso via both the limma and GaGa methods with an adjusted p value ≤0.001 and setting a cut off of 1.4-fold differential expression (see Methods section). This identified a list of 134 transcripts, representing 105 genes transcribed at higher levels in the Tiassalé strain and 16 genes with higher transcription in the susceptible N’Gousso strain (Additional file 6: Table S3). Eleven of these transcripts were selected for qPCR validation yielding a positive correlation with the array (Pearson’s correlation (r = 0.870)) (Figure 4).

Mosquito rearing conditions

The An. gambiae used in these experiments were all reared under standard insectary conditions at 27°C and 70-80% humidity under a 16:8 hour photoperiod. The N’Gousso strain is originally from Cameroon and is susceptible to all classes of insecticide [17]. N’Guosso is the M molecular form of A gambiae, recently re-classified as a separate species, Ano coluzzi[24]. In contrast, the Tiassalé strain from Côte D’Ivoire is resistant to all classes of insecticide [2, 3]. This strain was colonised from the field site in 2012 and is a mixture of the M and S molecular forms. At the time of the study, the LD₅₀ for the Tiassalé strain was 68 and 81 fold higher than the corresponding value for the N’Gousso strain for permethrin and deltamethrin respectively. Further details of the resistance profile of this strain are contained within references by Edi et al. [2, 3].

Microarray experiments

RNA was extracted from three dissected body parts: the malpighian tubules, the midgut and the abdomen integument (containing the fat body but also epidermal, neuronal, muscle and oenocyte cells) with the remaining undissected body parts forming a fourth sample group. Mosquitoes were collected between the hours of 8AM and 2 pm and dissected immediately on a CO₂ block. Post dissection, each body part was added to extraction buffer from the PicoPure RNA extraction kit, heated for 30 minutes at 42°C and frozen at -80°C as per manufacturers instructions. Each biological replicate for each strain consisted of RNA, extracted using PicoPure RNA Isolation kit (Arcturus), from 12 3–5 day old non-blood fed, presumed mated females. The quantity and quality of the RNA was assessed using a nanodrop spectrophotometer (Nanodrop Technologies UK) and Bioanalyser (Agilent) respectively. Four biological replicates were prepared for each body parts per strain. RNA from the four dissections was pooled according to the proportion of RNA extracted from each body part to reconstitute the ‘whole organism’ (7%, 6%, 24% and 63% RNA from abdomen integument, malpighian tubules, midgut, and remaining material respectively). The use of a reconstituted reference sample minimised potential sources of bias that could have arisen from circadian changes in gene expression and changes in the proportion of the M or S molecular form in the different biological replicates. 100 ng of RNA was amplified and labelled with Cy3 and Cy5, using the Two colour low input Quick Amp labelling kit (Agilent) following the manufacturers instructions. Samples were then purified (Qiagen) with the cRNA yield and quality assessed using the nano-drop and Bioanalyser respectively. RNA from each Tiassalé body part was competitively hybridised with the respective N’Gousso body part, as well as each body part from the resistant and susceptible strain being compared to the re-constituted whole organism (Additional file 1: Figure S1). Dye swaps were performed on two out of four technical replicates for each array, to correct for dye bias.

Labelled cRNAs were hybridised to the whole genome 8×15k An gambiae array (ArrayExpress accession number A-MEXP-2211). Microarray hybridisation, washing and scanning were performed according to previously described protocols [8].

Microarray analysis

The resulting data were analysed using R. Within-array normalisation was carried out by Loess, and between array normalisation by Aquantile. Signals were corrected for dye automatically. The limma package [25] was used to fit linear models to the normalised data. In the case of complete loop designs a design matrix was used to infer the appropriate contrast matrices for each array. All parameters used were default. A bespoke pipeline using the GaGa package [26] was used to fit gamma-gamma models of variation to normalised corrected signals, in order to assign probes to one of two patterns of expression X equals Y or X does not equal Y, where X represents the resistant population arrays and Y the susceptible arrays. These data were subsequently used to assess enrichment in each expression pattern, through GO term analysis using the TopGO package [27]. A standard FDR adjusted p value cut off of p ≤0.05 was applied to all data describing localisation of detoxification candidates. A second stringent selection method was used to reduce the probe list based on previously published methodology [19], requiring that the following criteria were met: adjusted p-value ≤0.001, raw fluorescence intensity > median, and Tiassalé vs. N’Gousso ±0.485 Log₂ fold change between the strains. All candidates selected also demonstrated a positive GaGa analysis fold change, indicative of higher transcript localisation in the resistant tissue to the susceptible, thereby utilising all available array data.

RT-qPCR

RNA (4 μg) from each biological replicate was reverse transcribed using Oligo dT (Invitrogen) and Superscript III (Invitrogen) according to manufacturers instructions. Quantitative real-time PCR was performed using SYBR Green Supermix III (Applied Biosystems) using an MX3005 and the associated MxPro software (Agilent). Primer Blast (NCBI) [28] was used to design primer pairs (Additional file 7: Table S4). Where possible, primers were designed to span an exon junction but this was not possible for six of the P450 genes (CYP325A1, CYP6P3, CYP4G17, CYP6Z3, CYP12F2 and CYP6Z2) due to the high degree of polymorphisms in their DNA sequence. Each 20 μl reaction contained 10 μl SYBR Green Supermix, 0.3 μM of each primer and 1 μl of 1:10 diluted cDNA. Standard curves were produced using whole N’Gousso cDNA, in 1, 1:5, 1:25, 1:125 dilutions, (48.3 ng/μl to 0.386 ng/μl). qPCR was performed with the following conditions: 3 minutes at 95°C, with 40 cycles of 10 seconds at 95°C and 10 seconds at 60°C. All amplification efficiencies of designed primers were within acceptable range (90-120%), following MIQE guidelines [29].

Preparation of antibodies

Fragments encoding unique peptides for CYP6Z1 and CYP4G17 were cloned into the pET 16b vector. Upon expression, the resulting His-tagged peptide was purified to homogeneity by Ni-NTA affinity chromatography and used to raise rabbit polyclonal antibodies. The CYP6Z1 peptide sequence was: VALRDLNNPDSFINNIRTAGVFLCPGLLKFTGINSLSPPMKKFTTEVISSHLHQRETGQVTRKDFIQMLTDLRRKAGSSGEETLTDA and the CYP4G17 peptide: KRQLKIHLRLDPLFNLTGVKKEQERLLQIIHGLTRKVVREKKQLYERQMAEGKMPSPSLTEIIGKEEKPGEGQLGGSPAFISQ. The antibody for CYP6Z2 was a gift from Dr Mark Paine (LSTM, UK). Rabbit antibodies to CYP6Z2 were prepared to the C terminal peptide sequence MRIDHRK by Moravian Biotechnology, Brno, Czech Republic.

Immunofluorescence and microscopy

Female mosquitoes (3–5 days old) were fixed in cold solution of 4% PFA (methanol free, Thermo scientific) in phosphate-buffered saline (PBS) for 4 h and then were cryo-protected in 30% sucrose/PBS at 4°C for 12 h. Finally, mosquitoes were immobilized in O.C.T. (Tissue-Tek, SAKURA) and stored at -80°C until use.

Immunofluorescent analysis, followed by confocal microscopy, was performed to longitudinal sections of frozen pre-fixed mosquito specimens. More detailed, 10 mm sections, obtained in Leitz kryostat 1720 digital, were washed (3 × 5 min) with 0,05% Tween in PBS and blocked for 3 h in blocking solution (1% Fetal Bovine Serum, biosera, in 0,05% Triton/PBS). Then, the sections were stained with rabbit primary antibodies in 1/500 dilution, followed by goat anti-rabbit (Alexa Fluor 488, Molecular Probes) (1/1000) that gave the green colour. Also To-PRO 3-Iodide (Molecular Probes), which stains DNA specifically (red colour), was used, after RNAse A treatment. As controls, pre-immune serums (in 1/500 dilutions) and anti-rabbit (Alexa Fluor 488, 1/1000) were tested, in parallel with a-P450’s to check specificity of each primary antibody. Finally images were obtained on Leica TCS-NT Laser Scanning microscope using the 40-objective.

Supporting information

The data sets supporting the results of this article are included within the article and its additional files. All microarray datasets are MIAME compliant [30] and deposited in ArrayExpress and VectorBase (Accession numbers E-MTAB-2808).