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
© Bonizzoni et al. 2015
Received: 13 May 2015
Accepted: 10 September 2015
Published: 17 September 2015
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.
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.
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.
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.
KeywordsAnopheles gambiae insecticide resistance RNA-seq metabolic detoxification
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 . However, malaria still kills more than half a million people a year  and weakening of malaria control programs could favor malaria resurgence . 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 . 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 . 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 . The understanding of the mechanisms of insecticide resistance and the monitoring of resistance are two of the five pillars of this plan . Historically, the identification of the mechanisms of resistance has been important for developing molecular monitoring tools of resistance .
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 . 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 . 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 . The N1575Y mutation is thought to either enhance resistance due to the L1014F mutation or compensate for its fitness costs . 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 . 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) . 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 .
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 , 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.
Phenotypic resistance to deltamethrin
Resistance to deltamethrin was tested by the standard WHO tube test  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 . 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 .
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) . 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 . 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  and 1575-F (5′ TAAACAGCCTATACGGGAAACG) and 1575-R (5′ CGAGGAATTGCCTTTAGAGGTTTCT), respectively.
RNA-seq library preparation and sequencing
RNA-seq sample summary
Total aligned reads3
Western Province, Kenya
Raw read quality and contamination was evaluated using the Bioconductor package qrqc . Scythe v.0.990 and Sickle v.1.200  were used for Illumina adapter and quality-based trimming. Differential expression analyses followed the Tuxedo pipeline , which was run in Blacktie , 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 . 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  and AnoXcell . Functional enrichment analyses was done using g: Profiler .
SNPs were called with Freebayes v. 0.9.4 and SnpEff [41, 42], following a previously described pipeline . 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  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 , 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 .
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 . The odds ratio test was applied to determine if the odds of each analyzed SNP differed significantly between resistant and susceptible mosquitoes .
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 . 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
Differential expression of genes between field-caught resistant and susceptible mosquitoes
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 .
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 . 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 . 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.
Candidate resistance genes
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 . 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 , 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).
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 . 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 . 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 . Due to RNA degradation in dead specimens , 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  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  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 . 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 .
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 , 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 . In Kenya, two rainy seasons generally occur between April-August and November-December . 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 , 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 . AGAP004177 encodes a heat shock protein with 23S rRNA (uridine2552–2′-O)-methyltransferase activity . 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)  and SMART databases . 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 . 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 , 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 . 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 . 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 . Recent genotyping data showed that the L1014F does not account for high levels of resistance in An. coluzzii mosquitoes from Burkina Faso  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 , and/or changes to gene expression. For instance, if metabolic detoxification mechanisms or penetration resistance are prevailing , 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”  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 . Sulfotransferases catalyze the sulfate conjugation of hormones and xenobiotic compounds in humans and plants . Sulfonation can lead to either enhanced secretion of the compounds or the production of bioactive metabolites . Primary catabolic pathways identified for deltamethrin include the CYP-mediated hydroxylation of deltamethrin to 4′-hydroxydeltamthrin ; 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 . 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) .
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.
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.
Open AccessThis 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.
- The Global Fund [http://www.theglobalfund.org/en/about/diseases/malaria/] Accessed April 21st 2015.
- President’s Malaria Initiative [http://www.pmi.gov/] Accessed April 21st 2015.
- WHO. World malaria report. Geneva: World Health Organization; 2013.Google Scholar
- 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.View ArticleGoogle Scholar
- 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
- WHO. Global plan for insecticide resistance management in malaria vectors (GPIRM). Geneva: World Health Organization; 2012.Google Scholar
- WHO. Test procedures for insecticide resistance monitoring in malaria vector mosquitoes. Geneva: World Health Organization; 2013.Google Scholar
- 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.View ArticlePubMedGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.Google Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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
- 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.View ArticlePubMedGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.View ArticlePubMedGoogle Scholar
- Reid WR, Zhang L, Liu F, Liu N. The transcriptome profile of the mosquito Culex quinquefasciatus following permethrin selection. PLoS One. 2012;7:e47163.PubMed CentralView ArticlePubMedGoogle Scholar
- Vontas JG, Small GJ, Hemingway J. Glutathione S-transferasesas antioxidant defence agents confer pyrethroid resistance in Nilaparvata lugens. Biochem J. 2001;357:65–72.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.View ArticlePubMedGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.View ArticlePubMedGoogle Scholar
- 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.View ArticlePubMedGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.View ArticlePubMedGoogle Scholar
- 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.
- 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.PubMedGoogle Scholar
- DNA Technologies and Expression Analysis Core, UC Davis Genome Center [http://expression.genomecenter.ucdavis.edu/].
- Bioconductor package qrqc [http://bioconductor.org/packages/release/bioc/html/qrqc.html].
- Scythe v.0.990 and Sickle v.1.200 [https://github.com/ucdavis-bioinformatics].
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- Blacktie [https://github.com/xguse/blacktie].
- CommeRbund [http://www.bioconductor.org/packages/release/bioc/html/cummeRbund.html].
- VectorBase [www.vectorbase.org].
- AnoXcel [http://www.niaid.nih.gov/LabsAndResources/labs/aboutlabs/lmvr/Pages/TranscriptomeResources.aspx#transcriptomes %20].
- G:Profiler [http://biit.cs.ut.ee/gprofiler/].
- Freebayes v. 0.9.4 [http://bioinformaticsbcedu/marthlab/FreeBayes].Google Scholar
- SnpEff [http://snpeff.sourceforge.net/].
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.View ArticlePubMedGoogle Scholar
- 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 CentralPubMedGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- Pearson Correlation Coefficient Calculator [http://www.alcula.com/calculators/statistics/correlation-coefficient/]
- Sanger F, Coulson AR. A rapid method for determining sequences in DNA by primed synthesis with DNA polymerase. J Mol Biol. 1975;94:441.View ArticlePubMedGoogle Scholar
- Odds Ratio Calculator [http://vassarstats.net/odds2x2.html].
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.View ArticlePubMedGoogle Scholar
- 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.View ArticlePubMedGoogle Scholar
- Liu N. Insecticide resistance in mosquitoes: impact, mechanisms, and research directions. Annu Rev Entomol. 2015;60:537–59.View ArticlePubMedGoogle Scholar
- 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.View ArticlePubMedGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- Weetman D, Donnely M. Evolution of insecticide resistance diagnosis in malaria vectors. Trans R Soc Trop Med Hyg. 2015;109:291–3.View ArticlePubMedGoogle Scholar
- Gallego Romero I, Pai AA, Thung J, Gilad Y. RNA-seq: impact of RNA degradation on transcript quantification. BMC Biol. 2014;12:42.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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 CentralPubMedGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.PubMedGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- https://courseware.e-education.psu.edu/courses/earth105new/content/lesson07/03.html, accessed May 1st 2015.
- Hager J, Staker BL, Bugl H, Jakob U. Active site in Rrmj, a heat shock-induced methyltransferase. J Biol Chem. 2002;277:41978–86.View ArticlePubMedGoogle Scholar
- Vazquez-Laslop N, Thum C, Mankin AS. Molecular mechanism of drug-dependent ribosome stalling. Mol Cell. 2008;30:190–202.View ArticlePubMedGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- Eukaryotic Orthologous Groups (KOG) [http://genome.jgi-psf.org/help/kogbrowser.jsf].
- SMART databases [http://smart.embl-heidelberg.de/].
- Nakamura MT, Nara TY. Structure, function, and dietary regulation of Δ6, Δ5, and Δ9 desaturases. Annu Rev Nutr. 2004;24:345–76.View ArticlePubMedGoogle Scholar
- 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.View ArticlePubMedGoogle Scholar
- Kooi C, Sokol PA. Differentiation of thermolysins and serralysins by monoclonal antibodies. J Med Microbiol. 1996;45:219–25.View ArticlePubMedGoogle Scholar
- 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.View ArticlePubMedGoogle Scholar
- 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.Google Scholar
- 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.Google Scholar
- 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.Google Scholar