Skip to main content
Download PDF

RNA-seq analyses of changes in the Anopheles gambiae transcriptome associated with resistance to pyrethroids in Kenya: identification of candidate-resistance genes and candidate-resistance SNPs

Parasites & Vectors volume 8, Article number: 474 (2015) Cite this article

Abstract

Background

The extensive use of pyrethroids for control of malaria vectors, driven by their cost, efficacy and safety, has led to widespread resistance. To favor their sustainable use, the World Health Organization (WHO) formulated an insecticide resistance management plan, which includes the identification of the mechanisms of resistance and resistance surveillance. Recognized physiological mechanisms of resistance include target site mutations in the para voltage-gated sodium channel, metabolic detoxification and penetration resistance. Such understanding of resistance mechanisms has allowed the development of resistance monitoring tools, including genotyping of the kdr mutation L1014F/S in the para gene.

Methods

The sequence-based technique RNA-seq was applied to study changes in the transcriptome of deltamethrin-resistant and -susceptible Anopheles gambiae mosquitoes from the Western Province of Kenya. The resulting gene expression profiles were compared to data in the most recent literature to derive a list of candidate resistance genes. RNA-seq data were analyzed also to identify sequence polymorphisms linked to resistance.

Results

A total of five candidate-resistance genes (AGAP04177, AGAP004572, AGAP008840, AGAP007530 and AGAP013036) were identified with altered expression between resistant and susceptible mosquitoes from West and East Africa. A change from G to C at position 36043997 of chromosome 3R resulting in A101G of the sulfotransferase gene AGAP009551 was significantly associated with the resistance phenotype (odds ratio: 5.10). The kdr L1014S mutation was detected at similar frequencies in both phenotypically resistant and susceptible mosquitoes, suggesting it is no longer fully predictive of the resistant phenotype.

Conclusions

Overall, these results support the conclusion that resistance to pyrethroids is a complex and evolving phenotype, dependent on multiple gene functions including, but not limited to, metabolic detoxification. Functional convergence among metabolic detoxification genes may exist, with the role of each gene being modulated by the life history and selection pressure on mosquito populations. As a consequence, biochemical assays that quantify overall enzyme activity may be a more suitable method for predicting metabolic resistance than gene-based assays.

Background

Over the past decade, several initiatives including the Global Fund, the President’s Malaria Initiative, private foundations and national governments supported a massive scale-up of antimalarial interventions in Africa [1, 2]. These control programs targeted malaria vectors, through insecticide-treated nets (ITN) and indoor residual spraying (IRS), as well as human hosts by improving diagnosis and implementing artemisinin-combination treatments (ACT). As a result, the annual number of malaria-related deaths in sub-Saharan Africa decreased by 49 % between 2000 and 2012 [3]. However, malaria still kills more than half a million people a year [3] and weakening of malaria control programs could favor malaria resurgence [4]. As a consequence, the identification of elements that could threaten the sustainability of malaria control strategies is critical to continue the fight against this disease. Currently, the core strategies of vector control (IRS and ITN) rely on insecticides [5]. The World Health Organization (WHO) recommends the use of four classes of insecticides in IRS (pyrethroids, organochlorine, organophosphate and carbamate), while only pyrethroids are approved for use on ITNs [6]. Extensive use of insecticides imposes selection pressure on mosquito populations for increased resistance. Resistance to insecticides is widespread geographically in Africa and involves primarily, but is not limited to, pyrethroids [6–8]. As a response to this situation, the WHO formulated an action plan to support the sustainability of control programs focused on the use of insecticides [6]. The understanding of the mechanisms of insecticide resistance and the monitoring of resistance are two of the five pillars of this plan [6]. Historically, the identification of the mechanisms of resistance has been important for developing molecular monitoring tools of resistance [9].

One of the main malaria vectors in sub-Saharan Africa is Anopheles gambiae sensu strictu (hereafter referred as An. gambiae). An. gambiae exists as two molecular forms (M and S), emerging as distinct species, mainly due to reduced fitness of hybrids in nature [10, 11]. The S form, named An. gambiae, has the broadest distribution, occurring in West and East Sub-Saharan Africa; the M form, named Anopheles coluzzii, is primarily confined to West Africa, with the exception of northern Zimbabwe [10]. The same types of resistance mechanisms have been identified in An. gambiae and An. coluzzii, including behavioral and physiological resistance. Physiological resistance includes target site mutation, metabolic detoxification, and penetration resistance. The target site for pyrethroids is the para voltage-gated sodium channel [12]. Two alternative amino-acids changes at position 1014 of the para gene (L1014S and L1014F) have been associated with resistance [12, 13]. Originally detected in West (L1014F) and East (1014S) Africa, both mutations are now found across mosquito populations from sub-Saharan Africa, but at different frequencies. The L1014F mutation is rare in East Africa, but approaches fixation in many countries of West Africa; the L1014S allele has increased in frequency across all sub-Saharan Africa in the past decade [8, 9]. An additional mutation (N1575Y) was detected recently in An. coluzzii mosquitoes from Burkina Faso that harbored the L1014F mutation [14]. The N1575Y mutation is thought to either enhance resistance due to the L1014F mutation or compensate for its fitness costs [14]. The frequency of this mutation is unknown in An. gambiae.

Metabolic detoxification occurs when pyrethroids are catabolized or sequestered and eliminated before they reach the voltage-gated sodium ion channel [12]. Enzymes encoded by three large gene families (esterases, P450 mono-oxygenases [CYPs] and glutatione-S transferases [GSTs]) play important roles in insecticide detoxification [11, 15]. Gene-expression studies identified several detoxification genes over-expressed in resistant vs susceptible An. gambiae and An. coluzzii laboratory strains and wild populations, with the products of two genes (CYP6P3 [AGAP002865] and CYP6M2 [AGAP008212]) shown to be able to metabolize pyrethroids [16, 17]. A state of stress, manifested primarily through lipid peroxidation, is also seen upon insecticide exposure [18, 19]. Target site mutation and metabolic detoxification may co-occur with varying frequencies in natural populations and their relative contribution to resistance phenotype may be influenced by the intensity of vector control strategies and the breeding environments [20, 21]. Furthermore, other physiological mechanisms of resistance have been proposed, including thickening of the cuticle and modification of the digestive tract lining, which may help reduce insecticide penetration and absorption (penetration resistance) [8, 22]. Behavioral modifications, such as earlier biting and outdoor feeding, which results in mosquitoes avoiding insecticides, have also been recognized as important factors contributing to increased resistance (behavioral resistance) [23]. However, current data are insufficient to determine whether behavioral resistance traits are genetic or adaptive and the difficulty in accurately measuring mosquito behavior in the wild has limited the understanding of behavioral resistance [23].

We examined the gene expression profile of deltamethrin-resistant and -susceptible mosquitoes from the Western Province of Kenya by RNA-seq to further the understanding of resistance mechanisms and possibly characterize markers for resistance surveillance. The application of this sequence-based approach allowed us to interrogate absolute changes in transcript accumulation profiles and investigate single nucleotide polymorphisms (SNPs) that may be associated with insecticide resistance. A step-wise filtering approach, including comparison to literature data [24], was used to identify a total of five candidate resistance genes. Additionally, for SNPs, the kdr L1014S mutation was detected at similar frequencies in phenotypic pyrethroid-resistant and susceptible-mosquitoes and the non-synonymous change A101G in sulfotransferase gene AGAP009551 was found to be associated with the resistance phenotype.

Methods

Mosquito samples

Anopheles gambiae larvae were collected in the Western Province of Kenya, in breeding sites around the localities of Bungoma, Busia and Emutete in 2012 (Fig. 1). Multiple mosquito samplings in different localities within the Western Province allows capturing the mosquito genetic variability of this region and makes gene expression analysis results more broadly applicable. A previous survey identified wide-spread resistance to deltamethrin in mosquitoes from this area, with mortality ranging between 66.5 and 78.7 % [25]. Larvae were transported to the insectary of the Centre for Global Health Research, Kenya Medical Research Institute (KEMRI) in Kisumu, and reared to adulthood for the resistance bioassay. The prevalent mosquito species in this area is An. gambiae [24]. Mosquitoes of the Kisumu strain also were used. The Kisumu strain originated in the early 1990s and is the pyrethroid susceptible standard strain [26].

Fig. 1
figure 1

Map of the Western Province of Kenya. Localities in the Western Province of Kenya around which multiple larvae collections occurred are shown with a red circle. Kagamega and Kisumu are the capitals of the Western and Nyanza Provinces, respectively. Nearby countries are shown with a square and in purple. Main roads are in yellow

Full size image

Phenotypic resistance to deltamethrin

Resistance to deltamethrin was tested by the standard WHO tube test [7] on single mosquitoes reared from field-collected larvae. Mosquitoes alive 24 h after the 60-min insecticide-exposure time were classified as resistant. Susceptible mosquitoes were those, which were knocked-down early after the insecticide exposure. Knocked-down mosquitoes were observed for 2–3 min, and the tubes were moved gently to confirm mosquitoes showed no signs of recovery; morphological signs of distress were also checked (i.e. loss of legs) before considering mosquitoes as susceptible. Upon collection, mosquitoes were stored in RNA-later (Ambion) to avoid RNA degradation. This strategy for collecting resistant- and susceptible-mosquitoes allows sampling mosquitoes from the same locality and minimizes the impact of heterogeneous genetic background and environmental conditions on gene accumulation profiles. Additionally, early knock-down has already been shown to be a valid approximation for susceptibility [25, 27, 28]. However, this phenotyping method will not allow differentiation between insecticide-induced and constitutive differential expression between resistant- and susceptible- mosquitoes. Constitutive differential expression will be investigated by comparing data from field-collected mosquitoes and mosquitoes of the susceptible Kisumu strain, which is expected to be highly inbred [26]. The WHO tube tests were conducted at the same hours of the day on all experiments to avoid influence of the circadian clock on transcript accumulation profiles [29].

Nucleic acids extractions

DNA was extracted from mosquito legs using the Fast Tissue-to-PCR kit (Fermentas). Total RNA was isolated from single mosquitoes using Trizol (Invitrogen) [30]. RNA concentration was measured by NanoDrop. RNA quality was analyzed on an Agilent 2100 Bioanalyzer.

Species discrimination between An. gambiae and An. arabiensis and genotyping of the para-sodium channel gene

Each deltamethrin-phenotyped mosquito was identified as An. gambiae or An. arabiensis by amplifying the species-specific rDNA [31]. Codons 1014 and 1575 of para sodium channel gene (AGAP004707) were analyzed on a total of 324 An. gambiae mosquitoes by direct sequencing of the fragments obtained by the polymerase chain reaction (PCR) with primers Adg1 and Adg2 [13] and 1575-F (5′ TAAACAGCCTATACGGGAAACG) and 1575-R (5′ CGAGGAATTGCCTTTAGAGGTTTCT), respectively.

RNA-seq library preparation and sequencing

A total of nine RNA-seq libraries were prepared from pooled RNA of 12 mosquitoes each (Table 1). Illumina paired-end RNA-seq libraries were prepared following a standard procedure and sequenced for 100 bp from each end on an Illumina HiSeq2500 at the DNA Technologies and Expression Analysis Core at the Genome Center of the University of California Davis [32]. Raw sequencing data have been deposited at NCBI’s Sequence Read Archine (SRA) under study accession number SRP052073.

Table 1 RNA-seq sample summary
Full size table

Data analyses

Raw read quality and contamination was evaluated using the Bioconductor package qrqc [33]. Scythe v.0.990 and Sickle v.1.200 [34] were used for Illumina adapter and quality-based trimming. Differential expression analyses followed the Tuxedo pipeline [35], which was run in Blacktie [36], using the VectorBase An. gambiae assembly P3 (and associated annotation AgamP3.7). Relationships among conditions and replicates were explored with Multiple Dimensional Scaling (MDS). MDS was computed using CummeRbund, a program within the Tuxedo pipeline [37]. MDS is a linear transformation method of variance stabilized count where the directions that maximize the separation (or discrimination) between the different samples are visualized. Transcript functional annotation was conducted using the biomart function in VectorBase [38] and AnoXcell [39]. Functional enrichment analyses was done using g: Profiler [40].

SNPs were called with Freebayes v. 0.9.4 and SnpEff [41, 42], following a previously described pipeline [43]. Prior to variant finding, multiple mapping reads (those with a mapping quality of zero) were removed from the Tophat output. Programs for estimating allele frequency from DNA-seq data of large pooled samples (>50) have been developed [43] and used successfully on RNA-seq data from highly inbred and not biologically-challenged samples [44, 45]. In our case, we expect RNA-seq data to include population variation and reflect differential expression of transcripts because we are using mosquitoes from the field and our experimental design compares mosquitoes of two different phenotypes (resistant and susceptible). Differential coverage across genes will result in unbalanced pool composition and biased allele frequency estimates when treating RNA-seq data as DNA-seq [43], as a consequence, we abstain from estimating SNP allele frequency using RNA-seq data through a DNA-seq-focused pipeline.

qRT-PCR validation of RNA-seq data

The accumulation levels of 18 transcripts were analyzed by quantitative real-time reverse transcription PCR (qRT-PCR). Specifically, cDNA was prepared using SuperScript III (Invitrogen) and random primers from pooled RNA of 3 resistant or 3 susceptible mosquitoes, including a total of 27 resistant and 27 susceptible mosquitoes. qRT-PCR reactions were run and analyzed as described previously using the S7 ribosomal protein gene as internal control [25, 46]. RNA from different phenotyped mosquitoes was used for RNA-seq library preparation and qRT-PCR experiments, providing for biological replicates. The Pearson correlation coefficient was calculated between fold changes in transcript accumulation levels between resistant and susceptible mosquitoes as obtained by qRT-PCR and RNA-seq, respectively [47].

Genotyping of candidate resistance genes and association with the resistance phenotype

The coding sequences of 39 genes harboring SNPs identified from RNA-seq libraries were analyzed in individual resistant and susceptible mosquitoes (Additional file 1). Specifically, genomic DNA was extracted from 54 resistant and 54 susceptible mosquitoes using the Fast Tissue-to-PCR kit (Fermentas). Genomic DNA was used as template in a PCR reaction with 11.5 μl of Master mix (Fermentas) and 10 μM of each forward and reverse primers (Additional file 1). PCR reactions were run in a MyCycler (Biorad) under the following conditions: 94 °C for 3 min followed by 40 cycles of 94 °C for 30 s, 50–54 °C for 45 s and 72° for 45–90 s. PCR products were sequenced directly, using the standard Sanger method [48]. The odds ratio test was applied to determine if the odds of each analyzed SNP differed significantly between resistant and susceptible mosquitoes [49].

Results

para gene polymorphism and “deltamethrin -resistant or -susceptible” phenotype

A total of 324 An. gambiae adults were genotyped at the para gene after their phenotype was established as “resistant” or “susceptible” to deltamethrin using the standard WHO tube test [7]. All resistant and >90 % of the susceptible mosquitoes were homozygous for the codon encoding for serine at position 1014 of the para gene. One susceptible and one resistant mosquito were homozygous for the codon encoding for phenylalanine at position 1014. The N1575Y mutation was not detected in any tested mosquito.

RNA-seq libraries of deltamethrin-resistant and susceptible mosquitoes: quality-control

Paired-end Illumina RNA-seq libraries were constructed and 100 bp were sequenced from each end from four RNA pools of 12 resistant or susceptible mosquitoes from the Western Province of Kenya and one pool of 12 mosquitoes of the susceptible Kisumu strain. The total number of reads per library aligned to the An. gambiae genome (AgamP3.7 gene set) ranged from 124,918,668 to 154,463,158 and showed no significant differences across samples (Table 1, Additional file 2). Variation in gene expression among libraries was visualized by MDS, which resulted in suggested two distinct clusters for the resistant and susceptible replicate samples (Fig. 2a). The sample consisting of pooled mosquitoes of the Kisumu strain was markedly distant from either of the two other groups (Fig. 2a). The greater distance between field mosquitoes and mosquitoes of the Kisumu strain may have resulted from different life history (e.g., exposure to insecticides and fluctuating temperature) or variation in larval rearing conditions. Distance among libraries within each cluster reflects variations among samples. Variations in gene expression profiles, which are suggestive of a broad spectrum of transcription within “resistant” and “susceptible” phenotypes, was expected because the WHO tube assay used to phenotype mosquitoes is a binary test that discriminates based on a 60 min exposure to a standard insecticide-discrimination dose (0.05 % deltamethrin). Finally, we compared gene expression levels estimated from RNA-seq and qRT-PCR for 14 genes. Eight genes (AGAP012984, AGAP002724, AGAP003714, AGAP004779, AGAP009472, CPLCG3 [AGAP008446], CYP6M2 [AGAP008212] and CYP6P3 [AGAP002865] were found to be significantly differentially expressed between resistant and susceptible mosquitoes by qRT-PCR and a significant positive correlation (Pearson correlation = 0.8025, p < 0.01) was found between expression values as detected by qRT-PCR and RNA-seq (Fig. 1b).

Fig. 2
figure 2

Quality control of RNA-seq data. a multidimensional scaling (MDS) plot showing variation among RNA-seq libraries. RNA-seq libraries from resistant, susceptible and Kisumu mosquitoes. b Results of qRT-PCR. The level of expression of 14 genes was measured by qRT-PCR from different resistant and susceptible mosquitoes than those used for RNA-seq. An axterix indicate genes significantly differentially expressed between resistant (green) and susceptible (pink) mosquitoes. c Pearson correlation between fold-changes in gene expression between resistant and susceptible mosquitoes as determined by qRT-PCR (X-axe) and RNA-seq (Y-axe)

Full size image

Differential expression of genes between field-caught resistant and susceptible mosquitoes

Candidate resistance genes were identified based on the assumption that they are significantly differentially expressed (DE) between resistant- and susceptible-mosquitoes (Additional file 3). A total of 2457 genes (18.36 % of 13361 genes tested) were found to be significantly DE (Additional file 4). Gene functional annotation revealed an enrichment in transcription-related functions, RNA processing, regulation of metabolic processes, chromatin organization, protein digestion and cellular response to stress among the 1373 genes more highly expressed in resistant than susceptible mosquitoes. The 1083 genes overexpressed in susceptible vs. resistant mosquitoes included genes associated with signal transduction, transport and proteolysis and were enriched in carbohydrate metabolic processes and cuticle (Fig. 3).

Fig. 3
figure 3

Functional classifications of candidate-resistance genes. Genes expressed significantly more (R > S) or less (R < S) in resistant versus susceptible mosquitoes were classified based on their functions. A percentage was attributed to each function based on the total number of genes considered. Functional abbreviation is as follows: UNK (Unknown); TR (transport); TT (transcription and translation); STD (signal transduction); RTS (response to stress); REDOX (oxido-reduction processes); PROT (proteolysis); OBP (odorant binding proteins); MCT (microtubule-associated movement); MET (metabolism); DNA_R (DNA repair); DIV (diverse functions); CHR (chromosome or chromatin-related functions); CUT (cuticule); CC (cell cycle); CA (catalytic activity)

Full size image

Detoxification genes

Thirty-nine detoxification genes were significantly DE, with 12 genes (HPX15 [AGAP013327], HPX7 [AGAP004036], HPX15 (AGAP013327), CYP304B1 [AGAP003066], CYP306A1 [AGAP004888], CYP315A1 [AGAP000284], CYP4C26 [AGAP000192], CYP4C27 [AGAP009246), CYP6M2, CYP6N1 [AGAP008210], CYP9J3 [AGAP012291], COEAE2G [AGAP006723]) showing more than 2 fold differential expression (Additional file 5). Among the detoxification genes previously associated with insecticide resistance [15, 16, 24, 50–52], CYPM2, CYP9J15 (AGAP012296), GSTT2 (AGAP009016) and GSTE5 (AGAP009192) were found to be 2.12, 1.78, 1.57 and 1.54 fold expressed more in resistant than susceptible mosquitoes, respectively; while CYP325C2 (AGAP002205) was expressed 2.23 fold more in susceptible than resistant mosquitoes, contrary to what was recently detected in An. coluzzii mosquitoes from West Africa [24].

Cuticular protein genes

A total of 64 cuticular protein genes were significantly DE (Additional file 4). The majority (95.31 %) of these genes was more highly expressed in susceptible than resistant mosquitoes, a trend that is consistent with previous results [18, 25]. The RNA-seq based lower expression in resistant than susceptible mosquitoes of two cuticular protein genes (CPLG3 [AGAP008446] and CPLG4 [AGAP008447]), previously likned to insecticide resistance in An. gambiae/An. coluzzii mosquitoes from West Africa [46, 53], was lower in resistant than susceptible mosquitoes in the RNA-seq analysis. This was confirmed by qRT-PCR on independent samples of phenotypic resistant and susceptible mosquitoes (Fig. 1c).

Comparison with the Kisumu susceptible strain

To account for induction of gene expression during the insecticide exposure bioassay, we filtered the candidate-resistance genes by analyzing their expression profiles in mosquitoes of the Kisumu strain. We assumed constitutive candidate-resistance genes are significantly DE between field-resistant mosquitoes and mosquitoes of the Kisumu strain, but not between field-susceptible and Kisumu mosquitoes because resistance to insecticide is considered a pre-adaptive phenotype [54]. A total of 702 genes were significantly DE between resistant mosquitoes and mosquitoes of the Kisumu strain, while 12467 genes (93.30 % of the total number of tested genes teste) were not DE between susceptible mosquitoes from the field and from the Kisumu strain. This filtering approach reduced the number of candidate-resistance genes to 182 (Additional file 4), including 105 genes expressed more highly in susceptible mosquitoes and enriched in functions such as proteolysis, organic acid metabolic processes, transport and cuticle; and 78 genes expressed more highly in resistant mosquitoes and functionally related to endopeptidase activity, cytochrome P450s and nucleotide binding (Fig. 3).

Comparison to previously detected candidate-resistance genes

We compared the 182 DE genes with DE genes previously detected from mosquitoes from Emutete, a rural town in the Western Province of Kenya, in 2010 [25]. A total of 55 common DE genes were identified, including seven overexpressed in resistant versus susceptible mosquitoes and linked to transferase activity and metabolic/detoxification processes (Additional file 5). Genes expressed more highly in susceptible than resistant mosquitoes were enriched in functions such as proteolysis and cuticle and included HPX2 (AGAP009033) and six cuticle protein encoding genes.

A list of 295 candidate resistance genes has been recently identified in An. coluzzii mosquitoes from Burkina Faso and Côte d’Ivoire by comparing deltamethrin-resistant field-derived mosquitoes and mosquitoes of two susceptible strains [24]. When a similar comparison is done across our RNA-seq data, a total of 21 candidate resistance genes was identified (Additional file 6), including two genes (AGAP004177, AGAP004572) mapping within a QTL previously linked to resistance to pyrethroid [55] and associated with RNA methylation and lipid metabolism, respectively (Additional file 6). These two genes and three others (AGAP013036, AGAP007530 and AGAP008840) are also DE when considering the comparison to susceptible mosquitoes from Western Kenya (Table 2). qRT-PCR analyses confirmed changes in expression profile as detected by RNA-seq for AGAP013036, AGAP008840 and AGAP004177; AGAP004572 showed higher expression in resistant versus susceptible samples. Within samples variation tended to be large, suggesting none of these genes has a major effect, but indirect and/or additive effects on the resistant phenotype are more probable (Additional file 7).

Table 2 Candidate resistance genes
Full size table

SNPs identified in RNA-seq libraries

SNPs in the pyrethroid target site

RNA-seq analysis of the para gene confirmed the presence and absence of the L1014S and N1575Y mutations, respectively. An additional change from A to G was identified at nucleotide position 2391280 in the 9th exon of the para gene, leading to a change from lysine [K] to arginine [R]) at amino acid position 419, based on the Musca domestica Vss1 gene nomenclature [Accession N. AAB47604] (Additional file 8). The para gene codes for a protein with four homologous domains, each composed of six segments. Both the 419 and the 1014 mutations lie in the 6th segment of the first and second domains, respectively [56]. While various mutations at position 1014 are wide-spread in insects, this is the first report of a non-synonymous mutation at position 419. Genotyping data on 96 Kenyan mosquitoes found the K419R mutation only in insecticide-susceptible mosquitoes.

SNPs in other genes: non-synonymous SNP potentially associated with insecticide-resistance in the sulfotransferase encoding gene AGAP009551

We searched for SNPs across all genes, including those not DE between resistant and susceptible mosquitoes, and identified 310 genes with SNP coverage only in resistant mosquitoes; in 102 of these genes susceptible mosquitoes showed coverage only for the not wild-type mutant (reference) nucleotide base (Additional file 9). The majority (92.26 %) of SNPs was associated with non-synonymous mutations, followed by start gain/loss (7.42 %) and splice acceptor/donors (0.32 %) (Additional file 9). A total of 39 SNPs were chosen for further genotyping in single mosquitoes based on: 1) their location within previously identified pyrethroid resistance QTLs [55], or 2) RNA-seq coverage (Additional files 1 and 9). In general, SNPs associated with pyrethroid resistance were rare. We found a change from G to C at position 36043997 of chromosome 3R that results in an Alanine to Glycine substitution at codon 1010 of AGAP009551, and was significantly associated with the resistance phenotype (odds ratio [95 %]: 5.10 [1.30–19.99]) (Additional file 8).

Discussion

In this study, we compared the transcriptome of pyrethroid-resistant and susceptible An. gambiae mosquitoes from Western Kenya and we analyzed results using a step-wise approach starting from comparison among local mosquitoes and finishing with a comparison across An. gambiae and An. coluzzii originating from different ecosystems. We identified 5 genes consistently DE between resistant and susceptible An gambiae/An. coluzzii across Sub-Saharan Africa and 1 SNP strongly associated with the resistance phenotype.

Markers for resistance surveillance

WHO standard bioassays are generally applied to assess phenotypic resistance of mosquito populations [7]. However, these methods require a large number of field-caught mosquitoes and variation in the age and physiological conditions of the specimens can affect the consistency of bioassay results. Moreover, standardization of any bioassay method is often difficult across sentinel sites because of variations in local temperature, humidity, and testing conditions, which significantly confound the bioassay results. Historically, knowledge of the mechanisms of resistance led to the identification of markers for resistance. Most of the markers used to for monitoring resistance to pyrethroids are DNA markers [14, 54, 57]. These DNA markers include the kdr mutation L1014S/F and an additional mutation in the para gene that was recently identified (N1575Y) [14, 54, 57]. Markers based on differential expression between resistant and susceptible mosquitoes (RNA markers) are emerging, primarily for metabolic detoxification genes [24]. The identification of RNA markers is technically challenged by 1) the limited stability of RNA and 2) the fact that genes known to be associated with insecticide resistance (i.e. metabolic detoxification genes and cuticular protein genes) belong to large gene families.

WHO guidelines define the resistant and susceptible phenotypes as that of mosquitoes alive or dead twenty-four hours after insecticide exposure, respectively [7]. Due to RNA degradation in dead specimens [58], a number of alternative strategies have been used to sample RNA from susceptible mosquitoes, including analyzing the expression profile of mosquitoes of susceptible populations/strains [21, 24, 50], field-caught unexposed mosquitoes [24] or approximating susceptibility with early knock-down [25, 28]. Each strategy has its advantages and disadvantages. For instance, susceptible laboratory strain mosquitoes usually do not share the same genetic background as field-caught mosquitoes and are expected to have lower genetic variability than field-caught mosquitoes. However, for the same reasons, there is less within-sample variation in the gene expression profiles of laboratory strains. Using field-caught mosquitoes has the advantage of limiting the impact of environmental conditions on a heterogeneous genetic background. However, unexposed field-caught mosquitoes are a mixture of resistant and susceptible individuals, and the expression profile of early knocked-down mosquitoes will include both constitutively expressed and insecticide-induced DE genes. As a consequence, a comparison across studies to find common RNA markers may be influenced by the chosen experimental design. Additionally, DE genes may have a direct effect on the resistance phenotype or an indirect or additive effect and may belong to large gene families [59, 60]. Gene products with functions such as metabolic detoxification, proteolysis and cuticular proteins have been linked to resistance to insecticides [14–21, 24, 25, 50, 53]. However, all of these proteins are members of large gene families, and it is unclear whether one/few members of each family have a major role in resistance [54] or whether functional convergence and additive effects contribute to the resistance phenotype. Functional convergence with respect to pyrethroid metabolism among the 111 annotated CYPs from the An. gambiae genome could explain variation across populations in CYP-encoding genes identified as significantly over-expressed in resistant mosquitoes through qRT-PCR, microarray and RNA-seq approaches [15, 21, 24, 25, 50, 54, 61]. If multiple CYPs are able to metabolize pyrethroids [16, 17], a mosquito can become resistant when certain CYPs act synergistically without invoking significant changes in the expression of individual CYPs. Additionally, some CYPs may show tissue-specific expression [62]. Alternatively, local adaptation may lead to over-expression of different CYPs in various geographic populations. From the operational standpoint of resistance monitoring, the observed variation across populations in the over-expression of CYP-encoding genes supports the possibility that sensitive and field-deployable biochemical assays that quantify overall monooxygenase activity [63, 64], may prove suitable for predicting the resistance status of a mosquito population [65].

Candidate-resistance genes

Genes related to transcription and translation and response to stress were enriched among the 2457 genes that were DE between field-caught resistant and susceptible mosquitoes, suggesting that insecticide exposure is a stress for mosquitoes, leading to complex gene interaction mechanisms, including gene expression modulation. A comparison between these 2457 DE genes and the 182 constitutive DE genes found six detoxification genes constitutively and more highly expressed in resistant mosquitoes, including the previously identified CYPM2 and GSTE5 [15, 17, 24, 50], CYP303A1 (AGAP10077), CYP4C27 (AGAP009246), GSTD3 (AGAP004382) and HPX2. Fold-changes in differential expression ranged between 1.55 and 2.33, which suggests an additive effect on the resistant phenotype is more probable than a major effect by the product of one gene. An additive effect of multiple detoxification genes and/or functional convergence would explain the lack of detoxification genes among the 7 DE genes consistently and more highly expressed in resistant versus susceptible mosquitoes from Western Kenya in a span of 2 years (2010–2012).

Genes associated with proteolytic functions and cuticle appeared to be enriched among the constitutive genes expressed more in susceptible than resistant mosquitoes. Most of these cuticular protein genes were also more highly expressed in susceptible than resistance mosquitoes in 2010 [25], which is different than what was seen in An. gambiae and An. coluzzii mosquitoes from West Africa [24, 64]. Whether this result is related to the different ecosystem and more intense agricultural activities in West Africa in comparison to the Western Province of Kenya needs to be investigated. Areas near the equator in West and East Africa have a different seasonality. In West Africa, a long dry season is followed by a single intense rainy season between July and September [67]. In Kenya, two rainy seasons generally occur between April-August and November-December [67]. The longer dry season in West Africa may precondition mosquitoes with a thicker cuticle that could result in penetration resistance to insecticides.

When compared to a list of candidate-resistance genes recently identified from An. coluzzii mosquitoes from West Africa [24], 5 genes (AGAP004177, AGAP007530, AGAP008840, AGAP013036 and AGAP004572) showed similar significant differential expression, emphasizing their potential role as markers for resistance for closely related An. gambiae/An. coluzzii mosquitoes across different ecosystems. AGAP004177 and AGAP004572 are particularly interesting because they map within rtp1, a QTL previously associated with resistance to pyrethroid and including the para sodium channel gene [55]. AGAP004177 encodes a heat shock protein with 23S rRNA (uridine2552–2′-O)-methyltransferase activity [68]. Incorrect ribosomal methylation and/or impairment in ribosome maturation and function have been linked to drug resistance in prokaryotes [69, 70], suggesting the activity of this gene should be further investigated. AGAP004572 encodes a hypothetical protein with a fatty acid desaturase domain, based on the Eukaryotic Orthologous Groups (KOG) [71] and SMART databases [72]. Fatty acid desaturases are conserved proteins that create a double bond in long-chain fatty acids, the primary determinant of triglyceride melting temperature and cellular membrane fluidity [73]. Fatty acid desaturase activity and overexpression of cuticular proteins may provide an additional mechanism for insecticide penetration. AGAP013036 encodes a hypothetical protein, which the KOG database identifies as a probable phosphatidylethanolamine-binding protein (PEBP). PEBPs are highly conserved and have been associated with diverse functions such as lipid binding, signal transduction, serine protease inhibition and neuronal development [74], suggesting an indirect role in resistance. AGAP007530 encodes a hypothetical protein with a predicted zinc carboxypeptidase domain, based on KOG and SMART databases. In prokaryotes, zinc carboxypeptidases may be used to digest exogenous proteins and degrade tissues [75, 76], suggesting their potential role in insecticide catabolism should be investigated. AGAP008840 encodes a hypothetical protein linked to chromatin structure and dynamics. Its role in insecticide resistance is unclear.

SNPs associated with the resistance phenotype

The para sodium channel is the target site for pyrethroid. Electrophysiological studies of this protein from Drosophila melanogaster, Aedes aegypti, Musca domestica and Blatella germanica expressed in Xenopus oocytes showed that the L1014F/S mutation alters its function and provides protection from pyrethroids [56]. The frequency of the L1014S/F mutation increased significantly in the past 10 years across sub-Saharan Africa in both An. gambiae and An. coluzzii, concomitantly with the implementation of malaria control strategies [8, 9]. This coincidence was used to support the conclusion that there is causal relationship between L1014S/F mutations at the para sodium channel and the resistant phenotype [56]. In our experiments, RNA-seq data and genotyping in individually-phenotyped mosquitoes showed the L1014S mutation in both insecticide-resistant and susceptible An. gambiae mosquitoes, confirming that the L1014S mutation alone cannot fully account for the resistant phenotype [77]. Recent genotyping data showed that the L1014F does not account for high levels of resistance in An. coluzzii mosquitoes from Burkina Faso [24] nor it is linked with the resistant phenotype in An. gambiae and An. coluzzii mosquitoes from Nigeria and Benin [64, 66]. These results suggest that additional mechanisms of resistance have evolved, including novel mutations, such as N1575Y [14], and/or changes to gene expression. For instance, if metabolic detoxification mechanisms or penetration resistance are prevailing [66], the amount of insecticide available to bind the insecticide target will be reduced, explaining the observed limited correlation between L1014F/S and the resistant phenotype [14, 66]. In this scenario, the difference between resistant and susceptible mosquitoes could result from several possibilities, including barriers to insecticides penetration, catabolic pathways employed by resistant and susceptible mosquitoes, secondary pathways involved in the metabolism of insecticide intermediates and/or a more efficient “pyrethrome” [78] in resistant mosquitoes resulting in a quicker insecticide degradation and/or sequestration and elimination of toxic secondary metabolites.

RNA-seq data, followed by genotyping on singly phenotyped mosquitoes, identified the mutation A101G in the gene AGAP009551, which was associated with the resistant phenotype and should be further investigated. AGAP009551 is annotated as a sulfotransferase and maps within rtp3 [55]. Sulfotransferases catalyze the sulfate conjugation of hormones and xenobiotic compounds in humans and plants [79]. Sulfonation can lead to either enhanced secretion of the compounds or the production of bioactive metabolites [79]. Primary catabolic pathways identified for deltamethrin include the CYP-mediated hydroxylation of deltamethrin to 4′-hydroxydeltamthrin [17]; the final metabolites in the pathway include trans-hydroxymethyl deltamethrin and deltamethric acid. Alternatively, cleavage of the ester bond by carboxylesterases renders the fragment cyclopropane carboxylic acid and phenoxybenzylic alcohol inactive as insecticides [79]. Sulfotransferase-mediated conjugation could be an alternative degradation route for deltamehtrin, or act on the metabolites produced after oxidation, reduction or hydrolysis in Phase I of the degradation process. Sulfotransferases target electrophilic centers favoring conjugation to cellular sugars, glutathione or amino acids (biotransformation, Phase II of the degradation process) [12].

Conclusions

By comparing the expressional profile of deltamethrin-resistant and –susceptible mosquitoes from Western Kenya, we identified five genes with a similar expression profile in resistant An. gambiae and An. coluzzii mosquitoes from West and East Africa, suggesting these genes could be used as expression-based markers for resistance. We also identified a SNP in the sulfotraferase gene AGAP009551 that is strongly associated with insecticide resistance. No correlation with the resistant phenotype was seen for the kdr mutation L1014S.

Overall, our results support the idea that many genes may be involved synergistically in insecticide resistance, with their role being modulated by the life history and selection pressure of mosquito populations.

References

  1. The Global Fund [http://www.theglobalfund.org/en/about/diseases/malaria/] Accessed April 21st 2015.

  2. President’s Malaria Initiative [http://www.pmi.gov/] Accessed April 21st 2015.

  3. WHO. World malaria report. Geneva: World Health Organization; 2013.

    Google Scholar 

  4. Cohen JM, Smith DL, Cotter C, Ward A, Yamey G, Sabot OJ, et al. Malaria resurgence: a systematic review and assessment of its causes. Malar J. 2012;11:12.

    Article  Google Scholar 

  5. WHO. The Technical basis for coordinated action against insecticide resistance: preserving the effectiveness of modern malaria vector control. Geneva: World Health Organization; 2011.

    Google Scholar 

  6. WHO. Global plan for insecticide resistance management in malaria vectors (GPIRM). Geneva: World Health Organization; 2012.

    Google Scholar 

  7. WHO. Test procedures for insecticide resistance monitoring in malaria vector mosquitoes. Geneva: World Health Organization; 2013.

    Google Scholar 

  8. Ranson H, N’guessan R, Lines J, Moiroux N, Nkuni Z, Corbel V. Pyrethroid resistance in African anopheline mosquitoes: what are the implications for malaria control? Trends Parasitol. 2011;27:91–8.

    Article  CAS  PubMed  Google Scholar 

  9. Knox TB, Juma EO, Ochomo EO, Jamet HP, Ndungo L, Chege P, et al. An online tool for mapping insecticide resistance in major Anopheles vectors of human malaria parasites and review of resistance status for the Afrotropical region. Parasit Vectors. 2014;7:76.

    Article  PubMed Central  PubMed  Google Scholar 

  10. Lanzaro GC, Lee Y. Speciation in Anopheles gambiae – The Distribution of Genetic polymorphism and patterns of reproductive isolation among natural populations. In S. Manguin, editor. Anopheles mosquitoes – new insights into malaria vectors. INTECH 2013. DOI: 10.5772/56232.

  11. Norris LC, Main BJ, Lee Y, Collier TC, Fofana A, Cornel AJ, et al. Adaptive introgression in an African malaria mosquito coincident with the increased usage of insecticide-treated bed nets. Proc Natl Acad Sci U S A. 2015;112:815–20.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  12. Khambay BPS, Jewess PJ. Pyrethroids. In: Gilbert LI, Gill SS, editors. Insect control biological and synthetic agents. Oxford, UK: Elsivier; 2010. p. 1–29.

    Google Scholar 

  13. Martinez-Torres D, Chandre F, Williamson MS, Darriet F, Berge JB, Devonshire AL, et al. Molecular characterization of pyrethoid knowdown resistance (kdr) in the major malaria vector Anopheles gambiae s.s. Insect Mol Biol. 1998;7:179–84.

    Article  CAS  PubMed  Google Scholar 

  14. Jones CM, Liyanapathirana M, Agossa FR, Weetman D, Ranson H, Donnelly MJ, et al. Footprints of positive selection associated with a mutation (N1575Y) in the voltage-gated sodium channel of Anopheles gambiae. Proc Natl Acad Sci U S A. 2012;109:6614–9.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  15. David JP, Ismail HM, Chandor-Proust A, Paine MJI. Role of cytochrome P450s in insecticide resistance impact on the control of mosquito-borne diseases and use of insecticides on Earth. Phil Trans R Soc B. 2013;368:20120429.

    Article  PubMed Central  PubMed  Google Scholar 

  16. Müller P, Warr E, Stevenson BJ, Pignatelli PM, Morgan JC, Steven A, et al. Field-caught permethrin-resistant Anopheles gambiae overexpress CYP6P3, a P450 that metabolises pyrethroids. PLoS Genet. 2008;4, e1000286.

    Article  PubMed Central  PubMed  Google Scholar 

  17. Stevenson BJ, Bibby J, Pignatelli P, Muangnoicharoen S, O’Neill PM, Lian LY, et al. Cytochrome P450 6M2 from the malaria vector Anopheles gambiae metabolizes pyrethroids: sequential metabolism of deltamethrin revealed. Insect Biochem Mol Biol. 2011;41:492–502.

    Article  CAS  PubMed  Google Scholar 

  18. Reid WR, Zhang L, Liu F, Liu N. The transcriptome profile of the mosquito Culex quinquefasciatus following permethrin selection. PLoS One. 2012;7:e47163.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  19. Vontas JG, Small GJ, Hemingway J. Glutathione S-transferasesas antioxidant defence agents confer pyrethroid resistance in Nilaparvata lugens. Biochem J. 2001;357:65–72.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  20. Nkya TE, Akhouayri I, Kisinza W, David JP. Impact of environment on mosquito response to pyrethroid insecticides: facts, evidences and prospects. Insect Biochem Mol Biol. 2013;43:407–16.

    Article  CAS  PubMed  Google Scholar 

  21. Nkya TE, Akhouayri I, Poupardin R, Batengana B, Mosha F, Magesa S, et al. Insecticide resistance mechanisms associated with different environments in the malaria vector Anopheles gambiae: a case study in Tanzania. Malar J. 2014;13:28.

    Article  PubMed Central  PubMed  Google Scholar 

  22. Wood O, Hanrahan S, Coetzee M, Koekemoer L, Brooke B. Cuticle thickening associated with pyrethroid resistance in the major malaria vector Anopheles funestus. Parasit Vectors. 2010;3:67.

    Article  PubMed Central  PubMed  Google Scholar 

  23. Gatton ML, Chitnis N, Churcher T, Donnelly MJ, Ghani AC, Godfray HC, et al. The importance of mosquito behavioural adaptations to malaria control in Africa. Evolution. 2013;67:1218–30.

    Article  PubMed Central  PubMed  Google Scholar 

  24. Toe KH, N’Fale S, Dabire RK, Ranson H, Jones CM. The recent escalation in strength of pyrethroid resistance in Anopheles coluzzi in West Africa is linked to increased expression of multiple gene families. BMC Genomics. 2015;16:146.

    Article  PubMed Central  PubMed  Google Scholar 

  25. Bonizzoni M, Afrane Y, Dunn WA, Atieli FK, Zhou G, Xhong D, et al. Comparative Transcriptome Analyses of Deltamethrin-Resistant and –Susceptible Anopheles gambiae Mosquitoes from Kenya by RNA-Seq. PLoS One. 2012;7, e44607.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  26. Vulule JM, Beach RF, Atieli FK, Roberts JM, Mount DL, Mwangi RW. Reduced susceptibility of Anopheles gambiae to permethrin associated with the use of permethrin-impregnated bednets and curtains in Kenya. Med Vet Entomol. 1994;8:71–5.

    Article  CAS  PubMed  Google Scholar 

  27. Chuaycharoensuk T, Juntarajumnong W, Boonyuan W, Bangs MJ, Akratanakul P, Thammapalo S, et al. Frequency of pyrethroid resistance in Aedes aegypti and Aedes albopictus (Diptera: Culicidae) in Thailand. J Vector Ecol. 2011;36:204–12.

    Article  PubMed  Google Scholar 

  28. Zhong D, Chang X, Zhou G, He Z, Fu F, Yan Z, et al. Relationship between knockdown resistance, metabolic detoxification and organismal resistance to pyrethroids in Anopheles sinensis. PLoS One. 2013;8, e55475.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  29. Balmert NJ, Rund SS, Ghazi JP, Zhou P. Duffield GE. Time-of-day specific changes in metabolic detoxification and insecticide resistance in the malaria mosquito Anopheles gambiae. J Insect Physiol. 2014;64:30–9.

    Article  CAS  PubMed  Google Scholar 

  30. Trizol-based RNA extraction protocol [http://medicine.yale.edu/keck/ycga/microarrays/protocols/TRIZOLRNAIsolation_092107_tcm240-21453_tcm240-284-32.pdf] Accessed June 25th 2015.

  31. Scott JA, Brogdon W, Collins FH. Identification of single specimens of the Anopheles gambiae complex by the polymerase chain reaction. Am J Trop Med Hyg. 1993;49:520.

    CAS  PubMed  Google Scholar 

  32. DNA Technologies and Expression Analysis Core, UC Davis Genome Center [http://expression.genomecenter.ucdavis.edu/].

  33. Bioconductor package qrqc [http://bioconductor.org/packages/release/bioc/html/qrqc.html].

  34. Scythe v.0.990 and Sickle v.1.200 [https://github.com/ucdavis-bioinformatics].

  35. Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, et al. Differential gene and transcript expression analyses of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012;7:562.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  36. Blacktie [https://github.com/xguse/blacktie].

  37. CommeRbund [http://www.bioconductor.org/packages/release/bioc/html/cummeRbund.html].

  38. VectorBase [www.vectorbase.org].

  39. AnoXcel [http://www.niaid.nih.gov/LabsAndResources/labs/aboutlabs/lmvr/Pages/TranscriptomeResources.aspx#transcriptomes %20].

  40. G:Profiler [http://biit.cs.ut.ee/gprofiler/].

  41. Freebayes v. 0.9.4 [http://bioinformaticsbcedu/marthlab/FreeBayes].

  42. SnpEff [http://snpeff.sourceforge.net/].

  43. Bonizzoni M, Britton M, Marinotti O, Dunn WA, Fass J, James AA. Probing functional polymorphisms in the dengue vector, Aedes aegypti. BMC Genomics. 2013;14:739.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  44. Schlotterer C, Tobler R, Kofler R, Nolte V. Sequencing pools of individuals –mining genome wide polymorphism data without big funding. Nat Rev Genet. 2014;15:749–63.

    Article  PubMed  Google Scholar 

  45. Paris M, Marcombe S, Coissac E, Corbel V, David JP, Despres L. Investigating the genetics of Bti resistance using mRNA tag sequencing: application on laboratory strains and natural populations of the dengue vector Aedes aegypti. Evol Appl. 2013;6:1012.

    PubMed Central  CAS  PubMed  Google Scholar 

  46. Vannini L, Reed TW, Willis JH. Temporal and spatial expression of cuticular proteins of Anopheles gambiae implicated in insecticide resistance or differentiation of M/S incipient species. Parasit Vectors. 2014;7:24.

    Article  PubMed Central  PubMed  Google Scholar 

  47. Pearson Correlation Coefficient Calculator [http://www.alcula.com/calculators/statistics/correlation-coefficient/]

  48. Sanger F, Coulson AR. A rapid method for determining sequences in DNA by primed synthesis with DNA polymerase. J Mol Biol. 1975;94:441.

    Article  CAS  PubMed  Google Scholar 

  49. Odds Ratio Calculator [http://vassarstats.net/odds2x2.html].

  50. Djouaka RF, Bakare AA, Coulibaly ON, Akogbeto MC, Ranson H, Hemingway J, et al. Expression of the cytochrome P450s, CYP6P3 and CYP6M2 are significantly elevated in multiple pyrethroid resistant populations of Anopheles gambiae s.s. from Southern Benin and Nigeria. BMC Genomics. 2008;9:538.

    Article  PubMed Central  PubMed  Google Scholar 

  51. Matowo J, Jones CM, Kabula B, Ranson H, Steen K, Mosha F, et al. Genetic basis of pyrehtroid resistance in a population of Anopheles arabiensis, the primary malaria vector in Lower Moshi, north-eastern Tanzania. Parasit Vectors. 2014;7:274.

    Article  PubMed Central  PubMed  Google Scholar 

  52. Wilding CS, Weetman D, Rippon EJ, Steen K, Mawejje HD, Barsukov I, et al. Parallel evolution or purifying selection, not introgression, explains similarity in the pyrethroid detoxification linked GSTE4 of Anopheles gambiae and An. arabiensis. Mol Genet Genomics. 2015;290:201–15.

    Article  CAS  PubMed  Google Scholar 

  53. Awolola TS, Oduola OA, Strode C, Koekemoer LL, Brooke B, Ranson H. Evidence of multiple pyrethroid resistance mechanisms in the malaria vector Anopheles gambiae sensu stricto from Nigeria. Trans R Soc Trop Med Hyg. 2009;103:1139–45.

    Article  CAS  PubMed  Google Scholar 

  54. Liu N. Insecticide resistance in mosquitoes: impact, mechanisms, and research directions. Annu Rev Entomol. 2015;60:537–59.

    Article  CAS  PubMed  Google Scholar 

  55. Ranson H, Paton MG, Jensen B, McCarroll L, Vaughan A, Hogan JR, et al. Genetic mapping of genes conferring permethrin resistance in the malaria vector Anopheles gambiae. Insect Mol Biol. 2004;13:379–86.

    Article  CAS  PubMed  Google Scholar 

  56. Dong K, Du Y, Rinkevich F, Nomura Y, Xu P, Wang L, et al. Molecular biology of insect sodium channels and pyrethroid resistance. Insect Biochem Mol Biol. 2014;50:1–17.

    Article  PubMed Central  PubMed  Google Scholar 

  57. Weetman D, Donnely M. Evolution of insecticide resistance diagnosis in malaria vectors. Trans R Soc Trop Med Hyg. 2015;109:291–3.

    Article  PubMed  Google Scholar 

  58. Gallego Romero I, Pai AA, Thung J, Gilad Y. RNA-seq: impact of RNA degradation on transcript quantification. BMC Biol. 2014;12:42.

    Article  PubMed Central  PubMed  Google Scholar 

  59. Li T, Liu L, Zhang L, Liu N. Role of G-protein-coupled receptor-related genes in insecticide resistance of the mosquito, Culex quinquefasciatus. Sci Rep. 2014;4:6474.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  60. Gong Y, Li T, Zhang L, Gao X, Liu N. Permetrhin induction of multiple cytochrome P450 genes in insecticide resistant mosquitoes, Culex quniquefasciatus. Int J Biol Sci. 2013;9:863–71.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  61. Nkya TE, Poupardin R, Laporte F, Akhouayri I, Mosha F, Magesa S, et al. Impact of agriculture on the selection of insecticide resistance in the malaria vector Anopheles gambiae: a multigenerational study in controlled conditions. Parasit Vectors. 2014;7:480.

    PubMed Central  PubMed  Google Scholar 

  62. Ingham VA, Jones CM, Pignatelli P, Balabanidou V, Vontas J, Wagstaff SC, et al. Dissecting the organ specificity of insecticide resistance candidate genes in Anopheles gambiae: known and novel candidate genes. BMC Genomics. 2014;15:1018.

    Article  PubMed Central  PubMed  Google Scholar 

  63. Chen H, Githeko AK, Githure JI, Mutunga J, Zhou G, Yan G. Monooxygenase levels and knockdown resistance (kdr) allele frequencies in Anopheles gambiae and Anopheles arabiensis in Kenya. Med Vet Entomol. 2008;45:242–50.

    Article  CAS  Google Scholar 

  64. Gnanguenon V, Agossa FR, Badirou K, Govoetchan R, Anagonou R, Oke-Agbo F, et al. Malaria vectors resistance to insecticides in Benin: current trends and mechanisms involved. Parasit Vectors. 2015;8:223.

    Article  PubMed Central  PubMed  Google Scholar 

  65. Chouaïbou M, Zivanovic GB, Knox TB, Jamet HP, Bonfoh B. Synergist bioassays: a simple method for initial metabolic resistance investigation of field Anopheles gambiae s.l. Populations. Acta Trop. 2013;130C:108–11.

    PubMed  Google Scholar 

  66. Okorie PN, Ademowo OG, Irving H, Kelly-Hope LA, Wondji CS. Insecticide susceptibility of Anopheles coluzzii and Anopheles gambiae mosquitoes in Ibadan, Southwest Nigeria. Med Vet Entomol. 2015;29:44–50.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  67. https://courseware.e-education.psu.edu/courses/earth105new/content/lesson07/03.html, accessed May 1st 2015.

  68. Hager J, Staker BL, Bugl H, Jakob U. Active site in Rrmj, a heat shock-induced methyltransferase. J Biol Chem. 2002;277:41978–86.

    Article  CAS  PubMed  Google Scholar 

  69. Vazquez-Laslop N, Thum C, Mankin AS. Molecular mechanism of drug-dependent ribosome stalling. Mol Cell. 2008;30:190–202.

    Article  CAS  PubMed  Google Scholar 

  70. Monshupanee T, Johansen SK, Dahlberg AE, Douthwaite S. Capreomycin susceptibility is increased by TlyA-directed 2′-O-methylation on both ribosomal subunits. Mol Microbiol. 2012;85:1194–11203.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  71. Eukaryotic Orthologous Groups (KOG) [http://genome.jgi-psf.org/help/kogbrowser.jsf].

  72. SMART databases [http://smart.embl-heidelberg.de/].

  73. Nakamura MT, Nara TY. Structure, function, and dietary regulation of Δ6, Δ5, and Δ9 desaturases. Annu Rev Nutr. 2004;24:345–76.

    Article  CAS  PubMed  Google Scholar 

  74. Vallee BS, Coadou G, Labbe H, Sy D, Vovelle F, Schoentgen F. Peptides corresponding to the N- and C-terminal parts of PEBP are well-structured in solution: new insights into their possible interaction with partners in vivo. J Pept Res. 2003;61:47–57.

    Article  CAS  PubMed  Google Scholar 

  75. Kooi C, Sokol PA. Differentiation of thermolysins and serralysins by monoclonal antibodies. J Med Microbiol. 1996;45:219–25.

    Article  CAS  PubMed  Google Scholar 

  76. De Kreij A, Venema G, van den Burg B. Substrate specificity in the highly heterogeneous M4 peptidase family is determined by a small subset of amino acids. J Biol Chem. 2000;275:31115–20.

    Article  PubMed  Google Scholar 

  77. Donnelly MJ, Corbel V, Weetman D, Wilding CS, Williamnson MS. Black WC 4th: Does kdr genotype predict insecticide-resistance phenotype in mosquitoes? Trends Parasitol. 2009;25:213-219.

  78. Ismail HM, O'Neill PM, Hong DW, Finn RD, Henderson CJ, Wright AT. Pyrethroid activitybased probes for profiling cytochrome P450 activities associated with insecticide interactions. Proc Natl Acad Sci U S A. 2013;110:19766-19771.

  79. Weinshilboum RM, Otterness DM, Aksoy IA, Wood TC, Her C, Raftogianis RB. Sulfation and sulfotransferases 1: Sulfotransferase molecular biology: cDNAs and genes. FASEB J. 1997;11:3-14.

Download references

Acknowledgements

We thank two anonymous reviewers for their constructive comments. The project was supported by National Institutes of Health grants R21 AI098652-01, R01 AI050243 and D43 TW001505. This work was made possible, in part, through access to the Genomic High Throughput Facility Shared Resource of the Cancer Center Support Grant (CA-62203) at the University of California, Irvine. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Author information

Authors and Affiliations

  1. Program in Public Health, University of California, Irvine, CA, USA

    Mariangela Bonizzoni, Guofa Zhou, Ming-Chieh Lee, Jiabao Xu & Guiyun Yan

  2. Department of Biology and Biotechnology, University of Pavia, Pavia, Italy

    Mariangela Bonizzoni

  3. Centre for Global health Research, Kenya Medical Research Institute, Kisumu, Kenya

    Eric Ochomo & Andrew Githeko

  4. Department of Ecology and Evolutionary Biology, Yale University, New Haven, CT, USA

    William Augustine Dunn

  5. Bioinformatics Core of the UC Davis Genome Center, University of California, Davis, CA, 95616, USA

    Monica Britton & Joseph Fass

  6. School of Health Sciences, Jaramogi Oginga Odinga University of Science and Technology, Bondo, Kenya

    Yaw Afrane

  7. Delta-9 Technologies, San Diego, CA, 92101, USA

    Joshua Hartsel

Authors
  1. Mariangela Bonizzoni
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Eric Ochomo
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. William Augustine Dunn
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Monica Britton
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Yaw Afrane
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Guofa Zhou
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Joshua Hartsel
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Ming-Chieh Lee
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Jiabao Xu
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Andrew Githeko
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Joseph Fass
    View author publications

    You can also search for this author in PubMed Google Scholar

  12. Guiyun Yan
    View author publications

    You can also search for this author in PubMed Google Scholar

Corresponding author

Correspondence to Mariangela Bonizzoni.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’contributions

MB and GY designed the study; MB, EO, YA, JX and WAD performed the experiments; MB, WAD, MB, GZ, JH, JF, AG, ML and JM analyzed the data. MB and JH drafted the manuscript. All authors read and approved the final manuscript.

Additional files

Additional file 1:

List of primers. List of genes/primers used in qRT-PCR (A) and in genotyping (B). (XLSX 15 kb)

Additional file 2:

Quality control of RNA-seq data: Box Plot of transcript quantification levels. The distributions of FPKM scores across samples are visualized. (PDF 130 kb)

Additional file 3:

Output of the Tuxedo pipeline on gene differential expression. List of genes analyzed, gene quantification levels (output of Cufflinks) and analyses of differential expression (output of Cuffdiff). (XLSX 6499 kb)

Additional file 4:

Lists of DE genes. List of the 2457 genes significantly DE between field-caught resistant and susceptible mosquitoes (Sheet 2457 RvS). List of the 182 genes significantly DE between field-caught resistant and susceptible mosquitoes of the Kisumu strain, but not field-caught susceptible mosquitoes and Kisumu mosquitoes (sheet 182 constitutive DE gene). List of the 55 DE genes in mosquitoes from Western Kenya in 2010 and 2012 (Sheet 55 candidate resistance genes). (XLSX 394 kb)

Additional file 5:

Detoxification genes. Fold changes in gene expression across detoxification-related genes. Genes previously associated with insecticide resistance are marked with a red arrow. (PDF 173 kb)

Additional file 6:

Candidate-resistance genes. List of candidate resistance genes detected in An. gambiae mosquitoes from West and East Africa. (XLSX 18 kb)

Additional file 7:

qRT-PCR on candidate resistance genes. The level of expression of 4 candidate resistance genes was measured by qRT-PCR from different resistant and susceptible mosquitoes than those used for RNA-seq. All tested genes were significantly differentially expressed between resistant and susceptible mosquitoes. (PDF 35 kb)

Additional file 8:

Novel SNP in the para gene. Visualization of the change from A to G that was identified at nucleotide position 2391280 of the para gene (AGAP004707) from RNA-seq data. Gray horizontal lines represent RNA-seq reads from libraries of both resistant and susceptible mosquitoes. (PDF 17 kb)

Additional file 9:

SNP-genes. List of genes with SNPs identified only in resistant mosquitoes from RNA-seq data and Odd ratio analyses for the genes that were genotyped in single mosquitoe. (XLSX 49 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Bonizzoni, M., Ochomo, E., Dunn, W.A. et al. RNA-seq analyses of changes in the Anopheles gambiae transcriptome associated with resistance to pyrethroids in Kenya: identification of candidate-resistance genes and candidate-resistance SNPs. Parasites Vectors 8, 474 (2015). https://doi.org/10.1186/s13071-015-1083-z

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13071-015-1083-z

Keywords

Parasites & Vectors

ISSN: 1756-3305

Contact us