Distinct Genetic Footprint of Adaptation Observed Within Anopheles gambiae Taxa from The Gambia

Introduction

The Anopheles gambiae complex (Diptera: Culicidae) is a primary malaria vector in sub-Saharan (SSA).¹ The complex comprises nine sibling species, namely Anopheles arabiensis, Anopheles coluzzii, An. gambiae s.s, Anopheles melas, Anopheles merus, An. quadriannulatus, Anopheles bwambae, Anopheles amharicus and Anopheles fontenillei, that are morphologically identical but genetically dissimilar.^2–4 The dissimilarity is, in part, a consequence of local adaptation to survive and breed in diverse environments. When adapting to new environments, beneficial alleles may accumulate in the genome and become fixed in the mosquito population over time.⁵ Genomic regions harbouring these alleles, detected as signatures of selection (ie, selective sweeps), can offer insight into mosquito evolutionary history and genetic diversity.⁶ Indeed, selective sweeps have been explored to understand resistance mechanisms in malaria mosquitoes, underlying adaptation to vector control measures:⁷ The most reported are sweeps that host metabolic (ie cytochrome P450 and GSTe-2) and target site (ie Vgsc, Rdl and Ace-1) resistance genes.^8–10

In West Africa, the most commonly reported target site mutations in An. gambiae complex includes kdr-w (L1014F) and N1575Y in the Vgsc gene which are associated with resistance to pyrethroids and DDT.^11,12 The others include Ace-1 G119S mutation associated with resistance to carbamates and organophosphates, as well as the A302G, A302S and A296S-RDL mutations associated with resistance to dieldrin.^11,13–15 For metabolic resistance, commonly reported variants are in cytochrome P450 CYP6P3, CYP6M2, and CYP9K1 and these are associated with resistance to the pyrethroid deltamethrin.^16,17 Additionally, it is not unusual to observe the co-occurrence of multiple resistance mechanisms, including L119F-GSTe2 that confers resistance to DDT.^18,19 Importantly, additional and new resistance-associated variants continue to emerge across African malaria vector species driving phenotypic adaptation to pesticides, anthropogenic events and/or changing climate.

Geographically, emerging and primary vectors exist in sympatry, allowing both ecological and genomic interactions.²⁰ Genomic interactions between these species may include adaptive introgression of alleles that have important consequences, such as resistance to insecticides, immune response and/or vector competence.^21,22 Notably, change in mosquito behavior and ecological interactions can also be driven by chromosomal inversion for example the 2La inversion in An. gambiae and An. coluzzii that has been associated with change in resting preferences, desiccation tolerance, thermal resistance and resistance to insecticides.^23,24 These vector adaptations can impact malaria transmission dynamics and erode successes towards malaria elimination, underscoring the need to characterize genetic diversity and evolutionary mechanism in all molecular forms and cryptic taxa of An. gambiae s.l. across local populations to inform the design of effective and targeted vector control measures.

In the Gambia, primary analysis of the genomes of the An. gambiae Bissau molecular form (hereafter Bissau) identified mutations associated with metabolic resistance.^25,26 However, an exhaustive genomic scan of signatures of local adaptation in Bissau in relation to sibling species has not been conducted. As the taxa could be identified across different ecoregions of the Gambia, different populations may have taken parallel evolutionary trajectories. Moreover, it is unknown if there is adaptive gene flow between Bissau and other sibling taxa, important for designing and implementing vector control measures in The Gambia.

Our study seeks to unravel the genetic architecture underlying local adaptation in Bissau mosquitoes from The Gambia using whole-genome variation data from the An. gambiae genome project (AG1000G) Phase 3(Ag3). By examining patterns of selection and genetic exchange among taxa, we aim to identify selective sweeps while looking for evidence of gene flow across taxa.

Materials and Methods

Data Processing and Code Availability

The genomic data used in this analysis were obtained from the Anopheles gambiae genome project (AG1000G) phase 3(Ag3). Analysis was restricted to 2658 mosquito samples identified as An. gambiae s.l, collected between 2011 and 2021, across different sites in The Gambia (Table 1). To reduce background noise and ensure statistical robustness^27,28 in haplotype frequencies, selection scan analysis was performed for sites with a minimum sample size set at 10 or a subset of 100 maximum at 100 for larger sample size, thereby balancing the census size for comparability.

Table 1 Sampling Site and Number of An. gambiae Complex Mosquitoes Collected in the Gambia

Signals of Recent Selection

To detect selective sweeps, three complementary haplotype-based statistical methods, H12,²⁹ integrated haplotype score (iHS)³⁰ and Cross-Population, Extended Haplotype Homozygosity (XP-EHH)³¹ were used.³² Each chromosome arm, ie, 2 (R&L), 3 (R&L) and X were assessed separately for signal of selection using H12. Non-overlapping windows were computed and calibrated for each taxon and cohort to account for variation in demographic history.²² Thereafter, regions with significant peaks of selective sweep were identified for each taxa and population cohort. Since the sensitivity of H12 to selective signals depends on population demographic history, which can potentially lead to false positives when detecting true sweeps, a window size was calibrated for each cohort rather than set globally. This calibration was done independently for each cohort using calibrate_h12_window_size function implemented in the malariagen_data package.³³ The optimal window size was identified by plotting the distribution of H12 value across a range of candidate window sizes and selecting the value at which the 95^th percentile of the H12 distribution was at or below 0.1. Following H12 computation, true signals were identified as genomic regions containing a maximum H12 value that was an outlier relative to the genome-wide distribution, accompanied by a characteristic peak architecture with roughly exponential decay on both flanks. Peak boundaries were then defined as the genomic interval over which H12 values remained continuously elevated above background on each side of the maximum. To identify selective sweeps between taxa, XP-EHH (function xpehh_gwss) from the malariagen_data python package³³ was used. Additionally, the iHS statistic derived from extended haplotype homozygosity (EHH) was used to identify recent positive selection within population cohorts.³¹ This statistic detects incomplete selection sweeps where alleles have not reached fixation.^30,34 We employed the function ihs_gwss from the malarigen_data python package.³³ Candidate SNPs under selection were defined as those exceeding the 99th percentile of genome-wide iHS scores and XP-EHH, corresponding to the top 1% of outlier loci. To identify genome regions with unusual levels of haplotype sharing between population cohorts, the H1X statistic was computed using the function H1X_gwss.³⁵ H1X detects haplotype sharing occurring at high frequencies between population cohorts, unlike traditional haplotype homozygosity statistics (ie, H1 and H12) that analyse patterns within a single population.³⁵

Haplotype Analysis

To explore population structure and haplotype similarity, we computed pairwise haplotype distances within candidate sweep regions. Sub-sampling was done to reduce biases arising from uneven geographic representation. For An. coluzzii, up to 20 individuals were randomly selected per administrative region (admin1_iso), retaining all samples where fewer than 20 were available. All An. gambiae samples were retained. For Bisaau taxon, 50 individuals were randomly sampled from the heavily represented GM-N region, while samples from other regions (excluding GM-U) were retained. Random sampling was performed without replacement using a fixed seed for reproducibility. Haplotype distances were computed across concatenated sweep regions using the haplotype_pairwise_distances function in the malariagen_data Application Programming Interface (API). Hierarchical clustering (Ward’s method) was used to construct a condensed distance matrix and the results visualised as clustered heatmap. Furthermore, diplotype cluster analysis was performed in region 3R:31,000,000–31,080,000 coding for Odorant Receptors (OR) genes. Here, clustering of diplotypes was achieved with complete-linkage hierarchical clustering.³⁶

Fixation Index

Genetic divergence between different cohorts of An. gambiae was assessed using F_st following the approach of Hudson et al³⁷ F_st is a statistic that measures genetic differentiation based on allele frequency differences between populations.³⁸

Genetic Structure at Signature of Selection

We used Principal Coordinate Analysis (PCoA) and Uniform Manifold Approximation and Projection (UMAP) to analyze genetic structure. Analyses were restricted to genomic regions under sweep selection following H12 analysis. Biallelic SNPs with a minor allele frequency (MAF) below 0.05 were retained and used to compute pairwise genetic distances between individuals (allel.pairwise_distance function from the scikit-allel Python library). The resulting city-block distance matrix was used to compute PCoA (function PCoA in scikit-bio package). UMAP (umap-learn Python package), for non-linear dimensionality reduction, was computed using the filtered biallelic sites and configured as follows: n_neighbors=15, representing the number of neighbors, balancing local and global structure, min_dist=0.5, signifying the desired separation between close points in the embedding space while preventing excessive cluster compression n_components=2, indicating the target embedding dimension for visualization and metric=“euclidean” specifying the metric used to compute the distance.^39–41 This allowed fine-scale visualisation of population structure based on regions under selection. This is consistent with previously published applications of UMAP in population genomics.⁴²

Neighbor Joining Tree

Evolutionary relationship within the An. gambiae complex was investigated using an unrooted Neighbor-Joining tree based on genomic regions under selection. The tree was built using city block (Manhattan) distance metric and visualised using function plot_njt from malariagen_data Python API.

Pairwise Ancestral Relatedness Within and Between Populations

We used the Refined Identity-by-Descent (IBD) program from Beagle software suite to examine shared ancestry and detect segments of IBD. The Refined IBD program uses probability to detect IBD segments, both short and long, based on genotype data and Hidden Markov Models.⁴³ Briefly, 100,000 biallelic SNPs from each arm of chromosome 3 were selected and pruned for linkage disequilibrium (LD) using PLINK2. Genotype data were phased using Eagle v2.4.1, a fast and accurate phasing algorithm suited for large-scale genomic data. Refined IBD was used to detect IBD segments at a minimum length of 0.125 centiMorgans (cM) and logarithm of odds (LOD) threshold of 3. This filtering ensured accurate detection of recent shared ancestry across individuals in the dataset.

Results

Overall, our finding revealed contrasting patterns of sweeps among taxa, with loci hosting chemosensory genes being under strong selection in Bissau population compared to sibling taxa. Additionally, we observed weak genetic differentiation among taxa with gene flow occurring between Bissau and An. gambiae s.s. compared to other sibling pairs.

Genetic Relatedness Among Gambian Anopheles gambiae s.l. Populations

Identity-by-Descent (IBD) analysis on chromosome 3R revealed strong relatedness between Bissau and An. gambiae s.s. (supplementary Figure S1a and b) and partial affinity with An. coluzzii, while An. arabiensis and An. melas formed distinct genetic isolation (supplementary Figure s1c). In addition, PCoA revealed that An. arabiensis and An. melas formed distinct clusters separating from the remaining sibling species (supplementary Figure S2).

Common and Taxon-Specific Selective Signatures in Gambian Anopheles gambiae s.l. Populations

Genome-wide selection scans indicated insecticides as a major force behind adaptive evolutionary changes in the An. gambiae complex population: Selective sweeps were concentrated in regions harbouring metabolic resistance genes, including CYP450 monooxygenase (ie Cyp6aa/Cyp6p, Cyp9k) and glutathione S-Transferase (Gste) (Figure 1A). Strong signals were detected around the Voltage-gated sodium channel (Vgsc) locus on chromosome arm 2L in An. gambiae s.s (peak at 0.56) and An. coluzzii (peak at 0.28) (2–3 kb, Figure 1A). Downstream of Vgsc was a weak sweep (peak at 0.15) hosting carboxylesterase gene (Coeae2f) exclusive to An. gambiae s.s. In An. coluzzii, two prominent selective signatures (peak at respectively 0.37 and 0.31) hosting fumarylacetoacetate and E3 ubiquitin protein ligase genes were detected (26–30 kb, Figure 1A).

Figure 1 Genome-wide selection scan (GWSS) using H12 statistic (A) Bissau vs An. gambiae and An. coluzzii (B) Chromosome arm 3R shows locus hosting ORs genes as absent in sibling species An. arabiensis and An. melas. Genes of interest are highlighted with red, y-axis indicates H12 score per chromosome for different cohorts (species), x-axis indicates the chromosomal arm (C) Manhattan plot showing signature of selection as detected by XP-EHH across chromosome arm 3R which host the OR and Gste genes. The XP-EHH score represents that of Bissau in comparison to An. gambiae s.s. (D) Manhattan plot showing pairwise XP-EHH score comparison between Bissau and An. coluzzii. Y-axis indicates XP-EHH values: X-axis indicates chromosome regions. Candidate genes are detected in the top 1% of XP-EHH. Negative scores signify selection in Bissau while positive scores signify selection in sibling taxa.

On chromosome arm 2R, selective sweeps were centred on Cyp6aa/Cyp6p gene (shared across all taxa), with peack≈0.49 for An. gambiae s.s and Bissau and weak (< 0.1) for An. coluzzii, and alpha N-acetyltransferase gene restricted to An. coluzzii at 0.23 peak (Figure 1A). On chromosome arm X, the signal for Cyp9k1 (15 kb) overlapped with a cuticular protein gene (CPR125) (0.34, 0.57 respectively for An. gambiae s.s and An. coluzzii), suggesting combined metabolic and cuticular resistance mechanisms (Figure 1A). Whilst the selective sweeps were generally weak on chromosome arm 3L across all taxa (peak < 0.15) (Figure 1A), a clear signal at Gste (28 kb) was evident on chromosome arm 3R in An. gambiae s.s. (peak 0.26) and Bissau (peak 0.14), and the odorant receptor (OR) signal further downstream, restricted to Bissau only (peak 0.22) (Figure 1A). Interestingly, extending H12 analysis of sibling species (An. arabiensis and An. melas) failed to detect ORs sweep signal, though Gste and CYP450 signals were evident (Figure 1B) from further sweep selection analysis. Consequently, these two species were excluded.

Pairwise XP-EHH scans of chromosome arm 3R revealed distinct selective pressure driving genetic divergence amongst taxa. A strong signal of selection was consistently detected in the Gste genes cluster (28.5kb) in both An. gambiae s.s (XP-EHH score >2) and An. coluzzii (XP-EHH score < 2), indicative of recent insecticide-driven selective sweeps (Figure 1C and D). However, in Bissau population, the selection was most evident at loci associated with sensory perception, specifically olfactory receptor regions (ORs, 32 kb; OR29, 25.74 kb), suggesting species-specific local adaptation which could be related to environment or host-seeking behaviour (Figure 1C and D).

iHS Analysis Identified Candidate Genes Under Selection in Insecticide Detoxification, Immune and Sensory Reception

The iHS analysis performed on An. gambiae complex genome identified 159 genes under recent selection (see supplementary Table S1). These genes spanned multiple functional pathways, including sensory reception, immune response and detoxification of xenobiotics (insecticides), consistent with ongoing local adaptation in a heterogeneous environment (Figure 2). SNPs associated with insecticide resistance showed consistently elevated iHS scores (>5), highlighting persistent selective pressure from either vector control interventions or agricultural pesticides (Figure 2).

Figure 2 Chromosomal distribution maps and integrated haplotype scores (iHS) showing loci under selection in An. gambiae sibling species defined by specific colours. Each panel shows a chromosomal arm (labelled top left of each panel). The rows in each panel represent the karyotype (top row), iHS values for each SNP plotted on y-axis against the SNP position on the x-axis). The bottom row marks the regions with elevated iHS, with the gene labels shown on the kayotype in the top row.

Selection Signatures of Genes Involved in Insecticide Metabolic and Cuticular Resistance Pathway

Genome scan of An. gambiae species genomes with iHS, much like H12, did not detect selection at the Vgsc gene in chromosome arm 2L within the Bissau population (Figure 2 and supplementary Table S1). Nonetheless, strong selection signals were evident across multiple metabolic resistance genes, particularly the CYP450 family on multiple chromosomal arms. For example, in An. coluzzii, selection targeted Cyp4 genes on chromosome arms 2L (Cyp4J5, Cyp4J9, Cyp4J10) and 2R (Cyp4K2, Cyp4AR1, Cyp4D15, Cyp4D17, Cyp4D16, Cyp4D22) whilst in Bissau, Cyp4C genes (Cyp4C35, Cyp4C36, Cyp4C27) were under selection on chromosome arm 3R. Similarly, Cyp6 family genes (Cyp6AA1-2, Cyp6P1-5, Cyp6P15, Cyp6AD1), known mediators of pyrethroid resistance, were consistently under selection in all sibling taxa on chromosome 2R. The Gste cluster genes (Gste1-8), conferring resistance to DDT and pyrethroids, were under selection on chromosome arm 3R of An. gambiae s.s and An. coluzzii (Figure 2 and supplementary Table S1). Beyond CYP450 and Gste, selection targeted additional unique detoxification loci. In Bissau population, Choline acetyltransferase gene (ChAT, 20 kb) on chromosome arm 2R, a rarely studied resistance candidate compared to well-established Ace-1 gene associated with resistance to carbamates and organophosphates, showed strong selection.

Carboxylesterases were also implicated: Coeae5O on 2R was under selection in both Bissau and An. coluzzii, while Coeae2F, associated with organophosphate resistance, was selected in An. gambiae s.s. (Figure 3 and supplementary Table S1). Strikingly, candidate regions extended beyond classical detoxification pathways. The genes linked to cuticle formation, including CPR16 (2L, 15 kb) and CPLC (CPLCA1–3, CPLCX3, 26–27 kb) in Bissau, CPR135 (2L, 28 kb) and CPR144 (2L, 29 kb) in An. coluzzii, and CPR143 (on 3L) in An. gambiae s.s., were also under selection (Figure 3 and supplementary Table S1). Additionally, ABC (ATP-binding cassette) transporter family genes, associated with xenobiotic detoxification, were enriched in An. coluzzii (2L), including ABCA1, ABCA2, ABCB4, and ABCC7 (Figure 2 and supplementary Table S1).

Figure 3 Allele frequencies and haplotype diversity across OR genes in Bissau. (A) Heat map of allele frequencies for non-synonymous mutations in OR53 and OR13. (B) Diplotype clustering at OR15 showing a unique allele (C>T: T336I) in Bissau. Hierarchical clustering dendrogram aligned with observed heterozygosity, displayed as a horizontal bar. The leaves of the dendrogram are colored by the species of the individual. Substitutions in each individual are presented below the dendrogram with blocks of the same color indicating shared mutation.

Selection Signatures of Genes Involved in Immune Response Pathway

Besides insecticide-driven adaptation, iHS analysis revealed strong selection on immune-related genes across An. gambiae complex, highlighting immunity as an additional driver of ongoing evolution. The chromosomal distribution of genes under selection of these genes also varied among taxa. In Bissau, the selection targeted G-protein-coupled receptor neuropeptide (GPRNP1) on chromosome arm 2R, CLIPB14 on chromosome arm 3L, and glutathione peroxidase 2 (GPXH2) on chromosome X. These loci (GPRNP1 and GPXH2) were under selection in other taxa (Figure 2 and supplementary Table S1). In An. coluzzii, the selection signatures included CLIP-domain serine protease on chromosome 2R, serine protease (SP24D) on chromosome 2L, serine protease inhibitor (SRPN4) on chromosome arm 3R, peptidoglycan recognition proteins (PGRPs 2–3), TEPs genes (TEP1, TEP3, TEP6, TEP10, TEP11 and TEP18), short caspase (CASPS9-13), and fibrinogen on chromosome 3L; (Figure 2 and supplementary Table S1). In An. gambiae s.s., the selection of immune-related genes was more widespread, involving the clusters of C-type lectins (CTLMA7 and CTL2), HPX genes (HPX1, HPX14, HPX15), CLIPA8, and TEP4 on chromosome arm 3L (Figure 2). Additional selection peaks on 2R included CLIPC2–3 and salivary gland proteins, while clusters of CLIPB genes (CLIPB11–12, CLIPB18), Mitogen-activated protein kinase, and serine protease inhibitors (SRPN5–6, SRPN16) were detected on chromosome arm 3R (28–29 kb) (Figure 2 and supplementary Table S1).

Selection Signatures of Genes Involved in Sensory Receptor Pathway

iHS analysis revealed selection signatures on genes associated with gustation, olfaction and host seeking, highlighting a behavioural adaptation to environmental chemical cues. In Bissau, the candidate loci selected included gustatory receptor 56 (Gr56) on chromosome arm 2L (26–27 kb), a cluster of olfactory receptor (ORs) genes (OR60-61, OR63-64) on chromosome arm 3L and multiple ORs and odorant-binding proteins (Obp43, Or13, Or15, Or17, Or30, Or46, Or53, Or55) on chromosome 3R (31–32 kb) (Figure 2 and supplementary Table S1). In An. coluzzii, sweep selection signals were detected at gustatory receptor 32 (GR32) and a cluster of odorant binding proteins (Obp71, Obp49-55 region 39 kb) and Obp 69 (29–31 kb) on chromosome arm 2L (Figure 2 and supplementary Table S1). In An. gambiae s.s, the chemoreceptor genes (Or26-27, Or56-57, Obp19) and ionotropic receptors (IR41t.1–2), were under selection (Figure 2 and supplementary Table S1). Across taxa, selection sweeps were also detected in glutamate receptors (GLURIIa, GLURIIb) on the X chromosome, with GLURIIc additionally under selection in An. gambiae s.s. and An. coluzzii (Figure 2 and supplementary Table S1). The divergent but overlapping patterns of selection on chemosensory pathways reflect their role in ecological adaptation and host-seeking behavior in Anopheles gambiae s.l.

Whilst selective sweeps spanned multiple chemosensory gene families (ie IR, OR, GR and OBP), OR genes emerged as the most divergent with Bissau populations carrying distinct, high-frequency variants that were rare or absent in sibling taxa (supplementary Figure S3). For example, amino acid changes were more abundant in ORs 53 and 13 in the Bissau population than in sibling taxa (Figure 3A). These were particularly distinct as most variants had allele frequencies >60 in Bissau, contrary to very low frequencies in sibling taxa, except for the S237R and I382V variant that were nearing fixation in all taxa (Figure 3A). These OR mutations occurring at high frequency have not been validated for association with any phenotype. Diplotype clustering revealed a region with low heterozygozity within the Bissau population, suggesting positive selection (Figure 3B). Besides, distinct patterns of diplotype clustering was observed for other ORs across sibling taxa suggesting divergent adaptation mechanisms (supplementary Figure S3).

Adaptive Gene Flow and Parallel Evolution of Insecticide Resistance Reveal Closer Relatedness Between Bissau and An. gambiae s.s.

Haplotype cluster analysis revealed asymmetric adaptive gene flow across An. gambiae complex with the strongest haplotype sharing occurring between An. gambiae s.s. and Bissau than any other pair of sibling taxa (Supplementary Figure S4). Genome-wide clustering of sweep regions identified three major haplotype groups, one of which was mainly shared between Bissau and An. gambiae s.s. (Supplementary Figure S4a). Consistently, pairwise F_st analysis showed Bissau to be more genetically differentiated from An. coluzzii (F_st 0.168) than An. gambiae s.s. (F_st 0.124) (supplementary Figure s4b). Additionally, UMAP clustering grouped other sibling taxa distinctly with Bissau and An. gambiae s.s apart from the other sibling taxa (Supplementary Figure S4c).

At insecticide metabolic resistance loci, adaptive sharing was also most evident between Bissau and An. gambiae s.s. Haplotype cluster analysis at Cyp9K1 (15kb, chromosome X) revealed a haplotype shared between the two taxa, alongside a haplotype common to all taxa (Supplementary Figure S4d). Similarly, at the CYP450 (chromosome 2R, 28.8 kb) and Gste (chromosome 3R, 28 kb) metabolic resistance clusters, Bissau shared haplotypes with An. gambiae s.s. but not with An. coluzzii (supplementary Figure S4e and Figure S5).

Within the Bissau population, haplotype cluster analysis of the olfactory receptor genes region (chromosome 3R, 32 kb) revealed an isolated haplotype cluster in the North Bank region (Figure 4A). While gene flow was evident amongst Bissau population cohorts (Figure 4A), low-level genetic differentiation was observed in the North Bank population (F_st, 0.001 and 0.012) relative to the other sampling locations (Figure 4B). PCoA and Neighbor-joining analysis confirmed this sub-clustering, identifying the North_Bank population as a relatively separate cluster within the Bissau cohorts (Figure 4C and D). Additionally, haplotype clustering that focused solely on geographical distribution, revealed Bissau to be largely distributed in the North (Figure 4E).

Figure 4 (A) Haplotype clustering within a region containing OR genes. Individuals are ordered by sampling regions and dark blue colour indicate haplotype that are shared between populations (B) Fst analysis showing GM-N (North-bank) cohorts of Bissau population displaying some level of genetic differentiation (C) PCoA showing low genetic variation between An. gambiae lineages (D) NJ tree showing relationships across all cohort and taxa (E) Clustering of haplotypes by geographical origin.

Discussion

We analyzed An. gambiae complex genome from The Gambia for signatures of selective sweep to gain insights into their adaptive responses to selection pressures. There were strong selective sweeps, mainly associated with resistance to insecticide, vector competence and mosquito behavioral adaptation. The selection signatures of resistance encompassed well-known mechanisms including target, cuticular and detoxification enzymes. Interestingly, the patterns of selection varied among species whereas spatial variation was limited. This somehow reflects vector control strategies in The Gambia where insecticides are applied across the country than varying across regions. In contrast, the observed variation in adaptation to insecticides among species could be attributed to behavioral adaptation where species display temporal divergence in breeding peak and therefore a different impact of insecticides.

Detection of selective sweeps associated with metabolic resistance is consistent with other previous studies. For example, we detected several Cyp6 family^8,44–46 and Cyp4 family (ie Cyp4AR1, Cyp4D15, Cyp4D17, and Cyp4J5) aligning with recent findings from the Sahel region.⁴⁷ Another consistent finding was the detection of Cyp9K1 gene, which was recently linked to type I and II pyrethroids resistance in mosquitoes from East and Central Africa.^48,49 These further highlight Cyp9K1 as a potential driver of insecticide resistance in An. gambiae s.l across the continent. Interestingly, haplotype clustering analysis supported introgression rather than independent local evolution as a potential factor underlying CYP9K1 haplotype evolution among sibling taxa. Other metabolic resistance genes under selection among An. gambiae s.l and reported in several malaria endemic countries across the African continent are Gste^18,49–52 and carboxylesterase (Coeae2F) genes.⁵³ As is for Cyp9k1, haplotype clustering constructed on genomic regions that host Gste revealed adaptive sharing between taxa pairs, illuminating the potential of metabolic resistance genes spreading across sibling taxa through adaptive sharing. We note, however, that haplotypes sharing between sibling taxa within An. gambiae s.l complex at loci associated with insecticide resistance may also reflect incomplete lineage sorting (ILS) or more plausible parallel evolution, given their close evolutionary relationship, calling for a more detailed analyses across species and the continent.

Another notable finding was the detection of selective sweeps associated with individual species, further suggesting species-specific response to selective pressure. For example, sweeps that have not previously been associated with insecticide resistance were evident for An. coluzzii. The unique sweeps were detected at genomic regions around alpha N-acetyltransferase (AGAP002284) and E3 ubiquitin-protein ligase. Both genes have been detected in other mosquito populations, although their role in insecticide resistance remains largely unexplored.^54,55 However, alpha N-acetyltransferase has been implicated in detoxification of xenobiotics including agricultural pollutants some of which share a similar chemical class with some insecticides used in vector control.⁵⁶ Currently, An. coluzzii is among the most geographically widespread malaria vectors in the Gambia and could be in contact with some pesticides (eg neonicotinoids) used in agriculture, which may be driving such selection. Some of these are also commonly incorporated as synergists to increase toxicity against pyrethroid-resistant mosquitoes.⁵⁷ Possibly, these novel mutations could be allowing An. coluzzii to be more adapted to insecticide pressure, enabling its spatial spread across different ecologies. Indeed, resistance against neonicotinoids insecticides was recently detected in An. gambiae population from Cameroon following phenotypic bioassay.⁵⁸ Nonetheless, the role of these novel loci under selection in An. coluzzii should be functionally confirmed.

The selective signature in Bissau around the Choline acetyltransferase gene (ChAT) on chromosome arm 2R was unexpected as it is rather the well-known acetylcholinesterase (Ace-1) resistance gene that was previously found under selection, a common occurrence in parts of West Africa.⁵³ Both transferase and esterase enzymes are involved in the metabolism of acetylcholine (Ach) and it is possible that the higher frequency of ChAT associated alleles is being driven by selection pressures from sustained insecticide application and a favorable fitness cost. To preserve the effectiveness of available insecticides, we suggest further monitoring of this novel variant alongside other conventional insecticide resistance markers. Phenotypic analysis to help identify the class of insecticides driving selection at this locus is also important.

The genomic region around the ATP-Binding Cassette (ABC) genes being under selection was another noticeable finding in An. coluzzii. Although the connection between ABC transporters family and insecticides resistance is not well established for mosquitoes, they have been linked to pyrethroid resistance in some agricultural pests including Plutella xylostella.^59–61 Recently, ABCs were found to be highly expressed in a deltamenthrin-resistant strain of Anopheles sinensis, a primary malaria vector in China and southeast Asia.⁶² One of the mechanisms by which ABC transporters are thought to modulate resistance against deltamenthrin is through alteration of cuticle thickness, thus reducing insecticide penetration.⁶³ We detected selective sweeps spanning several cuticular proteins, including CPR125, commonly associated with rigid cuticle, potentially lending credence to cuticle thickness playing a role in enhancing insecticide resistance in An. gambiae s.l. Cuticle thickening mechanism of insecticide resistance and its confluence with immune response was further supported by evidence of selection of CLIP gene sub-family, an immune response gene previously linked to resistance phenotype.⁶⁴

For An. gambiae s.s., currently decreasing in The Gambia, there were multiple strong sweeps under selection compared to sibling taxa. Concerningly, the sweeps under selection code for genes linked to both target site (Vgsc) and metabolic insecticide resistance (Cyp6aa/Cyp6p, Cyp9k1 and Gste). This could be in response to intensive insecticidal interventions against malaria over the last couple of decades. If these variants are selected directionally to fixation, they may enable resistance against different classes of insecticides. The outcome could be An. gambiae s.s expansion and re-establishment across The Gambia, thereby undermining malaria elimination efforts.⁶⁵

Besides resistance, our sweep analysis detected immune-related alleles that are known to shape mosquito-plasmodium interaction in An. gambiae s.s. An example is HPX15 gene commonly associated with defence against pathogens and long-term fertility.⁶⁶ Additionally, under selection in An. gambiae s.s. was the Mitogen-activated protein kinase (AGAP029511) gene that plays a key role in P. falciparum male gametogenesis, also exflagellation, an important step for parasite development in mosquito mid-gut.⁶⁷ Several genes of the TEPs family, key components of the innate immune systems of mosquitoes, were under selection including TEP-1 that is commonly associated with malaria parasite clearance and works alongside CLIP serine protease.^65,67,68

For the recently described Bissau population, distinct and shared selective signatures were also revealed. For example, like An. coluzzii it had ChAT under selection but lacked signals at regions coding for target site genes, corroborating previous findings by Caputo et al²⁵ Comparable observations have been reported for Anopheles funestus where metabolic genes and not target site are frequently implicated in insecticide resistance.^69–71 Broadly, target site resistance is expected to succeed metabolic resistance as non-insecticide xenobiotics can exert selection pressure on metabolic genes well before mosquitoes encounter insecticides used for malaria control. An alternative explanation is that mutations conferring metabolic resistance are generally suggested to have lower fitness cost than target site alterations, thus favoring their fixation. However, the resistant-associated haplotypes of Vgsc genes can still spread to Bissau through adaptive introgression as reflected in shared haplotype across the three taxa in chromosome arm 2L.

As seen for the primary malaria vectors (An. gambiae s.s. and An. coluzzii), the CLIP serine protease (CLIP14) was under selection in Bissau. CLIP14 is associated with melanization response against Plasmodium, and this finding could suggest Bissau is involved in malaria transmission. The ecological interaction of Bissau with the human host was further suggested by the selective sweep signals in the region coding for odorant receptors such as OR30, an indole-sensitive odorant receptor in some mosquito species and Drosophila melanogaster.⁷² Besides, in the common fly Musca domestica, both ORs 30 and 46 are categorized as indole-sensitive receptors.⁷³ Indole serves dual role, in some mosquito’s species being an oviposition cue while in others inhibiting host attraction.⁷⁴ The OR53 associated with host-seeking in mosquitoes was found to have multiple non-synonymous mutations occurring at high frequency, possibly driven by an ongoing change in the behavior of Bissau.^75,76 Similarly, OR13 had appreciable frequencies of nonsynonymous substitutions, suggesting that it is under a strong positive selection in Bissau compared to other sibling species. As the alleles were not shared between Bissau and other sibling species, the selection may be an independent evolutionary event in this taxon.

A possible reason for OR evolution in Bissau could be the need to adapt and establish an ecological niche. For Anopheles mosquitoes, alleles associated with sensory perception are generally less studied. Yet local adaptation involving sensory genes may modify sensory perception and subsequently behavior. In Drosophila melanogaster (Fruit fly) for example, the OR13 is linked to pheromone detection although it is unclear whether mosquitoes retain comparable functions.⁷⁷ Nonetheless, comprehensive understanding of OR genes in Anopheles mosquitoes is needed as they can serve as targets for vector control strategies, through the disruption of important biological traits such as oviposition or host attraction.

The OR sweep signatures were detected across spatially widespread Bissau cohorts. Whilst some level of geographical segregation was evident among Bissau cohorts, the similar OR sweep patterns across all sub clusters was surprising. A similar sweep pattern was observed for loci involved in metabolic resistance (cyp6a/p and Gste), where no geographic location emerged as a resistance hotspot in The Gambia, although the Gambian river is expected to impose a strong barrier to gene flow between its banks. This may be due to the small size of The Gambia, less than 50 Km at its widest point and around 480 km extending inland, and recent studies show mosquitoes to be capable of long-distance migration.⁷⁸ With reference to vector control, these findings imply that for small countries like The Gambia resistance management may include using the same insecticides over a wide geographical scale and that insecticide rotation should follow periodic/temporal rather than spatial approach.

Interestingly, we failed to detect selection at loci commonly associated with response to climate conditions such as heat shock proteins, implying that adaptation to insecticide resistance more than climate change could be shaping evolution of emerging species. Overall, these results illuminate further the genetic landscape of An. gambiae s.l in The Gambia and reveal the need to continue genomic surveillance to understand and address the changing landscape of insecticide resistance in the region. Nonetheless, phenotypic analysis to further substantiate the genomic analysis are needed. A possible limitation on the interpretation of our findings, especially on gene flow barriers, could be the sampling scheme that was based on administrative regions rather than ecotypes.

Conclusions

This study constitutes a comprehensive scan of patterns of sweep selection and gene flow encompassing insecticide and non-insecticide resistance genes. The key findings are that An. gambiae sibling species in The Gambia display varying insecticide resistance mechanisms, with some resistance genes being under selection in one taxon while undetectable in others. Besides, the study validates the presence of key insecticide resistance gene previously detected in this region and reveals the likelihood of emergence of novel insecticide resistance markers. Indeed, phenotypic association and experimental functional validation are needed to support the association of genetic variants in these signatures of selection in adaptive process such as insecticide resistance. Selection of genes associated with sensory perception in the novel taxon, Bissau, as well as those implicated in insecticide resistance, immune system and parasite development especially CLIP14 commonly associated with melanization response against Plasmodium, calls for in-depth bionomic investigations of this species and its role in parasite-vector interactions and malaria transmission dynamics in the far west Africa. This work represents a major step towards advancing our understanding of the genomic structure and evolution of emerging malaria vectors, with the potential of informing vector control efforts in The Gambia.