Key gene modules and hub genes associated with pyrethroid and organophosphate resistance in Anopheles mosquitoes: a systems biology approach

Abstract

Indoor residual spraying (IRS) and insecticide-treated nets (ITNs) are the main methods used to control mosquito populations for malaria prevention. The efficacy of these strategies is threatened by the spread of insecticide resistance (IR), limiting the success of malaria control. Studies of the genetic evolution leading to insecticide resistance could enable the identification of molecular markers that can be used for IR surveillance and an improved understanding of the molecular mechanisms associated with IR. This study used a weighted gene co-expression network analysis (WGCNA) algorithm, a systems biology approach, to identify genes with similar co-expression patterns (modules) and hub genes that are potential molecular markers for insecticide resistance surveillance in Kenya and Benin. A total of 20 and 26 gene co-expression modules were identified via average linkage hierarchical clustering from Anopheles arabiensis and An. gambiae, respectively, and hub genes (highly connected genes) were identified within each module. Three specific genes stood out: serine protease, E3 ubiquitin-protein ligase, and cuticular proteins, which were top hub genes in both species and could serve as potential markers and targets for monitoring IR in these malaria vectors. In addition to the identified markers, we explored molecular mechanisms using enrichment maps that revealed a complex process involving multiple steps, from odorant binding and neuronal signaling to cellular responses, immune modulation, cellular metabolism, and gene regulation. Incorporation of these dynamics into the development of new insecticides and the tracking of insecticide resistance could improve the sustainable and cost-effective deployment of interventions.

Keywords: Anopheles arabiensis; Anopheles gambiae; Hub genes; Insecticide resistance; Molecular markers.

Fig. 1

Representation plots of the total read counts for An. gambiae (A) and An. arabiensis (B). Fig 1A represents a dendrogram plot for An. gambiae normalized read counts. DA = Djougou Alphacypermethrin, DD = Djougou Deltamethrin, BP = Bassila Permethrin, BA = Bassila Alphacypermethrin, and BP = Bassila Primiphos-methyl. Fig. 1B represents a dendrogram plot for An. arabiensis normalized read counts. MA0 = Migori Alphacypermethrin, MP0 = Migori Primiphos-methyl, MD = Migori Deltamethrin, SA0 = Siaya Alphacypermethrin, SD0 = Siaya Deltamethrin, SP0 = Siaya Primiphos-methyl,KIS = Kisumu strain An. gambiae, and DON = Dongola strain An. arabiensis

Fig. 2

Scale-free topology fitting index (R2) versus a soft threshold power plots and hierarchical clustering dendrogram plots for both An. gambiae and An. arabiensis

Fig. 3

Module trait relationship heat map plots and module preservation summary plots in An. gambiae and An. arabiensis.Each point represents a module labeled by color.The dashed lines indicate thresholds Z = 2 (no preservation) and Z = 10 (highly preserved modules)

Fig. 4

An. gambiae and An. arabiensis preserved module expression plots in resistant and susceptible anopheles populations.Higher values are represented with reds of increasing intensity, and lower values are represented with greens of increasing intensity

Fig. 5

Resistant anopheles gene co-expression networks generated using weighted gene co-expression network analysis based on TOM > 0.1 for visualization. This figure illustrates gene co-expression networks; each node (point) represents a gene and genes of the same color form modules. The edges (lines) connecting the nodes represent gene-to-gene relationships. Fig. 5A represents the An. gambiae gene co-expression network, whereas Fig. 5B illustrates the An. arabiensis gene co-expression network

Fig. 6

Panoramic view of molecular interactions in insecticide-resistant Anopheles mosquitoes. OBP-Odorant Binding Proteins; ORN-Olfactory Receptor Neurons; GPCRs-G coupled protein receptors; PKA-Protein Kinase A

Fig. 7

Gene expression profiles of resistant An. arabiensis and An. gambiae from (A) Kenya and (B) Benin exposed to deltamethrin, primiphos-methyl and alphacypermethrin compared to the susceptible An. arabiensis Dongola strain and An. gambiae Kisumu strain, respectively. The horizontal dotted line on the volcano plot denotes a P-value of 0.01, while the vertical dotted lines indicate twofold expression differences