Identification of Novel Gene Clusters Potentially Associated with Insecticide Resistance in Anopheles gambiae s.l.

1. Introduction

Malaria remains one of the leading causes of morbidity and mortality in sub-Saharan Africa, especially among children under five years and pregnant women [1]. Despite recent advances in malaria vaccine development, insecticide-based interventions, including long-lasting insecticidal nets (LLINs) and indoor residual spraying (IRS), remain the most effective prevention methods [2]. However, the success of these interventions depends on a good understanding of vector population dynamics, particularly their susceptibility to insecticides, to guide the deployment of these tools. In mosquitoes, two major mechanisms of insecticide resistance are target-site and metabolic resistance [3,4,5]. Target site resistance involves genetic changes that alter the protein targeted by the insecticide, while metabolic resistance results from changes in the sequence or expression of a complex array of enzymes and detoxification pathways that remove the insecticide from the system [6]. In addition to these mechanisms, behavioral and cuticular resistance, have also been described in malaria vectors [7]. The rise of resistance remains a major concern for malaria control and calls for continued resistance surveillance and rapid development of alternative control tools to accelerate malaria elimination.

Although target-site resistance has been extensively described in the African malaria vector An. gambiae s.l., metabolic resistance remains poorly characterized [8], likely due to the number of gene families involved in detoxification, redundancy among their members, and the multiple mechanisms by which metabolic resistance can arise [9]. However, investigations on this mechanism are becoming increasingly crucial, especially as recent studies demonstrated the growing contribution of the detoxification genes to the resistance phenotype observed in Anopheles gambiae s.l. [10,11,12]. This trend is likely driven by new mutations within genes encoding detoxifying enzymes that are associated with resistance phenotypes in members of the An. gambiae complex [13]. Enzymes that detoxify insecticides include cytochrome P450s, glutathione-S-transferases and carboxylesterases [14].

One of the major challenges in the ongoing fight against malaria is the identification of genomic bases of resistance, particularly in ways that can inform vector control programs in the field [15]. Traditionally, insecticide resistance has been monitored using phenotypic assays such as WHO tube tests or CDC bottle bioassays. While useful, these methods are labor-intensive, time-consuming, and require substantial logistical resources. In contrast, advances in genomics now offer powerful alternatives, enabling the identification and surveillance of molecular markers associated with resistance in natural mosquito populations [16]. To date, most genomic studies have focused on specific regions of the genome already known to be linked with resistance. However, given the dynamic and rapidly evolving nature of insecticide resistance, including the frequent emergence of novel resistance mechanisms, whole genome analysis represents an optimal investigative approach. Genome-wide scans for signatures of adaptive evolution, such as selective sweeps, provide a comprehensive framework to detect emerging resistance loci beyond known targets.

In this study, we used genome-wide scans for selection (GWSS) to detect genomic regions under signals of selection in Anopheles gambiae s.l. populations across sub-Saharan Africa. A literature review and genetic variability analyses were conducted to identify, within these regions and populations, aldehyde oxidase gene family, a novel candidate associated with insecticide resistance.

2. Materials and Methods

2.1. Data Collection and Analysis

The genomic data used in this study are from the Anopheles gambiae 1000 Genomes phase 3 Project (Ag1000G) called Ag3.0, published in February 2021 [17]. Mosquitoes, including three Anopheles species, Anopheles gambiae, Anopheles coluzzii and Anopheles arabiensis, were collected between 2004 and 2015 in 11 countries across sub-Saharan Africa (Figure 1, Table S1). Details on mosquito sample collection, storage and management of the genomic data, including access rights are available on the homepage of MalariaGEN [18].

Library preparation and sequencing were performed at the Wellcome Sanger Institute. All individuals were sequenced to a target coverage of 30× using the Illumina HiSeq 2000 and the Illumina HiSeq X platforms [18].

Alignment, SNP calling and sample quality control (QC), were performed by the MalariaGEN Resource Centre Team. Reads were aligned to the AgamP4 reference genome using BWA version 0.7.15. GATK version 3.7-0 RealignerTargetCreator and IndelRealigner were used for Indel realignment while single nucleotide polymorphisms were called using GATK version 3.7-0 UnifiedGenotyper. Given all possible alleles at all genomic sites where the reference base was not “N”, genotypes were called for each sample independently. After variant calling, both the samples and the variants have undergone a series of quality control analyses, to ensure data quality [18].

2.2. Detection of Genome Regions Under Recent Positive Selection and Identification of New Candidate Genes

Genome-wide scans for selection (GWSS) were performed on Anopheles gambiae s.l. to detect signals of recent selection using the Garud H12 statistic [16]. This statistic is sensitive to recent selection, making it ideal for detecting selective sweeps driven by recent insecticidal pressures. The analyses were first conducted at the country level, then at both country and species levels. For each analysis, the window size, defined by a specific number of SNPs, was calibrated to balance resolution and signal detection. Genomic regions under selective sweep were examined for insecticide resistance genes. Potential new candidates were identified through literature research using resources such as VectorBase, NCBI, and Google Scholar.

2.3. Genetic Variability

For each aldehyde oxidase gene, the transcript was identified. Single nucleotide polymorphisms (SNPs) were identified and grouped by their predictive effect. SNP frequencies within populations were determined, then filtered down to just non-synonymous SNPs that were at the highest frequencies. Copy number variation (CNV) was assessed by computing the frequencies of gene amplification and deletion across the aldehyde oxidase gene regions in the different populations.

Haplotype networks were built for the aldehyde oxidase genes and explored to find evidence of gene flow events between countries and species. Haplotypes carrying the main high frequency mutations were analyzed with a maximum genetic distance of two SNPs.

Population structure was analyzed using principal component analysis (PCA) and pairwise Fst estimates. Genetic diversity summary statistics were computed using SNPs called in the aldehyde oxidase genes region.

3. Results

3.1. Signals of Positive Selection - Identification of Potential New Insecticide Resistance Genes

The analyses included a total of 1908 individual genomes from An. gambiae s.s., An. coluzzii and An. arabiensis. Genome-wide selection was conducted using the Garud H12 statistics in windows ranging from 1000 to 2000 SNPs depending on the population. The results revealed several peaks of H12 values across different chromosome arms. These peaks correspond to the genomic regions of Vgsc, Ace1, Rdl and GSTe3 genes, indicating selective sweeps likely due to positive selection in these areas, as previously observed in An. gambiae populations.

Beyond these well-characterized genes, high H12 values were observed in a new locus on chromosome 2L around positions 28,510,000 and 28,590,000, indicating a signal of positive selection in this region. This genomic region corresponds to a cluster of genes encoding detoxification enzymes, including five aldehyde oxidases (AGAP006220, AGAP006221, AGAP006224, AGAP006225, AGAP006226), two glucosyl/glucuronosyl transferases (AGAP006222, AGAP006223), and two carboxylesterases (AGAP006227, AGAP006228). Additional population-specific scans were performed to determine the geographic and species-specific distribution of this selective sweep. In Anopheles gambiae s.s, the sweep was detected in Burkina Faso, Ghana and Guinea. In Anopheles arabiensis it was present in Tanzania and Malawi, but with relatively lower proportions. However, in An. coluzzii the sweep was completely absent (Figure 2).

We focused our further investigations on aldehyde oxidases, as three of them were located directly under the peak of the selective sweep and based on existing literature, appeared particularly promising as candidate insecticide resistance genes.

3.2. SNPs in the Aldehyde Oxidase Genes

Single Nucleotide Polymorphism (SNP) variation was analysed in aldehyde oxidase genes across the five countries, Burkina Faso, Ghana, Guinea, Tanzania and Malawi where these genes were found to be under positive selection. All five genes exhibited a high density of polymorphism across the different populations. Specifically, we identified 2,202 SNPs (including 569 non-synonymous coding SNPs) in AGAP006220, 2,661 SNPs (691 non-synonymous) in AGAP006221, 2,698 SNPs (1,433 non-synonymous) in AGAP006224, 2,515 SNPs (978 non-synonymous) in AGAP006225, and 2,104 SNPs (811 non-synonymous) in AGAP006226. The analyses identified a substantial number of non-synonymous coding SNPs at high frequencies across all five genes and in all three species, with some reaching 100% frequency in several cohorts. To focus on the most prevalent variants, the dataset was refined to include only those SNPs with minimum frequencies of 40, 50, or 75% in at least one sample cohort. The refinement showed 21 SNPs at frequencies higher than 40% in AGAP006220 and AGAP006225, and 56 SNPs above 75% in AGAP006221, AGAP006224 and AGAP006226 (Figure 3, Figures S1–S4). These patterns suggest that these variants might be under strong selective pressure and could play an important role in the adaptation of malaria vectors to environmental changes such as insecticide resistance.

3.3. CNVs in the Aldehyde Oxidase Genes

CNV frequencies were analysed in the five aldehyde oxidase genes across populations from Burkina Faso, Ghana, Guinea, Tanzania, and Malawi. The analysis revealed the presence of CNVs in all five genes. Deletions were detected only in AGAP006224, while amplifications were observed in all five genes. Overall, amplification frequencies were low (≤3%) across most of the genes and populations, with the exception of Anopheles arabiensis from Tanzania, where AGAP006225 and AGAP006226 exhibited amplification rates of 10.59% and 11.76%, respectively (Figure 4).

3.4. Gene Flow, Population Structure, and Genetic Diversity at the Aldehyde Oxidases Locus

To investigate gene flow at the aldehyde oxidases locus, haplotype network technique was used to visualize the relationships between haplotypes from different species and countries and show how variants are shared between populations. The results showed strong signs of adaptive gene flow between Anopheles gambiae s.l. populations in West Africa, as well as between populations of East Africa. Specifically, in West Africa, the largest network comprised 163 identical haplotypes shared across all the vector populations. Similarly, the same pattern is observed in East Africa, with 52 identical haplotypes shared between An. gambiae s.l. populations from Tanzania and Uganda (Figure 5). Inter-species gene flow was also evident with a large haplotype cluster of 1001 identical sequences shared between Anopheles gambiae (n = 609) and Anopheles coluzzii (n = 392), indicating strong evidence for adaptive gene flow between these two species. However, no shared haplotypes were detected between An. arabiensis and either An. gambiae or An. coluzzii (Figure 6).

These patterns were supported by population structure and genetic diversity analyses conducted at the aldehyde oxidases locus across the 11 countries. The principal component analysis (PCA) showed a distinct and homogeneous cluster of An. arabiensis, whereas An. gambiae and An. coluzzii clustered together, further supporting genetic exchange between these two species (Figure 8). The genetic differences (F_ST) between An. arabiensis and the two other species were generally higher than those observed between An. gambiae and An. coluzzii, indicating greater genetic differentiation (Figure S5). The results were confirmed by diversity indices, although some variation existed within species. Both measures of diversity (θπ and θw) showed similar results for An. gambiae and An. coluzzii compared to An. arabiensis (Figure 7). These results suggest a similar demographic history for An. gambiae and An. coluzzii.

Figure 7. Genetic diversity statistics of An. gambiae s.l. populations at aldehyde oxidases locus. A. Nucleotide diversity (θπ). B. Watterson theta (θw). The X-axis shows the different cohorts (country/sampling period/species): AGO=Angola, BFA=Burkina Faso, CIV=Côte d’Ivoire, CMR=Cameroon, COD=Democratic Republic of Congo, GHA=Ghana, GIN=Guinea, MLI=Mali, MWI=Malawi, TZA=Tanzania, UGA=Uganda.

Figure 7. Genetic diversity statistics of An. gambiae s.l. populations at aldehyde oxidases locus. A. Nucleotide diversity (θπ). B. Watterson theta (θw). The X-axis shows the different cohorts (country/sampling period/species): AGO=Angola, BFA=Burkina Faso, CIV=Côte d’Ivoire, CMR=Cameroon, COD=Democratic Republic of Congo, GHA=Ghana, GIN=Guinea, MLI=Mali, MWI=Malawi, TZA=Tanzania, UGA=Uganda.

Figure 8. Principal component analysis of An. gambiae s.l. populations at aldehyde oxidases locus. X-axis shows the first principal component (PC1) and Y-axis the second principal component (PC2).

Figure 8. Principal component analysis of An. gambiae s.l. populations at aldehyde oxidases locus. X-axis shows the first principal component (PC1) and Y-axis the second principal component (PC2).

4. Discussion

Malaria control remains a major challenge in sub-Saharan Africa. Despite ongoing interventions, the region still bears the highest global burden of the disease due to complicating factors such as insecticide resistance, changing vector behaviour and rapid urbanization [2,19]. However, strengthening surveillance, particularly through genomic tools, offers new opportunities to enhance the effectiveness of vector control strategies [20]. In particular, discovering new resistance genes is crucial for controlling malaria effectively by enabling early detection and monitoring within malaria vector populations. These insights help guide and inform the selection and rotation of insecticides, reducing the spread of the resistance variants [21]. Recent advances in genomics offer valuable opportunities to study the genomic architecture and evolution of malaria vectors, enabling high-resolution mapping of resistance genes, population structure, and adaptive evolution [21,22].

Our study detected a signal of positive selection within the aldehyde oxidase genes in An. gambiae and An. arabiensis populations across sub-Saharan Africa. The absence of positive selection in An. coluzzii supports recent findings on the carboxylesterases Coeae1f and Coeae2f, which, as previously mentioned, are located on the same locus [23]. Aldehyde oxidases are indeed known for their broad substrate specificity and marked species differences [24].

Our results suggest that aldehyde oxidase genes may play an important role in how malaria vectors adapt to environmental changes, such as those caused by insecticide-based vector control tools. Aldehyde oxidase is a cytosolic enzyme that belongs to the family of structurally related molybdo flavoproteins like xanthine oxidase. These enzymes have originated from gene duplication events of an ancestral eukaryotic xanthine oxidoreductase gene found in the flatworm, Caenorhabditis elegans [24]. Aldehyde oxidases are common enzymes involved in detoxifying and breaking down many substances across various organisms. Although references to aldehyde oxidase and its role in metabolism date back to the 1950s [24], the significance of this enzyme in insecticide metabolism has only emerged in recent years [25]. Like cytochromes P450, aldehyde oxidases contribute to oxidation by carrying out nucleophilic attacks using oxygen derived from water [26]. Previous studies have characterised these enzymes in Culex mosquitoes. Investigations into resistant Culex strains showed elevated levels of aldehyde oxidase activity, suggesting that these enzymes may confer a selective advantage in the presence of insecticides [27,28]. Beyond Culex species, growing evidence indicates that aldehyde oxidases may also contribute to insecticide resistance in An. gambiae s.l. In fact, previous studies have shown that aldehyde oxidases are co-amplified with insecticide resistance-associated esterases. These esterases have been recently demonstrated in association with resistance to pirimiphos-methyl in An. gambiae s.l. populations [23,27]. Furthermore, aldehyde oxidases might play a role in resistance to neonicotinoid insecticides. In vitro studies have shown their capacity to metabolize neonicotinoids such as imidacloprid, thiamethoxam, clothianidin, and dinotefuran [29],indicating their potential involvement in resistance to this newer class of insecticides. This is particularly relevant given the growing interest in neonicotinoids as alternative tools for vector control in areas where resistance to multiple insecticide classes is already widespread [30].

These observations are consistent with our results, which, in addition to positive selection, showed increased frequency of non-synonymous single nucleotide polymorphism (SNP) and copy number variations (CNVs), including gene amplifications in the genomic region of aldehyde oxidase genes. Indeed, SNPs and CNVs have been identified as playing significant role in the evolution of insecticide resistance [31,32,33,34]. These variants are particularly important in metabolic resistance because they can lead to increased gene expression [3,35]. However, given the low frequencies of CNVs observed, we hypothesize that the resistance ultimately induced by these detoxification genes would be more closely linked to SNPs, which were detected at high frequencies. This hypothesis is supported by previous studies showing that SNPs are a key genetic factor influencing variation in aldehyde oxidase expression and active protein levels [24].

The adaptive introgression observed between An. gambiae and An. coluzzii aligns with previous studies demonstrating a high potential for genetic exchange between these two closely related species [22]. Although no significant selective sweep was detected in An. coluzzii, the ongoing and increasing gene flow suggests that aldehyde oxidase genes may eventually come under selective pressure in this species as well. Given the spatial patterns of gene flow observed in both study regions, these putative resistance markers are expected to expand to large areas across the continent.

This study highlighted the ongoing evolution of the An. gambiae complex genome in response to environmental pressures and vector control interventions. Genomic changes, including mutations and gene duplications, are driven by factors such as exposure to insecticides, environmental changes and ecological shifts. These evolutionary processes can lead to increased resistance and altered behaviour, which complicate malaria control efforts. Therefore, continuous genomic surveillance is crucial to track and effectively respond to these changes.

5. Conclusions

This study highlighted the importance of genomic data in malaria surveillance. Genomic surveillance can identify how insecticide resistance develops in malaria-transmitting mosquitoes and help guide the deployment of new vector-control tools with multiple modes of action. The results of this study highlighted a cluster of genes that may be associated with insecticide resistance in An. gambiae s.l. populations in Africa. Indeed, the study identified a previously unreported locus, consisting of detoxifying genes encoding aldehyde oxidases, under recent positive selection in An. gambiae s.l. populations. While the involvement of aldehyde oxidases in insecticide resistance has been studied in Culex mosquitoes, our findings provide the first hypothesis supporting the potential role of these genes in the adaptation of An. gambiae s.l to the environmental changes. Genetic variation analyses further support this hypothesis by revealing both copy number variations (CNVs) and high-frequency non-synonymous SNPs within these genes.

Further studies, including gene expression and functional validation, are needed to confirm the role of aldehyde oxidases in insecticide resistance in An. gambiae s.l. If proven, these genes could serve as novel molecular markers to be integrated into resistance monitoring programs, contributing to more effective and adaptive vector control strategies.