Skip to main content
Download PDF

Whole transcriptomic analysis reveals overexpression of salivary gland and cuticular proteins genes in insecticide-resistant Anopheles arabiensis from Western Kenya

BMC Genomics volume 25, Article number: 313 (2024) Cite this article

Abstract

Background

Effective vector control is key to malaria prevention. However, this is now compromised by increased insecticide resistance due to continued reliance on insecticide-based control interventions. In Kenya, we have observed heterogenous resistance to pyrethroids and organophosphates in Anopheles arabiensis which is one of the most widespread malaria vectors in the country. We investigated the gene expression profiles of insecticide resistant An. arabiensis populations from Migori and Siaya counties in Western Kenya using RNA-Sequencing. Centers for Disease Control and Prevention (CDC) bottle assays were conducted using deltamethrin (DELTA), alphacypermethrin (ACYP) and pirimiphos-methyl (PMM) to determine the resistance status in both sites.

Results

Mosquitoes from Migori had average mortalities of 91%, 92% and 58% while those from Siaya had 85%, 86%, and 30% when exposed to DELTA, ACYP and PMM, respectively. RNA-Seq analysis was done on pools of mosquitoes which survived exposure (‘resistant’), mosquitoes that were not exposed, and the insecticide-susceptible An. arabiensis Dongola strain. Gene expression profiles of resistant mosquitoes from both Migori and Siaya showed an overexpression mainly of salivary gland proteins belonging to both the short and long form D7 genes, and cuticular proteins (including CPR9, CPR10, CPR15, CPR16). Additionally, the overexpression of detoxification genes including cytochrome P450s (CYP9M1, CYP325H1, CYP4C27, CYP9L1 and CYP307A1), 2 carboxylesterases and a glutathione-S-transferase (GSTE4) were also shared between DELTA, ACYP, and PMM survivors, pointing to potential contribution to cross resistance to both pyrethroid and organophosphate insecticides.

Conclusion

This study provides novel insights into the molecular basis of insecticide resistance in An. arabiensis in Western Kenya and suggests that salivary gland proteins and cuticular proteins are associated with resistance to multiple classes of insecticides.

Peer Review reports

Background

The main malaria vector control methods in Kenya include the use of insecticide treated nets (ITNs) containing pyrethroids and indoor residual spraying (IRS) using organophosphates and neonicotinoids [1, 2]. The continued use of these insecticide-based interventions has led to increased resistance among malaria vectors in Kenya, where resistance to all four traditional classes of public health insecticides (pyrethroids, organophosphates, organochlorines and carbamates) have been reported [3]. Recently, resistance to neonicotinoids – a new class of insecticide used in IRS in many countries in sub-Saharan Africa (SSA) – was reported with clothianidin in Central Africa, raising an alarm and highlighting the urgent need for close and timely monitoring of vector susceptibility [4].

Anopheles arabiensis is one of the principal vectors of malaria in Kenya, as well as An. gambiae s.s and An. funestus [5, 6] and recently, An. stephensi has been detected in Northern Kenya [7]. Behavioral plasticity of An. arabiensis has been shown to compromise the protective effects of ITNs, as they bite both indoors and outdoors throughout the night [8]. Several studies conducted in western Kenya have reported insecticide resistance in An. arabiensis as well as other malaria vector species, causing great concern to the National Malaria Control Program (NMCP) [6, 9,10,11,12].

Metabolic resistance and target site mutations are the most widely studied and reported mechanisms of resistance. Metabolic resistance involves large enzyme families including cytochrome P450s (CYP450s), carboxylesterases (COEs) and glutathione-s-transferase (GSTs) which are known to confer resistance in malaria vectors. These enzymes exist naturally in mosquitoes and their amplification or overexpression leads to heightened detoxification of insecticides making the mosquitoes resistant [13]. Target site resistance arises from point mutations that alter the insecticide binding sites or transportation channels inside the mosquito [13]. Other modes behind insecticide resistance include: cuticular modifications such as thickening of the cuticles or change in cuticle composition, which prevent or slow insecticide penetration [14]; and behavioral resistance which results in mosquitoes avoiding surfaces treated with insecticides or in mosquitoes that are not affected by spatial repellents [15]. Investigating gene expression profiles of insecticide resistant malaria vectors is important in understanding the underlying mechanisms and in identification of markers that can be used for monitoring purposes. So far, several studies have identified marker genes associated with both pyrethroid and organophosphate resistance in Anopheles mosquitoes due to their high expression levels [16,17,18]. In an An. arabiensis population from Ethiopia, genes including CYP9K1, CYP9L1, GSTE4 as well as COEs were associated with metabolic resistance [17]. In addition to detoxification genes, a study conducted by Isaacs, et al. [18] showed that salivary gland proteins were implicated in insecticide resistance.

In this study, we sought to explore insecticide resistance mechanisms using RNA-Seq to identify markers associated with resistance to DELTA, ACYP and PMM in An. arabiensis populations from Migori and Siaya counties in western Kenya. These results will inform the NMCP decision making around IRM to ensure the efficacy of vector control interventions.

Results

Phenotypic insecticide resistance of An. arabiensis from Migori and Siaya

A total of 2404 and 2424 mosquitoes from Migori and Siaya counties, respectively, were exposed to either ACYP, DELTA or PMM. Out of those, 120 mosquitoes from each site (30 mosquitoes surviving exposure to each insecticide and 30 that were not exposed) were pooled in groups of 10 and sequenced alongside three pools of 10 susceptible An. arabiensis from the Dongola reference strain. In Migori, 100% mortality was observed at the diagnostic dose (1X) of ACYP, 2X DELTA and 2X PMM, while in Siaya, complete mortality was observed at 2X ACYP, 2X DELTA and 5X PMM. The Migori samples had an average mortality of 91% to 1X DELTA, 92% to 0.5X ACYP and 58% to 1X PMM. In Siaya, the average mortality was at 85% to 1X DELTA, 86% to 1X ACYP and 30% to 1.5X PMM. PCR tests conducted on the legs of the bioassayed mosquitoes confirmed that they were all An. arabiensis (Fig. 1).

Fig. 1
figure 1

Determination of phenotypic insecticide resistance profiles of An. arabiensis from Siaya and Migori counties using CDC bottle bioassays. The average mortalities of mosquitoes after 30 min of insecticide exposure are shown as percentages on the y-axis with 95% confidence intervals

Full size image

Descriptive summary statistics of RNA-Seq data

Whole transcriptomic analysis was done on mosquitoes resistant to DELTA, ACYP, and PMM, the unexposed mosquitoes from both Migori and Siaya as well as the susceptible Dongola mosquito strain. The total number of paired-end reads generated for all the samples were 2.18 billion ranging from 26 – 107 million reads per library. Samples collected from Migori had a total of 1 billion reads ranging from 56 – 107 million. Siaya had 921 million reads ranging from 26 – 101 million while the susceptible Dongola An. arabiensis generated a total of 228 million reads ranging from 69 – 85 million. From all the samples, in average 98% of the reads were retained after filtering and removal of adapters and in average 52% of the reads were mapped to the An. arabiensis Dongola reference genome (Additional file 1). A typical low mapping rate of RNA-Seq against the reference genome is not surprising and can be attributed to two factors: inherent incomplete ribosomal depletion and the presence of large number of multi-mapped reads that are not reported here due to their meaningless in DEG analysis. The feature counts results showing the number of tags generated and the percentages assigned to exons in the sense orientation has been summarized in Additional file 2 and Additional file 3. The percentage of tags assigned to exons in the sense direction ranged from 43 to 64% (Additional file 2, Additional file 3).

Differential gene expression associated with alphacypermethrin resistance

Three pairwise comparisons were done to determine the genes which were differentially expressed for ACYP in both Migori and Siaya (Table 1). For Migori, a total of 1088 (795 up and 293 down regulated), 1101 (630 up and 471 down regulated), and 317 genes (160 up and 157 down regulated) were significantly differentially expressed (FDR ≤ 0.01, |FC|≥ 2) in Res-Sus (MA vs DO), Con-Sus (MU vs DO) and Res-Con (MA vs SU), respectively (Table 1, Additional file 4A). For Siaya, a total of 1225 (838 up and 387 down regulated), 1155 (815 up and 340 down regulated) and 63 genes (26 up and 37 down regulated) were significantly differentially expressed in the Res-Sus (SA vs DO), Con-Sus (SU vs DO) and Res-Con (SA vs SU), respectively (Table 1, Additional file 4B).

Table 1 Results showing summaries of differential gene expression patterns of An. arabiensis for alphacypermethrin, deltamethrin and pirimiphos-methyl. Genes that were significantly differentially expressed at p < 0.01 and fold change (FC) > 2 were considered as candidate genes of interest
Full size table

The list of differentially expressed genes (DEGs) that are shared between two or more comparisons (Res-Sus, Con-Sus, Res-Con), as described in Additional file 4, were extracted and mapped to their Fold Change expression and functional description (Additional file 5). Focusing on the significantly differentially expressed genes shared between Res-Sus and Res-Con, Migori had a total of 70 DEGs (66 up and 4 down regulated) (Additional file 4A). Among the top 5 genes with retrievable annotations included a serine protease 53-like, a microfibril-associated glyco 4-like, a synaptic vesicle glyco 2B-like, a senecionine N-oxygenase and an uncharacterized protein. This group also included 2 cuticular proteins (flexible cuticle 12-like and adult cuticle 1-like), a carboxylesterase-6 like, and 3 cytochrome P450s (CYP307A1, CYP4C36 and CYP6M2). Siaya had a total of 23 DEGs (17 up and 6 down regulated) (Additional file 4B). The top 5 up regulated genes with retrievable annotations included 3 cuticular proteins (1 cuticular protein CPLCG family and 2 cuticular protein CPLCG family), an ATP-binding cassette transporter and a basic proline-rich -like gene. No detoxification genes were detected in this group.

It is worth noting that in all Res-Sus ACYP comparisons from both sites, the DEGs associated with cuticular proteins (CPs) and salivary gland proteins (SGPs), were predominantly overexpressed at higher proportions compared to detoxification genes. In addition, CPs, SGs and GSTs had higher ratios of up regulated genes (Additional file 6). Among the annotated DEGs in Migori, 92% of cuticular (34/37), 97% salivary (28/29), 72% CYP450s (13/18), 86% GSTs (6/7), and 50% COEs (3/6) were up regulated in the Res-Sus comparisons. In Siaya, 90% of cuticular (37/41), 96% salivary (25/26), 65% CYP450s (13/20), 90% GSTs (9/10), and 50% COEs (2/4) were up regulated (Additional file 6).

Identification of core differentially expressed salivary gland proteins, cuticular proteins and the detoxification genes associated with ACYP resistance was performed by selecting genes commonly shared in the Res-Sus comparisons from both sites. Here, there was a total of 30 cuticular proteins (28 up and 2 down regulated) and 25 salivary gland protein genes (24 up and 1 down regulated). The up regulated cuticular protein genes included CPAP3-A1b, CPAP3-Az1c, CPR10, CPR15, CPR16 and CPR9 while salivary gland protein genes included D7L1 and D7L2 (Fig. 2A, Additional file 7). There was a total of 26 detoxification genes which were differentially expressed in the ACYP Res-Sus comparisons across both sites. These included the following 18 (10 cytochromeP450s, 6 GSTE and 2 COEs) up regulated genes: CYP307A1, CYP4C27, CYP4H15, CYP6M2, CYP6M3, CYP6P4, CYP6Z3, CYP9K1, CYP9L1, CYP9M1, GSTD12, GSTD7, GSTE4, GSTE7, GSTE8, GSTU1 & 2 COEs as well as the following 8 (5 cytochromeP450s, 1 GSTE and 2 COEs) down regulated genes: CYP302A1, CYP325D1, CYP325H1, CYP4C35, CYP9M2, GSTE1 and 2 COEs (Table 2).

Fig. 2
figure 2

Venn diagram of differentially expressed genes (log2FC > 1 and FDR < 0.01), showing the number genes that commonly differentially expressed in the resistant populations of the two sites for each insecticide. A Alpha-cypermethrin; B Deltamethrin; and; C Pirimiphos-methyl

Full size image
Table 2 Table showing differentially expressed genes in the resistant vs susceptible (Res-Sus) and unexposed vs susceptible pairwise comparisons. FDR =  < 0.05, FC =  > 2
Full size table

Differential gene expression associated with deltamethrin resistance

In Migori, a total of 906 genes (657 up and 249 down regulated) were significantly differentially expressed in mosquitoes resistant to deltamethrin compared to the susceptible strain (MD vs DO). In a pairwise comparison of deltamethrin resistant and the unexposed mosquitoes (MD vs MU), there was a total of 318 differentially expressed genes (117 up and 201 down regulated). A total of 1101 genes (630 up and 471 down regulated) were differentially expressed in the pairwise comparison of the unexposed and susceptible mosquitoes (MU vs DO) (Table 1, Additional file 4A). In Siaya, there was a total of 1187 significantly differentially expressed genes (856 up and 331 down regulated) in the deltamethrin resistant vs. susceptible pairwise comparison (SD vs DO). In the resistant and unexposed pairwise comparison (SD vs SU), there was a total of 60 differentially expressed genes (26 up and 34 down regulated). The pairwise comparison of unexposed and susceptible mosquitoes (SU vs DO) had a total of 1155 differentially expressed genes (815 up and 340 down regulated) (Table 1, Additional file 4B).

Migori had a total of 58 (40 up and 18 down regulated) significantly differentially expressed genes that were shared between Res-Sus and Res-Con (Additional file 4A). Some of the top up regulated genes with retrievable annotations included 3 cuticular proteins, a cytochrome P450 (CYP307A1), a carboxylesterase (AARA001215), probable chitinase and 2 serine proteases. In Siaya, there was a total of 24 differentially expressed genes (17 up and 7 down regulated) shared between Res-Sus and Res-Con (Additional file 4B). The up regulated genes included 4 cuticular proteins: cuticle 7-like, 2 adult cuticle 1-like and histidine-rich PFHRP-II (Additional file 5).

Consistent with the alphacypermethrin results, in all Res-Sus DELTA comparisons from both sites, most of DEGs associated with CPs and SGPs were overexpressed at higher proportions (Fig. 3). In Migori, 89% CPs (34/38), 96% SGPs (27/28), and 67% GSTs (2/3) DEGs were up regulated in the Res-Sus comparisons, respectively. In Siaya, 89% CPs (33/37), 96% SGPs (22/23), and 89% GSTs (8/9) were up regulated, respectively (Additional file 6).

Fig. 3
figure 3

Gene expression profiles of resistant An. arabiensis from (A) Migori and (B) Siaya exposed to deltamethrin, pirimiphos-methyl and alpha-cypermethrin in comparison to the susceptible An. arabiensis Dongola strain. The horizontal dotted line on the volcano plot denotes a P-value of 0.01 while the vertical dotted lines indicate twofold expression differences. On the x-axis of each plot, genes that are overexpressed in the population are > 0. The -log10FDR values greater than 40 are displayed as 40. Generally, differentially expressed genes associated with cuticular proteins (CPs) and salivary gland proteins (SGPs) were over-expressed at higher proportions in all the six comparisons. COE = carboxylesterases; CP = cuticular protein; CYP = cytochrome P450s monooxygenases; GST = glutathione S-transferases; SGP = salivary gland protein

Full size image

Identification of core genes associated with DELTA resistance was done by selecting genes commonly shared in the Res–Sus comparisons from both sites (Fig. 2B, Additional file 7). A total of 32 (30 up and 2 down regulated) cuticular proteins and 19 up regulated salivary gland proteins were associated with DELTA resistance at both sites. Some of the cuticular protein genes included: CPAP3-A1b, CPAP3-A1c, CPR10, CPR15, CPR16 and CPR9 while some of the salivary gland protein genes included D7L1 and D7L2. There was a total of 14 detoxification genes which were significantly differentially expressed in the DELTA Res-Sus comparisons. These included the following 8 up regulated genes: CYP307A1, CYP4C27, CYP9L1, CYP9M1, GSTE4 & 3 COEs as well as the following 6 down regulated genes: CYP325D1, CYP325H1, CYP4C35, CYP4J10, GSTE1 and a COE (Table 2).

Differential gene expression associated with pirimiphos-methyl resistance.

In Migori, a pairwise comparison between pirimiphos-methyl resistant and susceptible mosquitoes (MP vs DO) had a total of 1507 differentially expressed genes (1007 up and 500 down regulated). A pairwise comparison of pirimiphos-methyl resistant and unexposed mosquitoes (MP vs MU) had a total of 357 differentially expressed genes (193 up and 164 down regulated) while a pairwise comparison of unexposed to susceptible (MU vs DO) mosquitoes had a total of 1101 differentially expressed genes (630 up and 471 down regulated) (Table 1, Additional file 4A). In Siaya, a total of 1496 genes (987 up and 509 down regulated) were significantly differentially expressed in pirimiphos-methyl resistant as compared to susceptible mosquitoes (SP vs DO). A comparison between pirimiphos-methyl resistant and unexposed mosquitoes (SP vs SU) had a total of 68 differentially expressed genes (42 up and 26 down regulated) while a comparison of unexposed and susceptible mosquitoes (SU vs DO) had a total of 1155 differentially expressed genes (815 up and 340 down regulated) (Table 1, Additional file 4B).

Focusing on the significantly differentially expressed genes overlapping the Res-Sus and Res-Con comparisons, Migori had a total of 116 genes (94 up and 22 down regulated) (Additional file 4A). Among the top genes which were up regulated and had retrievable annotations, there were CYP450s (CYP6M2, CYP6P2, CYP6Z3), 2 cuticular proteins, a probable chitinase 10 and a synaptic vesicle glyco 2B-like. Siaya had a total of 41 genes (29 up and 12 down regulated) that overlapped the Res-Sus and Res-Con comparisons (Additional file 4B). The top 15 up regulated genes included a CYP450 (CYP9M1), a multidrug resistance-associated 1-like gene, a chitinase partial and a serine protease. Five DEGs were found to overlap all three pairwise comparisons, including CYP450s (CYP6M2 and CYP6Z3) and a carboxylesterase (AARA007309) (Additional file 5). The consistent overexpression of the CYP6M2 and CYP6Z3 in the primiphos-methyl resistant mosquitoes from both Siaya and Migori independent of which group they were compared to (Con or Sus), highlights their potential contribution to the detoxification of this insecticide.

Similar to the ACYP and DELTA results, the gene expression profiles of PMM resistant An. arabiensis from both Migori and Siaya showed a higher proportion of CPs and SGPs that were overexpressed in the Res-Sus comparisons (Fig. 3). In Migori, 88% (37/42), 97% (32/33), and 91% (10/11) of the DEGs associated with CPs, SGPs, and GSTs were up regulated in the Res-Sus comparisons, respectively. In Siaya, 74% CPs (32/43), 93% SGPs (27/29), and 90% GSTs (9/10) were up regulated, respectively, in the Res-Sus comparisons (Fig. 3, Additional file 6).

Identification of core genes associated with PMM resistance was done by selecting genes commonly shared in the Res–Sus comparisons from both sites (Fig. 2C, Additional file 7). A total of 32 (27 up and 5 down regulated) CP and 26 (25 up and 1 down regulated) SGP genes were associated with PMM resistance (Fig. 2C, Additional file 7). Some of the up regulated genes included CPAP3-A1b, CPR10, CPR15, CPR16 & CPR9 cuticular proteins as well as D7L1 and D7L2 salivary gland proteins. Among the detoxification genes present, 21 (11 CYP450s, 7 GSTs and 3 COEs) were up regulated while 7 (5 CYP450s, 1 GST and 1 COE) were down regulated. Some of the up regulated genes included CYP4C27, CYP6M2, CYP6M3, CYP9L1, CYP9M1, GSTE4 and 3 COEs. The down regulated genes included CYP306A1, CYP325D1, CYP325H1, CYP4C35, CYP9M2 GSTE1 and 1 COE (Table 2).

Shared genes associated with ACYP, DELTA and PMM insecticide resistance across collection sites

From both Siaya and Migori, comparing across ACYP, DELTA and PMM resistant and susceptible mosquitoes, 570 genes (446 up and 124 down regulated) were significantly differentially expressed (Fig. 4A). The following up regulated genes were among the top 20 which had retrievable annotations: 2 chitinases, 2 serine protease 53-like, a salivary gland protein, 2 thioester-containing partial and a neuronal acetylcholine receptor subunit.

Fig. 4
figure 4

Identification of DEGs associated with resistance to multiple insecticides. A Upset plot representing the intersection of DEGs between DELTA, ACYP, PMM resistant mosquitoes from Siaya and Migori when compared to the susceptible An. arabiensis Dongola strain (Res-Sus comparisons). The left horizontal bar plot (set size) reports the total number of DEGs in each comparison, the circles represent the set of comparisons associated with the intersection, while the vertical bar plot reports the number of unique and overlapped DEGs (intersection size) between the different combinations of the R-S comparisons. The highlighted bar plot represents core DEGs commonly shared across PMM, DELTA and ACYP, (red), DEGs specific to PMM (blue), DELTA (green) and ACYP (yellow). B Heatmap of the log2 fold change (log2FC) expression of all cuticular and salivary gland protein core DEGs that were differentially expressed in all the resistant vs susceptible pairwise comparisons from both sites. C Heatmap representing the log2FC expression of the top 10 detoxification genes shared between all insecticide resistant vs susceptible pairwise comparisons for both sites. The heatmaps are in a blue-red color gradient (red = over-expressed and blue = under-expressed). SA = Siaya alpha-cypermethrin, SD = Siaya deltamethrin, SP = Siaya pirimiphos-methyl, SU = Siaya unexposed, MA = Migori alpha-cypermethrin, MD = Migori deltamethrin, MP = Migori pirimiphos-methyl, MU = Migori unexposed, DO = Dongola (An. arabiensis susceptible strain)

Full size image

Among the 570 genes, there were 21 (20 up and 1 down regulated) cuticular proteins, 19 up regulated salivary gland proteins (Fig. 4B), 6 (3 up and 3 down regulated) CYP450s, 2 (1 up and 1 down regulated) GSTs and 2 up regulated COEs which were shared across sites and insecticides (Fig. 4C). Some genes of interest included a glutactin-like COE (AARA016468; FCs = 8.53, 8.68, 8.90, 10.13, 11.74 and 12.09), and cuticular protein genes such as a cuticle-like CPR16 (AARA002342; FCs = 6.41, 7.41, 8.73, 9.27, 10.46 and 11.82), an endocuticle structural glyco ABD-4 (AARA003903; FCs = 6.91, 8.73, 9.70, 11.19, 12.53 and 13.39) and an endocuticle structural glyco ABD-5-like (AARA016140; FCs = 6.08, 6.09, 6.55, 7.68, 8.28 and 15.14). Focusing on insecticide-specific genes, a total of 110 genes were associated only with PMM resistance while 11 genes were associated with DELTA and ACYP resistance (Fig. 4A). Genes such as CYP325K1 and AARA017332 whose ortholog in An. gambiae is COEJHE2E were specific to PMM resistance in addition to cuticular (AARA005785 and AARA007248) and salivary gland protein genes (AARA016215). The CYP4J10 gene was only specific to DELTA resistance. No gene with functional validation associated with insecticide resistance was present for ACYP.

Gene Ontology enrichment analysis

Gene Ontology (GO) enrichment analysis of the DEGs detected from each Res-Sus comparison (n = 6) was conducted using Goatools [19]. The list of enriched GO terms associated with the up and down regulated genes for each comparison is reported in Additional file 8. The overlap of the significantly enriched GO Biological Process (BP) terms across the 6 comparisons is depicted Fig. 5A, while the enriched GO terms of the overexpressed genes that overlapped all the six comparisons are shown in Fig. 5B. A total of nine GO terms were significantly enriched in all the comparisons, including “carbohydrate metabolic process” (GO:0005975), “energy derivation by oxidation of organic compounds” (GO:0015980), “generation of precursor metabolites and energy” (GO:0006091), “metabolic process” (GO:0008152), “purine-containing compound biosynthetic process” (GO:0072522), “regulation of proteolysis” (GO:0030162), “ribose phosphate metabolic process” (GO:0019693), “small molecule metabolic process” (GO:0044281), and “sulfur compound metabolic process” (GO:0006790). Not surprisingly, the number of enriched GO terms was positively correlated with the number of DEGs for the six comparisons, suggesting that the change in metabolic pathway level in mosquitoes is strongly associated with the changes in the expression of individual genes.

Fig. 5
figure 5

Gene ontology enrichment analysis of the overexpressed genes. The upset plot (A) depicts the number of unique and shared GO terms from the functional enrichment analysis. The horizontal bars (blue) indicate the number of enriched GO terms in each comparison (Res vs Sus), while the vertical bars (black) represent number of overlapping GO terms in the sets, indicated with black dot under each bar. The heatmap (B) represents the -log10(FDR) of the 9 enriched GO terms that overlap all the comparisons. GO enrichment analysis results are shown for only the overexpressed genes detected in the resistant group vs susceptible (Res vs Sus). The complete report of GO enrichment analysis is shown in Additional file 8

Full size image

Genetic variation analyses on Voltage-Gated Sodium Channel (Vgsc) and acetylcholinesterase (Ace-1) genes

RNA-Seq reads from pools of 10 mosquitoes in each representing the ACYP, DELTA and PMM experimental replicates were used for Single Nucleotide Polymorphism (SNP) analysis to approximate the allele frequencies at target site mutations in the kdr, ACE-1 and GSTE2 loci. With a focus on the non-synonymous SNPs on the Vgsc gene described by Clarkson, et al. [20], L995S was detected at a frequency of 33% in Migori and 17% in Siaya ACYP resistant mosquitoes and L995F was detected at 50% in both Siaya ACYP and DELTA resistant mosquitoes. No SNPs were detected in the ACE-1 and GSTE2 genes (Additional file 9). These results suggest that insecticide resistance in the study populations is mainly of a metabolic nature.

Validation of gene expression levels estimated by RNA-Seq using Quantitative Real-Time PCR (qRT-PCR)

To complement the RNA-Seq results, we conducted qRT-PCR to validate some of the differentially expressed genes including chitinase, salivary gland protein, NADH dehydrogenase and the D7 genes alongside two housekeeping genes, RPS7 and Actin5c. Most of the qRT-PCR data including the D7 and NADH endorsed the directionality of expression as estimated by RNA-Seq Fig. 6. However, qRT-PCR of the chitinase gene was not congruent with the RNA-Seq results, likely due to the low transcriptional signal of this gene as reflected in the RNA-Seq (Additional file 10).

Fig. 6
figure 6

qPCR validation of RNA-Seq gene expression levels in Siaya alpha-cypermethrin (SA), deltamethrin (SD), pirimiphos-methyl (SP) and unexposed (SU) mosquitoes, versus the susceptible Dongola strain (DO)

Full size image

Discussion

High throughput sequencing platforms have enabled genomic and transcriptomic-level studies of the genetic profiles of insecticide resistant malaria vectors. The overarching goal of this study was to investigate the transcriptomic profiles of DELTA, PMM and ACYP resistant An. arabiensis from western Kenya. Results from this study showed that An. arabiensis from the counties of Siaya and Migori had varied levels of phenotypic resistance to DELTA, PMM and ACYP. This resistance was associated with elevated expression of salivary gland and cuticular protein genes in addition to detoxification genes including CYP450s, COEs and GSTE.

The phenotypic resistance in the An. arabiensis populations from Siaya and Migori varied across the different insecticides. The higher frequency of PMM resistance compared to the two pyrethroids was unexpected. The use of insecticides in agriculture could contribute to selective pressure since some of the insecticides used in agriculture contain the same active ingredients as those used for vector control [21]. Whereas pyrethroid resistance had been previously reported in An. arabiensis from western Kenya at varying intensities, the vectors remained susceptible to organophosphate insecticides [12, 22]. The emergence of resistance to organophosphates may be associated with run off from insecticides used in agriculture into larval habitats [3, 23]. In addition, the use of organophosphates has now been introduced for IRS in western Kenya and continued usage will likely result in increasing resistance over time.

In both sites, there was a high level of overexpression of cuticular and salivary gland proteins, as well as detoxification genes, some of which have been previously reported [17, 18]. Salivary gland proteins are known to play important roles during blood feeding in mosquitoes such as releasing anti-coagulants and are thus important to transmission of malaria parasites [24]. In other cases, the SGPs aid in transmission of viruses by mosquitoes [25]. Here, we report overexpression of twelve salivary gland proteins present in insecticide resistant mosquitoes which included the D7 long (D7L1 and D7L2) and short form which have previously been associated with insecticide resistance [18]. Interestingly, eleven of these genes have been found to be present in An. arabiensis populations from Ethiopia, some of which were found to be overexpressed in the pyrethroid and organophosphate resistant group [17] and may point to their contribution to insecticide resistance in An. arabiensis across the East African region.

The ortholog of the D7 short form gene (AARA016237) in An. gambiae, D7r4, was found to be overexpressed in carbamate resistant mosquitoes [18]. The presence of D7r4 gene in both An. arabiensis and An. gambiae resistant to pyrethroid, organophosphate and carbamate insecticides could suggest the possibility of a cross resistance mechanism [17, 24]. A recent study by Freitas and Nery [24] identified several SGPs which could be considered as potential targets for novel vector control strategies. Further investigation will be necessary to ascertain the role of SGP genes associated with insecticide resistance.

Cuticular protein genes including CPAP3-A1b, CPR9, CPR10, CPR15 and CPR16 were significantly differentially expressed in mosquitoes showing resistance to all the three insecticides from both sites. The insect cuticle is mainly composed of chitin and cuticular proteins such as those of the CPR and CPAP families which protects the insect from harsh weather conditions. Any alterations to the cuticle composition could lead to changes such as thickening which may impede insecticide penetration and render the mosquitoes resistant. Some cuticular proteins such as the CPAP3, CPR and CPLCG families observed in the resistant mosquitoes in this study have previously been associated with resistance to different insecticides in Anopheles mosquitoes [26,27,28,29]. A study done by Zoh, et al. [29] showed that the CPAP3-A1b gene was linked to clothianidin (neonicotinoid) resistance in An. gambiae s.s from Côte d'Ivoire (Tiassale). This gene was also overexpressed in pyrethroid and organophosphate resistant An. arabiensis samples analyzed in this study and also in pyrethroid resistant An. arabiensis from Tanzania [30]. Although this gene can be considered as a potential marker for insecticide resistance, a further functional verification step is needed to confirm the role of the CPAP3-A1b gene in both An. gambiae and arabiensis populations resistant to different classes of insecticides. Additional studies by Yahouédo, et al. [28] Zhou, et al. [31] have provided evidence on the presence of genes belonging to the CPR and CPLCG families in insecticide resistant Anopheles mosquitoes.

Some members of the CYP450s, such as the CYP4G family, are involved in cuticle development by catalyzing the production of cuticular hydrocarbons, the most abundant lipid species in the epicuticle [32]. Studies conducted by Balabanidou, et al. [32] and Yahouédo, et al. [28] demonstrated that insecticide resistant mosquitoes had thicker cuticles compared to susceptible mosquitoes due to the overexpression of the cuticular protein genes that led to enriched deposition of CHCs. It will be important to further investigate the role of CYP540 genes overexpressed in resistant samples to determine their contribution to cuticle development.

In both Migori and Siaya, several detoxification genes were overexpressed in mosquitoes showing resistance to all insecticides. These detoxification genes included CYP9M1, CYP4C27, CYP9L1, GSTE4 and two COEs (AARA004790 and AARA016468). Previous studies have reported the presence of these genes, with the exception of CYP9M1, in mosquitoes resistant to insecticides including pyrethroids and organophosphates [17, 33,34,35]. Their presence in our study further supports their association with resistance and their potential for use as molecular insecticide resistance markers. Genes such as GSTE4, CYP9L1 and the two COEs (AARA004790 and AARA016468) have previously been identified in resistant An. arabiensis and might be considered as species-specific markers for resistance to pyrethroids and organophosphates [17]. GSTs belonging to the delta and epsilon classes have previously been demonstrated to metabolize insecticides including pyrethroids, organophosphates and organochlorines [33, 36,37,38,39,40,41]. Here, GSTE4 was found to be over-expressed across all samples from both sites resistant to pyrethroids and organophosphate. Although the role of GSTE4 in metabolizing insecticides is not fully understood, a previous study [42] suggests that this enzyme could be involved in sequestration of insecticides. Functional validation will be required to understand the role of GSTEs including GSTE4 in insecticide resistance in An. arabiensis.

Besides genes that were associated with multiple insecticide resistance, some were specific to the different insecticides. Focusing on genes with functional validation associated with insecticide resistance, CYP4J10 was found to be specific to only DELTA. This gene has previously been associated with resistance to permethrin which is also a pyrethroid [43]. Genes such as CYP325K1 which was specific to PMM resistance has been associated with pyrethroid resistance in An albimanus and Aedes aegypti [16, 44]. The gene AARA017332 whose ortholog in An. gambiae is COEJHE2E has been linked with permethrin resistance in An. arabiensis and here, PMM resistance suggesting that it might not be an insecticide specific gene [33]. In addition to the detoxification genes, there were some cuticular (AARA005785 and AARA007248) and salivary gland protein genes (AARA016215) which were specific to PMM but have not previously been associated with resistance.

Taken together these results describe the suite of marker genes involved in An. arabiensis resistance to key pyrethroid and organophosphate insecticides used in malaria vector control. Once validated, they can be considered as candidate molecular markers that National Malaria Programs can use to routinely monitor for the emergence of insecticide resistance. In addition, these results add to the growing body of knowledge on the molecular basis of insecticide resistance within key vector populations [45].

Conclusion

In this study, An. arabiensis populations from Siaya and Migori counties were found to be resistant to ACYP, DELTA and PMM insecticides at different levels. Furthermore, transcriptomic analysis revealed novel resistance genes in Kenyan An. arabiensis in addition to those that have previously been described in other countries. Gene expression profiles showed an overexpression of the salivary gland proteins belonging to both the short and long form D7 genes, and cuticular proteins (including CPR9, CPR10, CPR15, CPR16) that were shared between pyrethroid and organophosphate resistant mosquitoes. In addition, detoxification genes including CYP450s (CYP9M1, CYP325H1, CYP4C27, CYP9L1 and CYP307A1), 2 COEs and a GSTE (GSTE4) were also over-expressed. A functional validation of these markers is needed to confirm the role these markers in conferring insecticide resistance.

Methods

Larval sampling and rearing

Sampling of Anopheles larvae in Migori (1.0707° S, 34.4753° E) and Siaya (0.0626° N, 34.2878° E) Fig. 7 was conducted between August and October 2019. Larvae were sampled using dippers and placed in collection tins with labels indicating the collection date and site. The collected larvae were transported to and reared using standard methods described by Benedict [46] at the insectary at the Centre for Global Health Research, Kenya Medical Research Institute, Kisumu (KEMRI- CGHR). Resulting adult mosquitoes were sustained on a 10% sugar solution soaked in cotton balls and held for 3–5 days before exposure to insecticides for susceptibility testing. All the collected larvae from both sites were reared under temperatures of (30 ± 2 °C), while adults were reared at a temperature of 27 ± 2 °C, relative humidity of 80 ± 10%, and photoperiod of 12:12 light: dark cycle.

Fig. 7
figure 7

Map of Kenya (right) showing Migori and Siaya counties in expanded view where An. arabiensis were collected

Full size image

Insecticide resistance phenotyping

CDC bottle bioassays were conducted on the 3–5 days old F0 adult female mosquitoes following the U.S. Centers for Disease Control and Prevention (CDC) guidelines for evaluating insecticide resistance [47]. Mosquitoes from both Migori and Siaya were exposed to ACYP, DELTA and PMM insecticides at varying concentrations ranging from 0.5X to 5X the diagnostic doses. Technical grade insecticides were diluted using acetone in 50 ml falcon tubes to obtain stock diagnostic concentrations of 12.5 μg/bottle for ACYP, 12.5 μg/bottle for DELTA and 20 μg/bottle for PMM (Additional file 11).

Each experiment comprised of 1 control bottle treated with 1 ml of acetone and 4 test bottles each treated with 1 ml of the respective insecticide solution and left overnight to dry except for the PMM bottles which were dried for only 3 h prior to bioassay. Approximately 15–20 mosquitoes were introduced in each bottle using a mouth aspirator and exposed for the diagnostic time of 30 min. Mosquitoes were considered resistant if they were capable of standing and flying in a coordinated manner and susceptible if they died or were moribund. The unexposed mosquitoes were the field mosquitoes used in the control bottle that contained no insecticide. Resistant and unexposed mosquitoes were knocked down on ice to immobilize them then 2 -3 legs were cut and placed in a sterile 1.5 ml Eppendorf tube for species identification. The remaining carcass was preserved individually in a 1.5 ml Eppendorf tube containing RNA later and stored at 4 °C until shipment for sequencing at the U.S. Centers for Disease Control and Prevention (CDC) in Atlanta, USA.

Molecular species identification

Genomic DNA was isolated from the legs of the mosquitoes using the ethanol precipitation procedure [48]. 1μl of DNA from each sample was utilized as a template for the PCR process, along with known An. arabiensis DNA as a positive control [49]. The reactions were carried out using a BIORAD T100 thermal cycler under the following conditions: 95°C for 5 min, 30 cycles of: 95°C for 30 s, 56°C for 30 s, and 72°C for 30 s, and a final extension at 72°C for 5 min. Primers used were specific to An. arabiensis—5'AAG TGT CCT TCT CCA TCC TA 3', An. gambiae -5' CTG GTT TGG TCG GCA CGT TT 3' and a universal primer—5' GTG TGC CCC TTC CTC GAT GT 3'. A 2% agarose gel stained with ethidium bromide was used to visualize the 315bp amplicons for An. arabiensis.

RNA extraction, RNA-Seq library preparation and sequencing

RNA-Seq analysis included 3 groups of mosquitoes: susceptible An. arabiensis from the Dongola reference strain, resistant, and unexposed An. arabiensis from Migori and Siaya counties. The samples were labelled as DO (susceptible An. arabiensis Dongola strain), MA (Migori ACYP resistant), MP (Migori PMM resistant), MD (Migori DELTA resistant), MU (Migori unexposed), SA (Siaya ACYP resistant), SD (Siaya DELTA resistant), SP (Siaya PMM resistant), SU (Siaya unexposed). Total RNA isolation was carried out for 3 replicates of each group, each of which contained a pool of 10 mosquitoes. This was done using the Arcturus® PicoPure® RNA isolation kit (Life Technologies, USA) and quantification was carried out using the 4200 TapeStation System (Agilent Technologies, Palo Alto, CA, USA), according to the manufacturers’ protocols. RNA was treated using Baseline-ZERO™ DNase (Epicentre, Illumina) and removal of ribosomal RNA was done using the Ribo-Zero rRNA removal kit (Human/Mouse/Rat) (Epicentre, Illumina). The ScriptSeq v2 RNA-Seq Library Preparation Kit (Epicentre, Illumina) was used to prepare the individual libraries according to the manufacturer's instructions. Equimolar amounts of each library were pooled and sequenced (2 × 125 bp paired—end) on an Illumina HiSeq 2500 sequencer, using v2 chemistry. Sequencing was carried out at the Biotechnology Core Facility at CDC in Atlanta, USA.

RNA-Seq data analysis

Quality control filtering and mapping

For each sample, FastQC v0.11.5 was used to assess the quality of de-multiplexed paired end reads generated after sequencing [50]. Sequencing reads were trimmed and filtered using fastp v0.20.1 [51] to remove adapter and low-quality reads. Parameters used in fastp included a minimum base quality value of 20, required minimum length of 50, trimming of polyG tail and trimming of polyX in the 3’ ends to remove the low-complexity consecutive bases. Subsequently, the raw sequencing reads from similar sequencing lanes of R1 (forward) and R2 (reverse) were concatenated to increase the read depth to be used in subsequent downstream analysis. Trimmed reads were then aligned to the An. arabiensis Dongola reference genome assembly (genome assembly version: AaraD1, GeneBank assembly identifier = GCA_000349185.1) downloaded from VectorBase (release 48) using ‘subjunc’ v2.0.1, part of the subread aligner package with default parameters [52].

The resulting alignment was filtered to remove reads of low mapping quality (q < 10) and sorted using Samtools v1.10 [53]. Tag counting was done using ‘featureCounts’, part of the subread aligner package. Tags were defined as either a read pair or single, unpaired read. Aligned reads with at least 1 bp overlap in coding sequence (CDS) features in the sense orientation were counted, and the tabulated tag counts used as input for differential expression analysis with edgeR [54].

Differential gene expression analysis

Differential gene expression analysis was conducted for all groups. Prior to the analysis, genes that had a total tag count of more than 30 across all libraries were considered while the rest were filtered out to remove the lowly expressed genes. Variation in RNA-Seq data was modelled using a negative binomial distribution and a generalized linear model. The estimated log2 fold change (FC) for each gene was tested in edgeR using a Likelihood-Ratios (LR) test. P-values associated with log2FC were adjusted for multiple testing using the False Discovery Rate (FDR) approach and significance was defined as genes differentially expressed with an FDR-adjusted P-value < 0.01 [55].

In each site and for each insecticide, three pairwise comparisons were made between resistant vs susceptible (Res-Sus), unexposed vs susceptible (Con-Sus) and resistant vs unexposed (Res-Con). The unexposed mosquitoes are hereafter referred to as the control (Con) group. Res-Sus comparisons were to account for genes differentially expressed because of constitutive differences related to resistant phenotypes. Res-Con comparisons were to identify the genes which could have been induced because of exposure to the insecticides. Con-Sus comparisons were to explain overexpressed genes which are always present within a population since they are transcribed constitutively. Genes which were significantly differentially expressed at false discovery rate (FDR) < 0.01 and fold change (FC) > 2 were considered potential candidate markers for resistance. Of particular interest were the genes which were consistently and significantly differentially expressed across the different pairwise comparisons as these could potentially have been involved in cross resistance.

Gene ontology annotation and functional enrichment analysis

The gene ontology and functional annotation of the AaraD1.11 gene set (Genome version = GCA_000349185.1) were performed locally using blast2GO command line v1.4.4 [56] as follows. First, a local BLASTp (v2.9) search of the predicted protein coding sequences was conducted against the Arthropoda (taxid = 6665) category of the nr protein NCBI database with maximum e-value cut-off 10–3. Second, the protein sequences were searched against the InterPro database [57], using InterProScan v5 [58]. The Blastp and InterProScan outputs were simultaneously provided to Blast2GO command line as input, which map the RefSeq and InteProScan identifiers to the GO database as curated and updated in the Blast2GO database (August, 2022).

Subsequently, gene ontology enrichment (GOE) analysis was carried out on differentially expressed gene sets using GOATools [19]. The resulting annotated genes and their associated GO terms were used as reference for this GOE analysis. Fisher’s exact test was used to identify gene ontologies significantly enriched from the over- and under-expressed gene sets relative to the rest of the genome. An FDR adjusted P-value < 0.05 was used to determine the significantly enriched GO terms associated with the list of DEGs, while the redundant GO terms were eliminated using REVIGO (available at http://revigo.irb.hr).

Genetic variation analyses on Vgsc and Ace-1 genes

The RNA-Seq reads of the insecticide resistant, unexposed, and susceptible groups were mined for non-synonymous SNPs in the Vgsc gene (AARA017729) including V402L, D466H, M490I, G531V, Q697P, T791M, L995S, L995F, V1507I, I1527T, N1570Y, E1597G, K1603T, A1746S, V1853I, I1868T, P1874S, P1874L, A1934V, and I1940T, previously studied by Clarkson, et al. [20] and Messenger, et al. [17]. Additionally, we examined the I114T, L119F, L119V variants in GSTE2 gene (AARA008732) explored by Simma, et al. [59] and Messenger, et al. [17]. Furthermore, we investigated the G119S in ACE1 gene (AARA001814) previously studied by Weill, et al. [60], Simma, et al. [59] and Messenger, et al. [17]. The sorted BAM files used for differential gene expression analysis (see above) were used as input files for the SNP analysis using the SAMtools package v1.19 and BCFtools v.1.9 [53]. VCFtools v0.1.17 [61] was then used to generate the allele frequencies of the SNPs, and the variants of interest were extracted.

Gene validation using qRT-PCR

A quantitative real-time reverse transcription PCR (qRT-PCR) was carried out to validate a set of candidate marker genes that were significantly differentially expressed alongside two An. gambiae housekeeping genes: Actin (AGAP000651) and rPS7 (AGAP010592) that also worked with An. arabiensis samples. The housekeeping genes were used for normalization of the experiment. RNA was extracted from different mosquitoes of the F0 population used for RNA-Seq using the Arcturus™ PicoPure™ RNA Isolation Kit (Applied Biosystems, USA) following the manufacturer’s instructions. The extracted mosquitoes were from three replicates ACYP, DELTA and PMM resistant, unexposed and susceptible Dongola strain mosquitoes. cDNA was synthesized from 1ul of the extracted RNA using the High-Capacity cDNA Reverse Transcription kit (Applied Biosystems, USA) and oligo-d(T)23 (New England Biolabs, USA) in accordance with the manufacturer's instructions. A serial dilution of cDNA was used to create standard curves of Ct values for each gene, effectively nullifying any inaccuracies or deviations resulting from sample concentration. The qRT-PCR were performed using the PowerUp SYBR Green Master Mix (Applied Biosystems, USA) on a QuantStudio 3 Real-Time PCR system (Applied Biosystems, USA). Information about the primers including their sequences and their respective gene IDs is as documented in Additional file 12. The thermal cycling conditions were set at 50 ˚C for 2 min, 95 ˚C for 10 min, and thereafter, 40 cycles of: 15 s at 95 ˚C, 1 min at 60 ˚C, followed by 15 s at 95 ˚C, 1 min at 60 ˚C and a final step of 15 s at 95 ˚C. The 2 − ΔΔCT approach was used to calculate the relative expression level and Fold Change (FC) of each target gene from resistant field samples compared to the susceptible lab strain.

Availability of data and materials

Sequence data generated by this study is available at Sequence Read Archive (SRA) under the Bio Project accession number PRJNA986474.

Custom scripts used for all the analysis are available from the authors on request.

Abbreviations

ACE1:

Acetylcholinesterase

ACYP:

Alphacypermethrin

adjP:

P-value adjusted for multiple testing

An:

Anopheles

BP:

Biological process

CDC:

Centers for Disease Control and Prevention

COE:

Carboxylesterases

Con:

Control

CP:

Cuticular protein

CYP450s:

Cytochrome P450s

DEGs:

Differentially expressed genes

DELTA:

Deltamethrin

DNMP:

Division for National Malaria Program

FC:

Fold change

GO:

Gene ontology

GSTE:

Glutathione-S-transferase

IRM:

Insecticide resistance management

IRS:

Indoor Residual Spraying

Kdr:

Knock down resistance

KEMRI:

Kenya Medical Research Institute

LLINs:

Long lasting insecticidal nets

NMCP:

National Malaria Control Program

PMM:

Pirimiphos-methyl

qRT-PCR:

Real –time quantitative polymerase chain reaction

Res:

Resistant

SG:

Salivary gland

SNP:

Single nucleotide polymorphism

SSA:

Sub-Saharan Africa

Sus:

Suceptible

References

  1. DNMP K, ICF: Kenya Malaria Indicator Survey 2020. In. Nairobi and Rockville; 2021.

  2. The PMI VectorLink Project Kenya: Kenya end of spray report 2021: spray campaign: February 15 – April 13, 2021. Rockville; 2021.

  3. Ondeto BM, Nyundo C, Kamau L, Muriu SM, Mwangangi JM, Njagi K, Mathenge EM, Ochanda H, Mbogo CM. Current status of insecticide resistance among malaria vectors in Kenya. Parasit Vectors. 2017;10(1):1–13.

    Article  Google Scholar 

  4. Makoni M. Malaria fighters’ latest chemical weapon may not last long. Science. 2020;369(6508):1153.

    Article  CAS  PubMed  Google Scholar 

  5. Zhong D, Hemming-Schroeder E, Wang X, Kibret S, Zhou G, Atieli H, Lee M-C, Afrane YA, Githeko AK, Yan G. Extensive new Anopheles cryptic species involved in human malaria transmission in western Kenya. Sci Rep. 2020;10(1):1–13.

    Article  Google Scholar 

  6. Wanjala CL, Kweka EJ. Malaria vectors insecticides resistance in different agroecosystems in Western Kenya. Front Public Health. 2018;6:55.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Ochomo EO, Milanoi S, Abongo B, Onyango B, Muchoki M, Omoke D, Olanga E, Njoroge L, Juma E, Otieno JD. Detection of Anopheles stephensi Mosquitoes by Molecular Surveillance, Kenya. Emerg Infect Dis. 2023;29(12):2498–508.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Perugini E, Guelbeogo WM, Calzetta M, Manzi S, Virgillito C, Caputo B, Pichler V, Ranson H, Sagnon NF, Della Torre A. Behavioural plasticity of Anopheles coluzzii and Anopheles arabiensis undermines LLIN community protective effect in a Sudanese-savannah village in Burkina Faso. Parasit Vectors. 2020;13(1):1–10.

    Article  Google Scholar 

  9. Kawada H, Futami K, Komagata O, Kasai S, Tomita T, Sonye G, Mwatele C, Njenga SM, Mwandawiro C, Minakawa N. Distribution of a knockdown resistance mutation (L1014S) in Anopheles gambiae ss and Anopheles arabiensis in Western and Southern Kenya. PLoS ONE. 2011;6(9):e24323.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Ochomo E, Bayoh NM, Kamau L, Atieli F, Vulule J, Ouma C, Ombok M, Njagi K, Soti D, Mathenge E. Pyrethroid susceptibility of malaria vectors in four Districts of western Kenya. Parasit Vectors. 2014;7(1):1–9.

    Article  Google Scholar 

  11. Omondi S, Mukabana WR, Ochomo E, Muchoki M, Kemei B, Mbogo C, Bayoh N. Quantifying the intensity of permethrin insecticide resistance in Anopheles mosquitoes in western Kenya. Parasit Vectors. 2017;10(1):548.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Owuor KO, Machani MG, Mukabana WR, Munga SO, Yan G, Ochomo E, Afrane YA. Insecticide resistance status of indoor and outdoor resting malaria vectors in a highland and lowland site in Western Kenya. PLoS ONE. 2021;16(3):e0240771.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

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

    Article  CAS  PubMed  Google Scholar 

  14. Balabanidou V, Grigoraki L, Vontas J. Insect cuticle: a critical determinant of insecticide resistance. Curr Opin Insect Sci. 2018;27:68–74.

    Article  PubMed  Google Scholar 

  15. Chareonviriyaphap T, Bangs MJ, Suwonkerd W, Kongmee M, Corbel V, Ngoen-Klan R. Review of insecticide resistance and behavioral avoidance of vectors of human diseases in Thailand. Parasit Vectors. 2013;6(1):1–28.

    Article  Google Scholar 

  16. Mackenzie-Impoinvil L, Weedall GD, Lol JC, Pinto J, Vizcaino L, Dzuris N, Riveron J, Padilla N, Wondji C, Lenhart A. Contrasting patterns of gene expression indicate differing pyrethroid resistance mechanisms across the range of the New World malaria vector Anopheles albimanus. PLoS ONE. 2019;14(1):e0210586.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Messenger LA, Impoinvil LM, Derilus D, Yewhalaw D, Irish S, Lenhart A. A whole transcriptomic approach provides novel insights into the molecular basis of organophosphate and pyrethroid resistance in Anopheles arabiensis from Ethiopia. Insect Biochem Mol Biol. 2021;139:103655.

    Article  CAS  PubMed  Google Scholar 

  18. Isaacs AT, Mawejje HD, Tomlinson S, Rigden DJ, Donnelly MJ. Genome-wide transcriptional analyses in Anopheles mosquitoes reveal an unexpected association between salivary gland gene expression and insecticide resistance. BMC Genomics. 2018;19(1):1–12.

    Article  Google Scholar 

  19. Klopfenstein DV, Zhang L, Pedersen BS, Ramírez F, Warwick Vesztrocy A, Naldi A, Mungall CJ, Yunes JM, Botvinnik O, Weigel M, et al. GOATOOLS: A Python library for Gene Ontology analyses. Sci Rep. 2018;8(1):10872.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Clarkson CS, Miles A, Harding NJ, O’Reilly AO, Weetman D, Kwiatkowski D, Donnelly MJ. Consortium AgG: The genetic architecture of target-site resistance to pyrethroid insecticides in the African malaria vectors Anopheles gambiae and Anopheles coluzzii. Mol Ecol. 2021;30(21):5303–17.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Nkya TE, Akhouayri I, Kisinza W, David J-P. Impact of environment on mosquito response to pyrethroid insecticides: facts, evidences and prospects. Insect Biochem Mol Biol. 2013;43(4):407–16.

    Article  CAS  PubMed  Google Scholar 

  22. Orondo PW, Nyanjom SG, Atieli H, Githure J, Ondeto BM, Ochwedo KO, Omondi CJ, Kazura JW, Lee M-C, Zhou G. Insecticide resistance status of Anopheles arabiensis in irrigated and non-irrigated areas in western Kenya. Parasit Vectors. 2021;14(1):1–10.

    Article  Google Scholar 

  23. Munywoki DN, Kokwaro ED, Mwangangi JM, Muturi EJ, Mbogo CM. Insecticide resistance status in Anopheles gambiae (sl) in coastal Kenya. Parasit Vectors. 2021;14(1):1–10.

    Article  Google Scholar 

  24. Freitas L, Nery MF. Positive selection in multiple salivary gland proteins of Anophelinae reveals potential targets for vector control. Infect Genet Evol. 2022;100:105271.

    Article  CAS  PubMed  Google Scholar 

  25. Sun P, Nie K, Zhu Y, Liu Y, Wu P, Liu Z, Du S, Fan H, Chen C-H, Zhang R. A mosquito salivary protein promotes flavivirus transmission by activation of autophagy. Nat Commun. 2020;11(1):1–15.

    Google Scholar 

  26. Nkya TE, Akhouayri I, Poupardin R, Batengana B, Mosha F, Magesa S, Kisinza W, David J-P. Insecticide resistance mechanisms associated with different environments in the malaria vector Anopheles gambiae: a case study in Tanzania. Malar J. 2014;13(1):1–15.

    Article  Google Scholar 

  27. Nkya TE, Poupardin R, Laporte F, Akhouayri I, Mosha F, Magesa S, Kisinza W, David J-P. 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(1):1–12.

    Google Scholar 

  28. Yahouédo GA, Chandre F, Rossignol M, Ginibre C, Balabanidou V, Mendez NGA, Pigeon O, Vontas J, Cornelie S. Contributions of cuticle permeability and enzyme detoxification to pyrethroid resistance in the major malaria vector Anopheles gambiae. Sci Rep. 2017;7(1):1–10.

    Article  Google Scholar 

  29. Zoh MG, Bonneville J-M, Tutagata J, Laporte F, Fodjo BK, Mouhamadou CS, Sadia CG, McBeath J, Schmitt F, Horstmann S. Experimental evolution supports the potential of neonicotinoid-pyrethroid combination for managing insecticide resistance in malaria vectors. Sci Rep. 2021;11(1):1–16.

    Article  Google Scholar 

  30. Matowo J, Jones CM, Kabula B, Ranson H, Steen K, Mosha F, Rowland M, Weetman D. Genetic basis of pyrethroid resistance in a population of Anopheles arabiensis, the primary malaria vector in Lower Moshi, north-eastern Tanzania. Parasit Vectors. 2014;7(1):1–9.

    Article  Google Scholar 

  31. Zhou D, Duan B, Sun Y, Ma L, Zhu C, Shen B. Preliminary characterization of putative structural cuticular proteins in the malaria vector Anopheles sinensis. Pest Manag Sci. 2017;73(12):2519–28.

    Article  CAS  PubMed  Google Scholar 

  32. Balabanidou V, Kampouraki A, MacLean M, Blomquist GJ, Tittiger C, Juárez MP, Mijailovsky SJ, Chalepakis G, Anthousi A, Lynd A. Cytochrome P450 associated with insecticide resistance catalyzes cuticular hydrocarbon production in Anopheles gambiae. Proc Natl Acad Sci. 2016;113(33):9268–73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Abdalla H, Wilding CS, Nardini L, Pignatelli P, Koekemoer LL, Ranson H, Coetzee M. Insecticide resistance in Anopheles arabiensis in Sudan: temporal trends and underlying mechanisms. Parasit Vectors. 2014;7(1):1–9.

    Article  CAS  Google Scholar 

  34. David J-P, Strode C, Vontas J, Nikou D, Vaughan A, Pignatelli PM, Louis C, Hemingway J, Ranson H. The Anopheles gambiae detoxification chip: a highly specific microarray to study metabolic-based insecticide resistance in malaria vectors. Proc Natl Acad Sci. 2005;102(11):4080–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Nardini L, Christian RN, Coetzer N, Ranson H, Coetzee M, Koekemoer LL. Detoxification enzymes associated with insecticide resistance in laboratory strains of Anopheles arabiensis of different geographic origin. Parasit Vectors. 2012;5(1):1–12.

    Article  Google Scholar 

  36. Clark AG, Shamaan NA. Evidence that DDT-dehydrochlorinase from the house fly is a glutathione S-transferase. Pestic Biochem Physiol. 1984;22(3):249–61.

    Article  CAS  Google Scholar 

  37. Fournier D, Bride J-M, Hoffmann F, Karch F. Acetylcholinesterase. Two types of modifications confer resistance to insecticide. J Biol Chem. 1992;267(20):14270–4.

    Article  CAS  PubMed  Google Scholar 

  38. Grant DF, Hammock BD. Genetic and molecular evidence for a trans-acting regulatory locus controlling glutathione S-transferase-2 expression in Aedes aegypti. Mol Gen Genet MGG. 1992;234(2):169–76.

    Article  CAS  PubMed  Google Scholar 

  39. Kostaropoulos I, Papadopoulos AI, Metaxakis A, Boukouvala E, Papadopoulou-Mourkidou E. Glutathione S–transferase in the defence against pyrethroids in insects. Insect Biochem Mol Biol. 2001;31(4–5):313–9.

    Article  CAS  PubMed  Google Scholar 

  40. Vontas JG, Small GJ, Hemingway J. Glutathione S-transferases as antioxidant defence agents confer pyrethroid resistance in Nilaparvata lugens. Biochem J. 2001;357(1):65–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Ranson H, Paton M, Jensen B, McCarroll L, Vaughan A, Hogan J, Hemingway J, Collins F. Genetic mapping of genes conferring permethrin resistance in the malaria vector, Anopheles gambiae. Insect Mol Biol. 2004;13(4):379–86.

    Article  CAS  PubMed  Google Scholar 

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

    Article  CAS  Google Scholar 

  43. Sandeu MM, Mulamba C, Weedall GD, Wondji CS. A differential expression of pyrethroid resistance genes in the malaria vector Anopheles funestus across Uganda is associated with patterns of gene flow. PLoS ONE. 2020;15(11):e0240743.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Ishak IH, Kamgang B, Ibrahim SS, Riveron JM, Irving H, Wondji CS. Pyrethroid resistance in Malaysian populations of dengue vector Aedes aegypti is mediated by CYP9 family of cytochrome P450 genes. PLoS Negl Trop Dis. 2017;11(1):e0005302.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Mnzava AP, Knox TB, Temu EA, Trett A, Fornadel C, Hemingway J, Renshaw M. Implementation of the global plan for insecticide resistance management in malaria vectors: progress, challenges and the way forward. Malar J. 2015;14(1):1–9.

    Article  CAS  Google Scholar 

  46. Benedict M: MR4 Methods in Anopheles Research. Atlanta: CDC 2007.

  47. CDC CfDCaP: Guideline for evaluating insecticide resistance in vectors using the CDC bottle bioassay. Atlanta: The Centers; 2010.

  48. Collins FH, Mendez MA, Rasmussen MO, Mehaffey PC, Besansky NJ, Finnerty V. A ribosomal RNA gene probe differentiates member species of the Anopheles gambiae complex. Am J Trop Med Hyg. 1987;37(1):37–41.

    Article  CAS  PubMed  Google Scholar 

  49. Scott JA, Brogdon WG, Collins FH. Identification of single specimens of the Anopheles gambiae complex by the polymerase chain reaction. Am J Trop Med Hyg. 1993;49(4):520–9.

    Article  CAS  PubMed  Google Scholar 

  50. Andrews S: FastQC: a quality control tool for high throughput sequence data. 2010. 2016.

  51. Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):i884–90.

    Article  PubMed  PubMed Central  Google Scholar 

  52. Liao Y, Smyth GK, Shi W. The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote. Nucleic Acids Res. 2013;41(10):e108–e108.

    Article  PubMed  PubMed Central  Google Scholar 

  53. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, Subgroup GPDP. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Robinson MD, Stirzaker C, Statham AL, Coolen MW, Song JZ, Nair SS, Strbenac D, Speed TP, Clark SJ. Evaluation of affinity-based genome-wide DNA methylation data: effects of CpG density, amplification bias, and copy number variation. Genome Res. 2010;20(12):1719–29.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J Roy Stat Soc: Ser B (Methodol). 1995;57(1):289–300.

    Google Scholar 

  56. Conesa A, Götz S, García-Gómez JM, Terol J, Talón M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21(18):3674–6.

    Article  CAS  PubMed  Google Scholar 

  57. Hunter S, Apweiler R, Attwood TK, Bairoch A, Bateman A, Binns D, Bork P, Das U, Daugherty L, Duquenne L. InterPro: the integrative protein signature database. Nucleic Acids Res. 2009;37:D211–5.

    Article  CAS  PubMed  Google Scholar 

  58. Jones P, Binns D, Chang H-Y, Fraser M, Li W, McAnulla C, McWilliam H, Maslen J, Mitchell A, Nuka G. InterProScan 5: genome-scale protein function classification. Bioinformatics. 2014;30(9):1236–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Simma EA, Dermauw W, Balabanidou V, Snoeck S, Bryon A, Clark RM, Yewhalaw D, Vontas J, Duchateau L, Van Leeuwen T. Genome-wide gene expression profiling reveals that cuticle alterations and P450 detoxification are associated with deltamethrin and DDT resistance in Anopheles arabiensis populations from Ethiopia. Pest Manag Sci. 2019;75(7):1808–18.

    Article  CAS  PubMed  Google Scholar 

  60. Weill M, Malcolm C, Chandre F, Mogensen K, Berthomieu A, Marquine M, Raymond M. The unique mutation in ace-1 giving high insecticide resistance is easily detectable in mosquito vectors. Insect Mol Biol. 2004;13(1):1–7.

    Article  CAS  PubMed  Google Scholar 

  61. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, Handsaker RE, Lunter G, Marth GT, Sherry ST, et al. The variant call format and VCFtools. Bioinformatics. 2011;27(15):2156–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We would like to express our heartfelt gratitude to all individuals and organizations whose contributions and support have been invaluable in the successful completion of this research.

Our profound gratitude goes to the field and insectary team, consisting of Mathew Kipsum, Edward Esalimba, Duncan Omondi and Richard Amito who was the team lead, for their invaluable assistance in collecting larvae and raising them to adult stage. Their dedication and expertise greatly contributed to the success of this study. We also extend our sincere thanks to the entire Biotechnology Core Facility at CDC for sequencing our samples and specifically, Claudia Corredor for training and assisting in the wet lab analyses. We acknowledge the MR4 team at the Entomology Branch, CDC Atlanta for providing the susceptible colony for this study. We are immensely thankful to Kevin Owuor who generated the map used in this study. Furthermore, we deeply appreciate Grace Chumbe, Jenipher Adhiambo and Nicole Dzuris for facilitating all the administrative tasks pertaining to this study. Lastly, we gratefully acknowledge KEMRI Director General for the permission to publish this work.

Funding

This work was supported by the Bill and Melinda Gates Foundation, Grand Challenges Grant No. INV024969 and Investment Grant OPP1210769.

Author information

Authors and Affiliations

  1. Kenya Medical Research Institute (KEMRI), Centre for Global Health Research (CGHR), Kisumu, Kenya

    Diana Omoke, Stephen Okeyo & Eric Ochomo

  2. Entomology Branch, Division of Parasitic Diseases and Malaria, Centers for Disease Control and Prevention, Atlanta, USA

    Lucy Mackenzie Impoinvil, Dieunel Derilus, Audrey Lenhart & Luc Djogbénou

  3. University of Abomey Calavi (UAC), Abomey Calavi, Benin

    Helga Saizonou

  4. University of Cape Town ZA, Cape Town, South Africa

    Nicola Mulder

  5. School of Life Sciences, Arizona State University, Tempe, AZ, USA

    Nsa Dada

  6. Tropical Infectious Disease Research Center, University of Abomey- Calavi, Abomey Calavi, Benin

    Nsa Dada & Luc Djogbénou

Authors
  1. Diana Omoke
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Lucy Mackenzie Impoinvil
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Dieunel Derilus
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Stephen Okeyo
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Helga Saizonou
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Nicola Mulder
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Nsa Dada
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Audrey Lenhart
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Luc Djogbénou
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Eric Ochomo
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

EO, LD, AL, ND, LI, NM and DO: conceptualization of the study, methodology design, sample collection and revising the manuscript. DD and DO: data curation and software. DO, DD, LI, HS and SO: data analysis, interpretation, and validation as well as drafting and revising the manuscript.

Corresponding authors

Correspondence to Diana Omoke or Eric Ochomo.

Ethics declarations

Ethics approval and consent to participate

Approval for this study was granted by the Kenya Medical Research Institute (KEMRI), Scientific Ethics Review Unit steering committee (SERU 3275). Verbal informed consent was obtained from the head of each household prior to mosquito collections.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

Descriptive statistics of RNA-Seq reads and alignments.

Additional file 2.

Descriptive statistics of feature counts output.

Additional file 3.

The percentage of tags assigned to features in each replicate of alphacypermethrin, deltamethrin and primiphosmethyl experiments.

Additional file 4.

Venn diagrams showing the differentially expressed genes in the Res-Sus, Res-Con and Con-Sus pairwise comparisons from both Siaya and Migori.

Additional file 5.

Functional annotation of the overlapped DEG between Res-Sus, Con-Sus, or Res-Con comparisons for each experiment. All overlapped DEGs belonged to different intersection, for all the six experiments were mapped to their functional description and their fold change expression.

Additional file 6.

Bar plot comparing the number of up regulated and down regulated DEGs associated with CPs (Cuticular proteins), COEs (Carboxylesterases), CYPs (cytochrome P450s, GSTs (glutathione-S-transferases), and SGs (Salivary gland) proteins.

Additional file 7.

Functional description and fold change expression of the DEGs that are commonly shared between the insecticide resistant populations of the two sites for each insecticide.

Additional file 8.

Gene ontology enrichment analysis on the DEGs.

Additional file 9.

Heatmap showing the allele frequencies of kdr, ACE1 and GSTE2 genes in insecticide resistant samples.

Additional file 10.

Statistics of differentially expressed genes with annotations from VectorBase and Blast2go.

Additional file 11.

Concentrations of insecticides used for bioassays.

Additional file 12.

List of Oligonucleotide primers used for RT-qPCR.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. 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 in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Omoke, D., Impoinvil, L.M., Derilus, D. et al. Whole transcriptomic analysis reveals overexpression of salivary gland and cuticular proteins genes in insecticide-resistant Anopheles arabiensis from Western Kenya. BMC Genomics 25, 313 (2024). https://doi.org/10.1186/s12864-024-10182-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-024-10182-9

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us