Surveillance site selection

Following the detection of An. stephensi in Kenya in December 2022¹², the Kenya National Malaria Control Program (NMCP) and partners conducted additional sampling targeted at areas bordering locations where An. stephensi was detected, urban or peri-urban areas neighboring confirmed An. stephensi presence, sites with high habitat suitability based on environmental conditions and population density, points of entry from neighboring countries in the North of Kenya, towns along major transportation routes with considerable movement of people, goods, and animals, and counties reporting unexplained increases in malaria cases, especially outside usual seasonal patterns and areas reporting increased cases of both P. falciparum and P. vivax malaria. Within each of the selected counties, sampling teams prioritized areas with high animal ownership, and the presence of water reservoirs. In rural areas: oases, major water reservoirs such as small dams, streams, and rivers along known cattle grazing routes and sites of urban development and construction.

Sampling for Anopheles stephensi

Sampling for An. stephensi targeted both larval and adult stages between January 2023 and June 2024 across Kenya, with a focus on both indoor and outdoor collections under the coordination of the NMCP. In each of the selected counties, three sub-counties and two towns/villages per sub-county were selected, except in Mandera, where mosquito surveys were carried out in one sub-county, because of insecurity. The team conducted adult and larval surveys for 8 days per county, simultaneously. In areas where An. stephensi was already established, either monthly or quarterly sampling was conducted if resources allowed, though most data reported here was from one time sampling efforts.

Larval Sampling focused on potential larval habitats including man-made water containers, freshwater pools and pans, stream margins, discarded tires and plastic containers, irrigation ditches, water storage containers (metal and plastic tanks, concrete cisterns, barrels, clay pots), construction sites and areas near animal shelters. Collected larvae were preserved in 70% ethanol for species identification and teams were instructed not to carry any live material outside of the areas of collection to prevent accidental introductions to new areas.

Quarterly cross-sectional surveys led by the NMCP involved sampling of adult mosquitoes using indoor CDC light traps, outdoor UV light traps, and mechanical aspiration indoor and outdoor using Prokopacks.

Evaluation of adult trapping methods in Marsabit County

A host of adult collection methods were evaluated in Marsabit county where the vector was thought to have been established. Methods that were evaluated indoor and outdoor included: UV light traps, human landing catches (HLC), CDC light traps (CDCLT), and mechanical aspiration using Prokopacks. Host decoy traps (HDT), double bed net traps (DBT), BG—Sentinel and BG Pro traps (Biogents AG, Regensburg, Germany) with BG lure (synthetic human odor attractant) were only evaluated outdoors. The methods were evaluated for the density of An. stephensi trapped.

Sample processing and identification

Collected samples were morphologically identified using established keys, focusing on specific features such as palp speckling, wing vein patterns, and thoracic characteristics¹⁴. All Anopheles larvae were included in the molecular analysis.

Molecular analyses

Molecular identification was conducted through an initial Colorimetric Loop-Mediated Isothermal Amplification Assay (CLASS) for preliminary detection as previously described¹⁵. Briefly, Single mosquito legs from whole mosquito samples and DNA extracted from larvae samples were analyzed for the presence of target DNA indicated by a visible color change in the reaction mixture, typically from pink to yellow. This rapid detection method was followed by a confirmatory species-specific PCR to validate the results. To confirm the species and conduct genetic analyses, a portion of the mitochondrial cytochrome oxidase subunit I (COI) locus from morphologically identified An. stephensi, DNA was PCR amplified and sequenced using previously published protocol⁵. Briefly, to amplify the partial COI locus, we used LCO1490F (5′-GGTCAACAAATCATAAAGATATTGG-3ʹ) and HCO2198R (5ʹ-TAAACTTCAGGGTGACCAAAAAATCA-3ʹ) primers¹⁶ with the following thermal cycling conditions: 94 °C for 5 min, followed by 35 cycles of 94 °C for 30 s, 56 °C for 45 s, 72 °C for 1 min, and a final extension of 72 °C for 10 min.

Amplicons were visualized on 2% agarose gel to confirm the correct locus was amplified, then cleaned using ExoSap (Cytiva, Marlborough, MA) and sequenced using Sanger technology with BigDye chemistry (EdgeBio, San Jose, CA) and run on an ABI 3730 Genetic Analyzer (Thermo Fisher, Santa Clara, CA). Sequences were cleaned using CodonCode version 11.0.2 (CodonCode Corporation, Centerville, MA, USA) and submitted as queries to the National Center for Biotechnology Information’s (NCBI) Basic Local Alignment Search Tool (BLAST) (Altschul et al., 1990) against the nucleotide database in GenBank under default parameters for highly similar hits (98–100%) to confirm the species.

Genetic diversity and population structure

To determine the genetic diversity, structure and evolutionary relationships of the Kenyan An. stephensi COI sequences generated in this study, population genetic and phylogenetic analyses were performed.

Multiple sequence alignment was performed using MAFFT version 7 [1] and uneven ends were trimmed using BioEdit version 5.0.9. Population genetic statistics were generated using COI sequences for each collection site in DnaSP version 6¹⁷. The statistics generated included the number of polymorphic (segregating) sites (s), number of haplotypes (h), haplotype diversity (Hd), and nucleotide diversity (π). For population structure characterization and further comparative analysis, we performed a phylogenetic analysis with the Kenyan An. stephensi COI sequences and An. stephensi COI sequences from both the native and invasive range retrieved from NCBI Genbank.

The global An. stephensi COI sequences included eight from India (Genbank accession number: KT899888¹⁸, KX467337, MH538704, MK726121, MN329060, LR736015, LR736014, and LR736013), four from Pakistan (Genbank accession number: KF406694, KF406693, KF406701, and KF406680)^18,19, one from United Arab Emirates (Genbank accession number: MK170098)²⁰, seven from Sri Lanka (Genbank accession number: MF975729, MF975728, MF124611, MF124610, MF124609, MF124608, and MG970565)²¹, two from Yemen (Genbank accession number: OM865140, and PP387838)^22,23, one from Djibouti (Genbank accession number: KF933378)², three from Sudan (Genbank accession number: MW197100, MW197099, and MW197101)²⁴, two from southern Ethiopia (Genbank accession number: OQ865406, and OQ865407²⁵, nine from eastern Ethiopia (Genbank accession number: OK663480, OM801691, OM801703, OM801693, OM801697, MH651000, OK663481, OK663479, and OK663482^5,26,27, two from Somaliland (Genbank accession number: ON421572, and ON421574)²⁸, and two from Kenya (Genbank accession number: OR607950, and OR607949)¹². Anopheles maculatus was designated as an outgroup to be consistent with previous An. stephensi phylogenetic analysis^29,30. Phylogenetic relationships were inferred using Mr Bayes version 3.2.7³¹ which is based on Bayesian inferencing and relies on calculating the posterior probability distribution of phylogenies. The general time reversible (GTR) nucleotide substitution model³² with GAMMA rates of heterogeneity was used to determine how nucleotides evolve. The trees were visualized in FigTree version 1.4.4³³ and the percentage posterior probabilities included. We also mapped COI haplotype proportions across the study sites. For consistency and comparison, haplotype numbering was kept the same as in Carter et al.²⁵ and color coded as in Samake et al.³⁴.

Predicting the probability of occurrence across Kenya

Extensive entomological monitoring can be expensive, particularly when the species in question is rare, and difficult to collect with existing monitoring tools. A country-specific species distribution modeling approach was used to identify hotspots for the next phase of extensive entomological monitoring to inform control efforts for An. stephensi spread in Kenya. An ensemble species distribution model using the package ‘biomod2’ version 4.2.5.2 in R³⁵ version 4.1.2 was employed, borrowing from the first prediction of An. stephensi invasion¹³. Climatic variable selection was based on a recent habitat suitability modeling approach used to identify entomological surveillance points that led to the first detection of An. stephensi in Ghana³⁶.

The environmental and climatic data were downloaded at a resolution of 1 km square (km²) for explanatory variable pre-processing in R version 4.1.2. There were 19 bioclimatic raster layers, and an elevation raster layer from WorldClim platform^37,38, and normalized difference vegetation index (NDVI) from 2020, 2021, 2022, 2023³⁸, alongside population density based on the 2019 population census in Kenya³⁹. In total, there were 25 variables. Pre-processing involved cropping to the boundary extent of Kenya, raster resampling to match their resolutions, and encoding to American Standard Code for Information Interchange (ASCII) format.

The 19 bioclimatic raster layers and elevation were assessed for correlation using the ENMTools package in R to reduce multicollinearity. A correlation coefficient cutoff of ≥ 0.7 was used to identify high associations^13,36. Where two rasters had a high correlation, only one was chosen for inclusion in the model. To minimize bias related to arbitrary variable selection and ensure the selection of the most unique variables, a process that used a count function to generate the correlation frequency for each variable across other variables was devised. This identified bioclimatic variables with a correlation coefficient ≥ 0.7 across more than 50% of the 19 bioclimatic variables and elevation (frequency ≥ 10/20 of variables). These were removed, followed by an additional stepwise elimination of other highly correlated variables.

Both An. stephensi presence and absence data were included as response variables. A binary coding was used to ensure readability by the models, with presences coded as “1” and absences as “0”. These were based on entomological monitoring conducted between December 2022 and June 2024. To reduce the overlap between presence and absence, absence records were not included in sites where a presence was recorded in subsequent sampling rounds. Actual absence records were selected in all counties where An. stephensi was not detected throughout the entomological monitoring efforts, to June 2024. In total, there were 4,128 unique geographic coordinates with absence, and 32 unique geographic coordinates with presence across the counties.

The ensemble modeling included four models, Random Forest, Generalized Additive Model (GAM), Gradient Boosted Machines (GBM), and Extreme Gradient Boosting (XGBoost), with the “Bigboss” calibration option for each model. The “Bigboss” options are a set of model calibrations optimized for species distribution modeling⁴⁰. An additional step using K-fold cross-validation, with 3 partitions and 10 draws of cross-validation data (a total of 30 cross-validation runs for each model) reduce overfitting by each model. Additionally, each model was specified for an 80:20 split for test and training across the 30 cross validation runs. Model evaluation was based on the area under the receiver operating curve (ROC) and True Skill Statistic (TSS). The performance of each model was evaluated for inclusion in the ensemble using TSS.