Esterase-Mediated Pyrethroid Resistance in Populations of an Invasive Malaria Vector Anopheles stephensi from Ethiopia

1. Introduction

The rapid spread of the invasive malaria vector, Anopheles stephensi, across Africa poses a significant threat to malaria control efforts. Originating from Southeast Asia and the Arabian Peninsula, this mosquito species has rapidly colonized multiple African countries such as Djibouti, Ethiopia, Somalia, Kenya, Sudan, Nigeria and Ghana, leading to an increase in malaria cases [1,2]. Unlike the native African vectors, An. stephensi thrives in urban settings and has developed resistance to multiple classes of insecticides, including pyrethroids, organophosphates, and carbamates [3,4,5]. This widespread resistance undermines the effectiveness of traditional control measures such as long-lasting insecticidal nets (LLINs) and indoor residual spraying (IRS) [6,7]. To effectively manage insecticide resistance, it is crucial to monitor resistance patterns and understand the underlying mechanisms.

Insecticide resistance is primarily driven by two mechanisms: target-site insensitivity and metabolic resistance [8]. Target-site insensitivity, particularly mutations in the voltage-gated sodium channel (VGSC) gene, confers resistance (knockdown resistance, kdr) to pyrethroids and DDT. Mutations and duplications in the acetylcholinesterase 1 (Ace-1) gene have been associated with resistance to carbamates and organophosphates, while mutations in the GABA receptor (Rdl) gene confer resistance to phenylpyrazoles and organochlorines. Additionally, metabolic resistance, which often associated with overexpression of detoxification enzymes, can contribute to resistance to multiple insecticide classes. Other potential mechanisms, such as alterations in the mosquito microbiome and cuticle, may also play a role [9,10].

In An. stephensi, several target-site mutations associated with insecticide resistance have been detected. These include the kdr-West (L1014F) and kdr-east (L1014S) mutations in the voltage-gated sodium channel gene, the N177D mutation in the acetylcholinesterase 1 (Ace-1) gene, and the A296S mutation in the GABA receptor (Rdl) gene [11,12,13,14]. However, the frequency of these mutations is generally low. For instance, the kdr-West mutation has been reported at low frequencies in various populations [11,14,15,16]. Additionally, the V189L mutation in the GSTe2 gene has been detected at low frequency in Ethiopian mosquito populations [13]. These findings suggest that while target-site mutations contribute to insecticide resistance in An. stephensi, other mechanisms, such as metabolic resistance, may play a more significant role. Increased activity of detoxification enzymes, including glutathione S-transferases (GSTs), esterases, and cytochrome P450s, has been implicated in insecticide resistance [17]. Additionally, alterations in membrane transporters, such as ABC transporters, may contribute to reduced insecticide penetration [18,19]. Sanil et al. (2014) demonstrated a strong association between increased GST activity and resistance to multiple insecticide classes, including pyrethroids, organophosphates, organochlorines, carbamates, and biocides [20]. Recently, RNA-seq technology has been utilized to identify genes that are differentially expressed in An. stephensi larvae following permethrin exposure [21]. However, comprehensive whole-transcriptome analyses focused specifically on the molecular mechanisms of pyrethroid insecticide resistance in the newly invasive African field population of An. stephensi have not yet been conducted.

In this study, we performed standard World Health Organization (WHO) bioassays to evaluate insecticide susceptibility of a newly invasive field population of An. stephensi from Ethiopia. To investigate the molecular mechanisms underlying insecticide resistance, we conducted a whole-transcriptome analysis comparing gene expression profiles of resistant field populations with a laboratory-susceptible (STE2) strain. We further compared the transcriptomic profiles of untested field populations to identify shared patterns of gene expression associated with insecticide resistance

2. Materials and Methods

2.1. Mosquito Larval Sampling and Insecticide Susceptibility Bioassay

Anopheles mosquitoes were sampled from Awash Subac Kilo, a market town located in central Ethiopia at an altitude of approximately 1,000 meters. With a population of 34,669 (as of July 2021), this region experiences a hot, semi-arid climate characterized by minimal rainfall and a short-wet season from July to September. Both Plasmodium falciparum and P. vivax malaria parasites are present in this region. In June–July 2022, Anopheles larvae and pupae were collected from the field site and reared to the adult stage (3–5 days old). Prior to the insecticide susceptibility bioassay, the morphological identification of Anopheles individuals was performed while the mosquitoes were alive (a brief chill on ice if needed) following the taxonomic keys ) [22]. An. stephensi adult female mosquitoes were then used in the bioassay, where mosquitoes were exposed to pyrethroid insecticide-impregnated papers (0.05% deltamethrin and 0.75% permethrin) at diagnostic doses using WHO test tubes in a bio-secure insectary. Each insecticide was tested with batches of 25 female An. stephensi mosquitoes in both test and control tubes. At least 100 mosquitoes were exposed to each insecticide for 1 hour. Mortality was assessed 24 hours after exposure. Abbot's formula was applied to adjust for control mortality when needed. According to the WHO guidelines, mosquito populations were classified as 'resistant' if mortality rates were below 90%, as 'suspected resistant' if mortality rates were between 90% and 98%, and ‘susceptible” if mortality rates > 98% [8]. An. stephensi mosquitoes that survived the insecticide exposure test, along with unexposed samples, were promptly killed by placing them in a deep freezer for about 30 minutes until complete knockdown occurred. Samples were immediately stored in 0.5 ml Eppendorf tubes with RNALater for subsequent molecular and whole transcriptome analysis.

2.2. Molecular Identification of Anopheles Mosquito Species

DNA was extracted from a single leg of a subset of Anopheles mosquitoes to confirm species via qPCR. The highly sensitive TaqMan probe-based qPCR method was used to identify An. arabiensis and An. stephensi [23,24]. A total of 480 individuals were tested for molecular confirmation of species identification. Additionally, a subset of 20 individuals from each species underwent Sanger sequencing of ITS2 PCR products to further validate species identity [25].

2.3. RNA Extraction and Transcriptome Sequencing (RNA-seq)

RNA was extracted individually from An. stephensi female mosquito samples using column-based RNA extraction kits (cat # 2050, Zymo Research, CA, USA). The RNA was assessed for quality using agarose gel electrophoresis (1%) and quantified with the Qubit 3.0 RNA HS Assay (ThermoFisher). The RNA samples comprised three groups: field deltamethrin-resistant (AWR) mosquitoes, field unexposed control (AWK) mosquitoes, and a susceptible laboratory strain (STE2: BEI Resources, Cat# MRA-128). For each group, 10 RNA samples were pooled in equal amounts to create one composite sample. Two replicates/pools were prepared per group. In total, six pools, comprising 60 mosquitoes, were generated for RNA library preparation and RNA-seq. The IDT xGen RNA library prep kit with Illumina Ribo-zero Plus was used for library preparation. Final quality control was performed using Kapa qPCR and Agilent TapeStation. The libraries were sequenced on an Illumina® NovaSeq S4-PE150.

2.4. Data Filtering, Mapping Reads, and Transcriptome Analysis

The generated sequences underwent an initial filtering process. Raw paired-end sequence reads were assessed for quality using FASTQC (v0.12.1) [26], and the Trimmomatic (v.0.39) [27] was employed to remove Illumina adapters. After quality control and filtering, the sequence reads were analyzed using CLC Genomics Workbench v23 (CLC Bio, Cambridge, MD), with the reference genome VectorBase-68_AstephensiUCISS2018. The RNA-seq analysis was conducted using the “RNA-seq tool” and the reference genome annotated with genes and transcripts, applying the following parameters: Mismatch cost = 2, Insertion cost = 3, Deletion cost = 3, Length fraction = 0.8, Similarity fraction = 0.8, Maximum number of hits per read = 10, Count paired reads as two = No, Ignore broken pairs = Yes, and Expression value = RPKM. The Differential Expression in Two Groups tool in CLC was used to conduct a statistical differential expression analysis for a set of expression tracks and a control. The experiments included two biological replicates for both the field population group and the lab strain control group. The raw gene expression matrix was transformed into a normalized gene expression matrix. The biological coefficient of variation (BCV), representing biological variability, was estimated using edgeR, which models RNA-seq data variability through a negative binomial distribution [28]. Differentially expressed genes (DEGs) were identified using the "Differential Expression in Two Groups" tool in CLC, with the STE2 group used as the control. All RNA-seq data were screened for False Discovery Rate (FDR), and results were deemed valid if the FDR was less than 0.05 [29]. DEGs were determined based on a |fold-change| > 1.5 and false discovery rate (FDR) < 0.05. ShinyGo bioinformatics tool was used for functional annotation and pathway analysis of the DEGs [30].

2.5. Gene Co-Expression Network Construction with WGCNA

Gene co-expression networks analysis were conducted using the WGCNA package (version 1.72) in R based on differentially expressed genes (DEGs) [31]. The soft threshold (power) was determined using the “pickSoftThreshold” function to generate a scale-free network for both population datasets. The cluster dendrogram and co-expression gene module tree for DEGs were created using the "Dynamic Tree Cut" and "Merged Dynamic" methods. The module eigengene (ME), representing the expression profiles of the module genes, was calculated using WGCNA. Gene co-expression networks were then independently constructed with the blockwise Modules function using these parameters: power = 28, minModuleSize = 30, deepSplit = 2, and mergeCutHeight = 0.3. Edge and node data for the target module were exported with the "export Network To Cytoscape" function to build a co-expression network, which was visualized using Cytoscape software [32].

3. Results

3.1. Species Composition and Pyrethroid Susceptibility of An. Stephensi

Molecular identification of the Anopheles species indicated that 88.3% were An. stephensi and 11.7% were An. arabiensis (n=480), with no other Anopheles species identified. High levels of pyrethroid insecticide resistance were detected in the An. stephensi field population, with mortality rates falling below 60% for both tested insecticides. The mortality rate for permethrin was 59.3% ± 26.3% (n = 176), while for deltamethrin, it was 39.4% ± 18.9% (n = 250) (Figure S1).

3.2. RNA Sequencing and Preliminary Analysis of the Raw Data

RNA-seq was performed from three different groups: field resistant (AWR), field unexposed (AWK), and lab strain (STE2), with two independent biological replicates per group. The RNA used for library construction had to pass quantity and quality control criteria. Six libraries were constructed from the three different groups, and raw reads were produced from the six libraries using the Illumina sequencing platform. RNA-seq generated a total of 70–85 million reads per sample, with an average of 80.6 ± 23.6 million reads. After trimming, 57.5% ± 0.7% of paired end reads, 19.5% ± 0.2% of forward reads, and 9.4% ± 0.3% of reverse reads were retained, while 14.1% ± 0.5% of reads were discarded (Table 1).

3.3. Differentially Expressed Gene Analyses

The cleaned, high-quality RNA-seq reads were used to evaluate differences in gene expression (Figure 1), with FPKM values calculated to quantify expression levels across all genes in the three groups (6 pooled samples). The raw gene expression matrix was transformed into a normalized gene expression matrix (Figure S2). A total of 10,568 protein-coding genes (83.2% of the 12,705 total genes) were expressed at a minimum exon read count of 5 across the six pooled samples (Table S1). Principal component analysis (PCA) of the RNA-seq data revealed greater variation between groups than within groups (Figure 2A). The estimated coefficient of variation was 0.0517, corresponding to a biological coefficient of variation (BCV) of 0.2273. This BCV value is lower than that typically observed in other studies for well-controlled experiments (BCV = 0.4)(Figure 2B) [33], indicating the stability and reliability of the RNA-seq data.

Differential gene expression analysis identified 885 significantly differentially expressed genes (DEGs) of the two groups (Table S2). Specifically, 560 DEGs were found between AWR and STE2, with 168 downregulated and 392 upregulated genes (Figure 3A). Similarly, 490 DEGs were identified between AWK and STE2 (Figure 3B). These results provide insights into the molecular mechanisms underlying phenotypic differences between the two groups.

3.4. Classification of the Differentially Expressed Genes

Cluster analysis of the 885 differentially expressed genes revealed distinct patterns of gene expression of the two groups (Figure S3). Approximately one-third of the genes (165) were shared between the two groups, while the remaining two-thirds (395 and 325) were unique to each group, indicating significant genetic divergence among the populations (Figure S4).

A total of 20 genes involved in metabolic detoxification (P450s, GSTs, and esterases) were consistently differentially expressed between resistant, untested field, and susceptible laboratory An. stephensi populations (Table 2). Notably, genes from the esterase family, including venom carboxylesterase-6 and esterase E4, exhibited the most significant upregulation (up to 70-fold) in resistant populations compared to the susceptible strain. Additionally, four P450 genes and three UGT genes were consistently upregulated. In contrast, a few genes, such as chorion peroxidase and protein arginine N-methyltransferase 5, were downregulated (Table 2).

3.5. Functional Annotation of the Differentially Expressed Genes

Gene Ontology (GO) enrichment analysis was conducted to functionally annotate the differentially expressed genes (DEGs) based on the An. stephensi SDA-500 annotation. The GO analysis identified a total of 27 enriched GO categories, which include nine biological processes, one cellular components, and 17 molecular functions (Table S3). The most significantly enriched categories were in the catalytic activity under molecular functions and metabolic process under biological process (Figure 4). Furthermore, the pathway network analysis of functional linkages highlighted the top 20 GO enrichments, including GO:0003824 catalytic activity, GO:0016491 oxidoreductase activity, GO:0016787 hydrolase activity, GO:0004252 serine-type endopeptidase activity, GO:0008236 serine-type peptidase activity, GO:0017171 serine hydrolase activity, GO:0005576 extracellular region, GO:0008252 nucleotidase activity, GO:0008061 chitin binding, GO:0004175 endopeptidase activity, GO:0016788 hydrolase activity acting on ester bonds, GO:1901564 organonitrogen compound metabolic proc., GO:0008152 metabolic process, GO:0052689 carboxylic ester hydrolase activity, GO:0004497 monooxygenase activity, GO:0055085 transmembrane transport, GO:0005506 iron ion binding, GO:0022857 transmembrane transporter activity, GO:0005215 transporter activity, GO:0006575 cellular modified amino acid metabolic proc (Figure 5).

3.6. Construction of Gene Co-Expression Network

A soft threshold (power) of 28 was applied to construct a scale-free network for the differentially expressed genes (Figure S5A and Figure S5B). From the 560 differentially expressed genes between the resistant field population and the susceptible lab strain, five gene co-expression modules were identified from the original 20 modules (Figure 6). The blue and green modules exhibited consistent expression patterns across replicates (Figure S6). Most genes in the green module were upregulated, while those in the blue module were downregulated. A total of 48 genes in each module showed co-expression (Table S4 and Table S5). Figure 7 presents the co-expression network of genes associated with pyrethroid insecticide resistance in An. stephensi, with edge weight thresholds of 0.4 and 0.68, respectively.

4. Discussion

The recent invasion of An. stephensi in the Horn of Africa poses a significant threat to malaria control efforts. To effectively combat this invasive vector, it is crucial to assess it susceptibility to insecticides and understand the mechanisms underlying its insecticide resistance. This study aimed to investigate the molecular basis of insecticide resistance in An. stephensi populations in Ethiopia. Through RNA-seq analysis, we identified significant upregulation of esterase E4 and venom carboxylesterase-6 in the resistant field populations compared to a laboratory-susceptible strain. While the kdr mutation was reported at low frequencies, our findings suggest that metabolic resistance, mediated by detoxification enzymes, could be a primary driver of insecticide resistance in An. stephensi. These insights emphasize the need for implementing insecticide resistance management strategies which include integrated vector management to combat this growing challenge.

Esterases are enzymes that catalyze the hydrolysis of ester bonds. Esterase E4 has been previously implicated in insecticide resistance in other insects [34]. Recent studies have highlighted the potential role of venom carboxylesterase-6-like genes in insecticide resistance. For instance, Wang et al. (2023) identified a potential insecticide resistance marker within this gene in Bemisia tabaci [35]. Additionally, Kang et al. (2021) identified venom carboxylesterase-6-like as a critical odorant-degrading enzyme (ODE) in Aphidius gifuensis [36]. In our study, we observed significantly higher expression of venom carboxylesterase-6-like (a 52- to 70-fold increase) and Esterase E4 (a 26- to 33-fold increase) in the field population. This suggests that these genes may play a role in the detoxification of pyrethroid insecticides. Additionally, several other metabolic detoxification genes were found to be overexpressed, including monooxygenases (e.g., CYP6D3), glutathione S-transferases (e.g., GST1), and immunity genes (e.g., scaf, SP24D, SND1).

Next-generation sequencing (NGS), particularly RNA-seq, has revolutionized transcriptomic analysis by providing a comprehensive and unbiased view of gene expression [37,38]. Despite advancements, challenges persist, especially with organisms that have high rRNA content, such as archaea and certain insects [39,40,41]. Although RNA-seq protocols have greatly improved, residual rRNA contamination can still affect data quality, particularly in studies using kits developed for mammalian samples, where rRNA removal may be less specific in non-target organisms like mosquitoes. In our study, we successfully employed RNA-seq to identify differentially expressed genes in An. stephensi mosquitoes despite these limitations. While the Illumina Ribo-Zero Plus rRNA Depletion Kit proved effective, residual rRNA contamination was observed, highlighting the need for optimized protocols for non-standard organisms.

In Ethiopia, An. stephensi populations developed resistance to various insecticides. Several target-site mutations, including kdr-West (L1014F) in the voltage-gated sodium channel, N177D in acetylcholinesterase 1 (Ace-1), and A296S in the GABA receptor (Rdl), were detected in these populations [11,12,13]. While our RNA-seq analysis provides valuable insights into gene expression patterns, it was not optimal for SNP detection due to the limitations of pooled RNA-seq, such as reduced genomic coverage and potential allele frequency bias. To address these limitations, future studies should focus on high-quality RNA extraction and increased sequencing depth to enable accurate SNP detection and analysis.

5. Conclusions

The rapid spread of the invasive malaria vector, Anopheles stephensi, poses a significant threat to malaria control efforts in Africa. The emergence of insecticide resistance in this species further complicates vector control strategies. Our RNA-seq analysis revealed that An. stephensi populations in Ethiopia exhibited elevated expression of genes associated with pyrethroid resistance, including esterases (e.g., esterase E4 and venom carboxylesterase-6), P450s, and GSTs. These findings suggest that metabolic resistance plays a crucial role in the development of insecticide resistance in this species. Functional validation of key genes involved in insecticide resistance, coupled with the development of improved resistance markers, will be essential for monitoring resistance trends and guiding future control efforts.