Extensive new Anopheles cryptic species involved in human malaria transmission in western Kenya

Overview of molecular determination of Anopheles species

Out of the 3556 Anopheles mosquitoes, 87.1% (3099/3556) were determined by species-specific PCRs or multiplex-PCRs and sequencing as major species An. gambiae sensu stricto (hereafter referred to as An. gambiae) (1440), An. arabiensis (718), and An. funestus sensu stricto (hereafter referred to as An. funestus) (941) in the five study sites (Fig. 1, Table 1, Supplementary Fig. S1). A subset of 21 randomly selected individuals from each major species identified by PCRs were confirmed by ITS2 sequencing based on similarity (> 98%) to the sequences of anopheline voucher species retrieved from NCBI GenBank database (Supplementary Fig. S2).

The remaining 457 collected anophelines (12.9%) were classified into 18 rare species groups based on ITS2 sequence homology. Except for two species groups (An. sp.18 and An. sp.19), the ITS2 sequences of all the species were identified as different species based on their similarity (> 98%) to the sequences of Anopheles voucher species retrieved from NCBI GenBank database (Supplementary Fig. S2). The ITS2 sequences of two species could not match with similarity > 98% threshold to reference anopheline sequences or known vector species in GenBank databases, suggesting the existence of novel cryptic species.

Pairwise comparison of ITS2 sequence similarities of the 21 Anopheles species indicated that except for one pair with 98.5% identity between An. gambiae and An. arabiensis, all pairs showed a similarity of 90% or less with confirmed species classifications (Supplementary Table S1). Phylogenetic tree analysis indicated that the 21 species belong to two different Subgenus (Subgenus Cellia Theobald and Subgenus Anopheles Meigen) in five species series groups, including Myzomyia, Neocellia, Pyretophorus, Cellia, and Myzorhynchus series (Fig. 2, Supplementary Table S2). The two new species An. sp.18 and An. sp.19 as well as An. sp.17 (a recently reported species¹³) belong to two different series groups, and An. sp.18 belongs to a different Subgenus (Subgenus Anopheles Meigen). The ITS2 sequence of An. sp.9 is homogenous with that of An. theileri (GenBank acc. JN994151) and An. sp. 9 BSL-2014 (GenBank acc. KJ522821)¹⁴. The ITS2 sequences obtained in the study are available in GenBank with accession numbers: MT408564-MT408584.

Comparison of morphological and molecular identifications

Of the 3556 mosquitoes tested, 3226 (90.7%) samples with available morphological data were used to evaluate the accuracy of morphological identification as compared to molecular identification. Among the 3226 mosquito samples, 2192 (67.9%) individuals were morphologically identified as An. gambiae s.l., 938 (29.1%) as An. funestus, 94 (2.9%) as An. coustani, and the remaining 2 (0.1%) as An. pharoensis (Table 2). The An. gambiae s.l. complex and An. funestus group (An. funestus, An. cf.rivulorum, and An. leesoni) had a similar percentage of matches (gambiae complex: 85.8%, 1881/2192, funestus group: 85.2%, 799/938) between molecular assay and morphological identification, while only 53.2% of specimens morphologically identified as An. coustani were confirmed by the molecular assay (50/94). Based on molecular assays, An. gambiae s.l. complex had the lowest misidentification (4.3%), followed by An. funestus group (6.8%), while 18.0% (11/61) An. coustani specimens were morphologically misidentified as An. gambiae complex (9) or An. funestus (2). Based on morphological identification, less than 15% of specimens morphologically assigned to An. gambiae complex (14.2%) or An. funestus group (14.8%) were identified as other species, while there were 44 specimens morphologically assigned to An. coustani (46.8%) that were classified by molecular assay into 9 anopheline species, including An. rufipes (17), An. funestus (7), and An. gambiae complex (6). Overall, more than 60% (264/427) of the rare species were morphologically misidentified as An. gambiae s.l. Specifically, nearly 20% (81/427) and 7.2% (31/427) of rare species were misidentified as An. funestus and An. coustani, respectively, whereas only 11.7% (50/427) of the rare species were correctly identified as An. coustani. Altogether, 84.0% (2710/3226) identification alignment was observed between the morphological and molecular analysis (Table 2).

Comparison of Anopheles species distributions and diversity

Overall, the three major species were found in all five study sites but at varying proportions. Anopheles funestus accounted for a large proportion (44.7–98.2%) of species observed throughout the five study sites (Fig. 1A, Table 1). Anopheles gambiae was the predominant species in two highland sites (56.7% in Emutete and 61.3% in Iguhu) and one lowland site (47.73% in Kombewa), whereas An. gambiae was nearly absent in Homa Bay (0.4%), a lowland site, and Kisii (0.6%), a highland site. An. arabiensis, was observed in high proportion in lowland areas (Homa Bay: 70.8% and Kombewa: 24.9%) than in highland areas, which ranged from 1.8% (Emutete) to 5.2% (Iguhu).

Seventeen of 18 rare species were identified in the highland areas, whereas only six rare species were detected in the lowland areas, suggesting that cryptic species might be more related to the sympatric An. gambiae than An. arabiensis. In lowland sites, the most abundant rare anopheline species was An. sp.15 (n = 17), followed by An. rufipes (n = 14) and An. cf.rivulorum (n = 14), whereas multiple rare species (such as An. christyi, An. sp.1, and An. sp.17) were identified in the highlands (Fig. 1B, Table 1).

A significantly higher species diversity was observed in the highland areas than in the lowland areas (Shannon index H, t-test, t =  − 6.59, df = 3419, p < 0.001). From 2015 to 2019, the average observed species richness (S) per year was 12.0 ± 0.27 in highland and 6.8 ± 0.25 in lowland (Table 3). In lowland, the highest species abundance (n = 10) was observed in 2018, whereas in highland, the highest species abundance of 15 was detected in 2016 and 2017. A significantly decreased diversity was observed from 2017 to 2018 (Shannon index H, t-test, t = 3.37, df = 340, p < 0.001) and 2018 to 2019 (Shannon index H, t-test, t = 3.83, df = 799, p < 0.001) in lowland. However, in highland, a significant increased diversity was observed from 2016 to 2017 (Shannon index H, t-test, t =  − 7.06, df = 833, p < 0.001), then decreased diversity from 2017 to 2018 (Shannon index H, t-test, t = 1.97, df = 787, p < 0.05), and remaining high species richness in 2019 (n = 10) (Table 3).

Molecular determination of Anopheles mosquito host blood meal source

A total of 1,372 blood-fed female mosquitoes from 16 Anopheles species were successfully genotyped for blood meal sources (Table 4). Of these, 41.6% females were identified to have had human blood meals, 53.6% of mosquito blood meals were identified as bovine, whereas the remained 4.8% individuals had blood meals originating from other animals, e.g., goat, pig, and dog. For the major vector species, the highest human blood index was found in An. funestus (0.72), followed by An. gambiae (0.51), whereas very few samples of An. arabiensis had human blood meals (Pearson’s Chi-squared test: χ² = 532.4, df = 8, p < 0.0001). The majority (> 90%) of An. arabiensis blood meals originated from cows (Fig. 3). Similar patterns of blood meal source were observed in the highlands and lowlands. Human blood meal sources were detected in a total of ten rare Anopheles species, including eight rare species from highland and two from lowland. The blood meal source of the other three rare species (An. leesoni, An. maculipalpis, and An. pretoriensis) were identified as bovine (Table 4).

Multiplex-qPCR identification of Plasmodium sporozoite infection

Of the total 3,556 female anopheline mosquitoes tested, 348 (9.8%) were positive for Plasmodium sporozoite infections (Table 5). Highest sporozoite rate was observed in An. funestus (17.3%; 95% CI 14.9–19.7%), followed by An. gambiae (9.9%; 95% CI 8.4–11.4%), and An. arabiensis (4.2%; 95% CI 3.6–4.8%). Sporozoite rates were two-fold higher in the lowland than that in highland sites for An. funestus (26.4% vs. 13.1%) and 1.5-fold higher in An. gambiae (14.3% vs. 7.8%). In contrast, An. arabiensis showed a higher sporozoite rate in the highlands (8.3%) compared to the lowlands (3.6%). Surprisingly, Plasmodium sporozoites were also detected in eight out of 18 rare anopheline species, confirming the role of rare cryptic species on malaria transmission in western Kenya (Table 5). Between 2015 and 2019, an increasing trend in the annual sporozoite rate was observed in the three major vector species both in the lowland (Chi-squared test: χ² = 9.04, df = 1, p < 0.01) and highland (Chi-squared test: χ² = 8.08, df = 1, p < 0.05) study sites (Fig. 4). From the outdoor collection of mosquito samples in 2016 and 2017, four rare species were found positive for Plasmodium infection in the two highland sites (Iguhu and Emutete) (Supplementary Table S4). No statistically significant difference in Plasmodium infection was detected between indoor and outdoor both in Iguhu (Pearson’s Chi-squared test: χ² = 0.78, df = 1, p = 0.37) and Emutete (Pearson’s Chi-squared test: χ² = 0.07, df = 1, p = 0.78). Three species of Plasmodium were detected with the majority being P. falciparum (87.6%, 305/348), and the remaining were P. malarae and P. ovale (Supplementary Table S5). Mixed species infections were common, accounted for 5.7% (20/348) of the total Plasmodium sporozoite positive samples. No statistically significant difference between major vector species in the component of parasite species was observed both in lowland (Pearson’s Chi-squared test: χ² = 7.35, df = 8, p = 0.49) and highland (Pearson’s Chi-squared test: χ² = 14.84, df = 12, p = 0.25).

Mitochondrial DNA barcode diversity of cryptic species

Eighteen of 21 COX1 sequence groups or species had a 1:1 relationship with the ITS2 sequence groups. Two species, An. sp.18 and An. sp.19 did not show reliable COX1 sequencing results due to failure of the PCR amplification using the universal COX1 barcoding PCR primers, indicating possible sequence variation in primer regions. However, high COX1 sequence variation was observed in the newly discovered species An. sp.17 (Supplementary Fig. S2). Among the 21 tested individuals, which showed identical ITS2 gene sequences, 13 haplotypes (An. sp.17_H01–An. sp.17_H13) of COX1 gene were identified with 17 separating sites. Higher haplotype diversity (Hd = 0.929, Pi = 0.0081) and nucleotide diversity were found in this cryptic species, as compared to An. gambiae (Hd = 0.700, Pi = 0.0029) and An. funestus (Hd = 0.710, Pi = 0.0015) (n = 21 for each species). Phylogenetic tree analysis indicated that at least four haplotype groups with three of them in different clades matching the reference sequences (An. sp.C, An. sp.D, and An. sp.F, respectively) from NCBI GenBank database (Fig. 5, Supplementary Table S3). No matching references in GenBank was found for the clade with haplotypes An. sp.17_H01, An. sp.17_H02, An. sp.17_H04, and An. sp.17_H06. The sequences obtained in the study were deposited in NCBI GenBank with the accession numbers: MT375202-MT375229.

Source: Ecology - nature.com