The invasion of An. stephensi into the Horn of Africa is a significant threat to malaria control and elimination efforts in sub-Saharan Africa. Understanding the origins of this invasive species and ongoing gene flow can provide greater insight into its emergence in this region and improve predictions of future spread. Here we explore the genomic landscape of An. stephensi and demonstrate, using population structure analysis, that our Ethiopian field isolates are genetically distinct from Indian and Pakistani insectary colony populations. Ancestry analysis indicated that each field population constituted a distinct ancestral population. The Indian field isolates, despite the large geographical distance between the regions (Bengaluru and Mangaluru, around 350km), shared an ancestral population (K4) in the presence of several minor ones. The Indian colony samples were divided in two groups, based on their origin (Walter Reid or laboratory-reared strains). All samples from Ethiopia inherited the K1 ancestry, which was present at a minor proportion in some wild-caught Indian samples, and a single Pakistan SDA500 colony sample (intermediate form). This observation suggests a possible South Asian origin of the Ethiopian samples, as previously reported2. Previous studies using candidate gene analysis have identified high levels of genetic diversity across Ethiopian populations, implying either a one large invasion incident, or multiple smaller colonisation events31. The samples analysed here were collected in the same market town, which could have resulted from a population expansion in this area, resulting in the homogenous ancestry observed.
The ancestry observed in Ethiopia might have been dominant in other areas in South Asia at the time of the importation event. A similar relationship is observed between the Indian insectary-reared and field isolates, where the former was collected in 2016 in Chennai, India and have only K5 ancestry. This K5 ancestry also appears in 14 Indian field samples with the K1 Ethiopian ancestry, suggestive of shared ancestry between the Indian populations, despite the large geographical distance between Chennai, Bengaluru, and Mangaluru. Other studies have indicated that Anopheles spp. population structure remains stable over time, and physical distance is a larger driver of genetic variation36,37. To further understand the origins of the An. stephensi in the Horn of Africa, further WGS should be performed on more field isolates from An. stephensi’s native geographical regions, including India, Pakistan, Afghanistan, Iran, as well as elsewhere in the Horn of Africa. One of this studies key limitations is the over-abundance of colony samples analysed; the local field populations would likely exhibit far greater genetic diversity not represented in the samples used here.
The kdr L1014F insecticide resistance mutation was detected at a low frequency in this Ethiopian population and was absent in the Indian field samples. This kdr mutation confers resistance to pyrethroids and DDT which has been found previously in Indian populations of An. stephensi collected in 2016, and in Afghan populations, collected in 201824,25. The kdr L1014F mutation is absent in the Mangalore and Bangalore populations analysed here, despite evidence of extensive pyrethroid resistance in both these cities3,38. The low proportion of Ethiopian isolates with this SNP, along with its presence as a heterozygous genotype, implies it has recently arisen in this population. The kdr L1014F mutation has been previously identified in Ethiopian samples at a similar frequency to our study21. The low allelic frequency of the resistance marker, in the context of phenotypic resistance in the sampling location, indicates mechanisms other than target site modifications are responsible for this observed phenotype.
The other known target site mutation identified in the Ethiopian isolates was gaba A296S, which confers resistance to dieldrin20,22. This substitution was also identified in the Indian field populations, highlighting its prevalence, despite the insecticide being banned since the 1990s. The presence of this mutation corroborates the historical shared ancestry of the isolates analysed here. The gaba V327I mutation was also identified, and has a strong association to the A296S alteration, with five Indian field and colony samples carrying both SNPs39,40. A further 22 putatively novel SNPs were identified in genes associated with insecticide resistance. Six of these SNPs appeared in two Pakistani colony samples (Ste32 and Ste58) known to be insecticide susceptible. Two of these SNPs occurred in the ace-1 gene (N177D and V59A), and four in GSTe2 (F196Y, C146S, H97A, and G66A). The remaining 16 missense SNPs occur in populations where the resistance status is unknown, so it is not possible to infer any impact on insecticide susceptibility.
Of the genes identified as having ongoing directional selection, five could be potentially involved in insecticide resistance. The nAChR receptor subunit beta was found under selection by both within (iHS) and between (XP-EHH) population analysis. The nAChR receptor subunit alpha was also identified using the iHS metric as being under ongoing selection in Indian and Ethiopian populations. These two subunits of the receptor are the result of splicing of the nAChR gene; the presence (alpha) or absence (beta) of two cysteines determines their type41. Mutations within nAChR have been reported to result in resistance to neonicotinoids, which are pesticides that mediate synaptic transmission via nAChR, resulting in insect mortality42,43. This type of pesticide usage has been reported to result in neonicotinoid resistance in other Anopheles species44. Worryingly, neonicotinoids are considered an alternative to pyrethroids, for vectors with high levels of resistance to pyrethroids3,45–47. In Ethiopia, the President’s Malaria Initiative have using SumiShield, which contains a neonicotinoid since 2021. However, this is after the collection of these samples, so it is unknown whether they had been exposed to the chemical48. In addition, significant directional selection was detected in glutamate-gated channel genes in the Pakistan colony isolates. Mutations in these loci have been associated with ivermectin resistance in Drosophila49–51. Ivermectin is an anti-parasitic/endectocide, often used in mass drug administrations, which kills both the Plasmodium parasite, and Anopheles spp.mosquito when they ingest the blood of a treated host52–54. Ivermectin has been trialled as a vector control method using mass drug administration to help reduce malaria cases55.
Distinct signals of selection in a gaba gene were identified in India field and Pakistan colony samples, the insecticide fipronil that acts on gaba causing neurotoxicity56. This phenylpyrazole that has been suggested as a one health approach to vector control57. Similarly to ivermectin, fipronil can be used in mass drug administration, particularly in livestock (e.g., cattle targeting zoophilic vectors (e.g., An. stephensi), where this approach has been trialled successfully58–60. Resistance to fipronil has been reported in Iranian An. stephensi isolates, with both kdr and rdl mutations hypothesised to be involved61. Mutations within GABA receptors have been associated with reduced insecticide efficacy of fipronil, although this has not been observed in Anopheles spp62–64. Further surveillance of this gene could provide valuable insights into its potential involvement in fipronil resistance.
Another notable gene exhibiting strong directional selection by iHS was CYP307a1, identified in the Indian field populations. A 494bp deletion was detected in the intronic region of this locus in one Ethiopian isolate, and may impact on gene expression or result in alternative splicing65. CYP307a1 is a member of the cytochrome P450 gene family, and has previously been linked to resistance to both DDT and pyrethroids in An. funestus66,67. Similarly, in other insect species (Cydia pomonella), the upregulation of CYP307a1 has been associated with deltamethrin (pyrethroid) resistance68. Typically, CYP307a1 is involved in ecdysteroid hormone biosynthesis, these hormones control mosquito behaviour, nervous system development, and reproduction69. This gene has not been confirmed to directly result in insecticide resistance but warrants further investigation.
Other structural variants were detected that resulted in amino acid alterations, including a 67bp deletion in the CYP6a1 gene, found in a single Indian colony sample. This deletion was one of two detected in this gene, where the other was a 3’ UTR variant found in both Indian colony and field samples. Altered expression levels of CYP6a1 have previously been associated with deltamethrin resistance in D. melanogaster and C. pipiens70,71. With the absence of kdr mutations in these Indian field samples, but phenotypic pyrethroid resistance reported near the collection sites, it is likely other mechanisms contribute to resistance; such as metabolically mediated resistance3,19,27,72. This CYP6a1 gene identified here is not in the CYP6a cluster that was observed to have elevated coverage proportional to the genome median in Ethiopian isolates. To further investigate this gene cluster, read orientations and breakpoints would need to be identified to confirm whether this increased coverage was due to a copy number variation in the population.
Further investigations, using WGS or targeted amplicon sequencing, in tandem with susceptibility bioassays, are needed to investigate the impact of these mutations on insecticide response. The novel missense SNPs potentially linked with resistance should be used as targets in high-throughput molecular assays, to support surveillance and assist functional work to understand and validate underlying mechanisms associated with resistant phenotypes.
In conclusion, this study gives greater insight into the population genetics of cross-continental An. stephensi. Applications of WGS analysis to larger An. stephensi sample cohorts, across different geographical regions, will be key to understanding gene flow and identifying insecticide resistance markers. Such insights will enable public health authorities to make informed choices about vector surveillance and insecticide usage.