A population genomic unveiling of a new cryptic mosquito taxon within the malaria-transmitting Anopheles gambiae complex

Abstract

The Anopheles gambiae complex consists of multiple morphologically indistinguishable mosquito species including the most important vectors of the malaria parasite Plasmodium falciparum in sub-Saharan Africa. Nine cryptic species have been described so far within the complex. The ecological, immunological and reproductive differences among these species will critically impact population responses to disease control strategies and environmental changes. Here, we examine whole-genome sequencing data from a longitudinal study of putative A. coluzzii in western Burkina Faso. Surprisingly, many specimens are genetically divergent from A. coluzzii and all other Anopheles species and represent a new taxon, here designated Anopheles TENGRELA (AT). Population genetic analysis suggests that the cryptic GOUNDRY subgroup, previously collected as larvae in central Burkina Faso, represents an admixed population descended from both A. coluzzii and AT. AT harbours low nucleotide diversity except for the 2La inversion polymorphism which is maintained by overdominance. It shows numerous fixed differences with A. coluzzii concentrated in several regions reflecting selective sweeps, but the two taxa are identical at standard diagnostic loci used for taxon identification, and thus, AT may often go unnoticed. We present an amplicon-based genotyping assay for identifying AT which could be usefully applied to numerous existing samples. Misidentified cryptic taxa could seriously confound ongoing studies of Anopheles ecology and evolution in western Africa, including phenotypic and genotypic surveys of insecticide resistance. Reproductive barriers between cryptic species may also complicate novel vector control efforts, for example gene drives, and hinder predictions about evolutionary dynamics of Anopheles and Plasmodium.

Keywords: Anopheles; admixture; cryptic taxa; reproductive barrier; selective sweep; vector.

Figure 1.

Genetic distinctiveness of AT. (A) In a PCA plot with Tengrela, GOUNDRY, and Ag1000G individuals, AT occurs as a distinct cluster close to GOUNDRY. Variants were chosen based on segregation in our data and may show ascertainment bias affecting the relationships of individuals within Ag1000G; the salient result is how AT and Tengrela A. coluzzii relate to these other individuals (B) AT remains distinct in a PCA after combining Tengrela samples with GOUNDRY and the A. coluzzii samples collected alongside GOUNDRY. To control for differences between studies, all reads were trimmed to the same length, and then alignment and genotyping were performed jointly.

Figure 2.

(A) Most common phylogenetic topology among sections of AT genome and nominal species of the A. gambiae complex. Numbers at branches are not bootstraps, but the percentage of 100 kb windows that support each clade (above branches: entire genome; below branches: X chromosome only). AT is sister to A. coluzzii across 42.1 % of the genome (55.7% of the X chromosome), more often than to any other species, and 95.5% of the genome (100.0% of the X chromosome) supports a clade with AT, A. coluzzii and A. gambiae to the exclusion of the other species. (B) The mitochondrial DNA phylogeny shows that AT shares a single haplotype that occupies a unique branch close to the A. gambiae PEST reference genome. GOUNDRY samples occur near AT or near A. coluzzii and A. gambiae samples from Tengrela and elsewhere in Burkina Faso (BF), consistent with an admixed origin.

Figure 3:

Relationships between AT, GOUNDRY, and A. coluzzii autosomes using jointly called genotypes. (A) Analysis with ADMIXTURE suggests two ancestral populations, closely approximated by contemporary AT and A. coluzzii, with GOUNDRY showing ancestry from both. (B) Analysis with dadi corroborates this model, with an AT/A. coluzzii split over one million generations ago, followed by ongoing gene flow and a recent admixed origin of GOUNDRY. Population sizes (heights of colored bars) and migration rates (widths of arrows) vary across three time periods (demarcated with dotted lines). (C) Analysis with TreeMix shows GOUNDRY as sister to AT but with in-migration from A. coluzzii.

Figure 4.

Unique genomic characteristics of AT. (A) F_ST across the genome between AT and Tengrela A. coluzzii. Tens of thousands of variants distributed across the genome are highly divergent between these taxa (F_ST > 0.8; orange dots at top), while over a thousand sites, concentrated in several clusters, are fixed or nearly so (“definitive differences”, F_ST > 0.99; red dots at top). Average F_ST in 100 kb windows is more modest (blue lines), but three regions stand out representing the 2La inversion, TEP1, and CYP9K1. (B) Nucleotide diversity (π), intertaxon divergence (Dxy), and site frequency spectra (Tajima’s D) in AT and A. coluzzii. Nucleotide diversity is low in AT except at the 2La inversion. Tajima’s D is mostly positive in AT, but three regions show low Tajima’s D, low π, and many high-F_ST sites (as shown in A), suggesting selective sweeps. (C) In most AT females, read coverage along most of the reference Y chromosome substantially exceeds the X/autosomal average. Coverage is negligible around sex-determining gene YG2. A. coluzzii females, in contrast, typically show negligible coverage across the entire Y except for a few repetitive sections.

Figure 5.

Regions of low heterozygosity in AT. (A) Heterozygosity (H) in blocks of 1 Mb across the genome for all 51 AT individuals. There are many long homozygosity tracts (red) in most individuals, but these are consistent with the relatively low genetic diversity observed across the population, except in the 2La inversion. (B) Histogram of the heterozygosity blocks depicted in A, based on number of heterozygous sites. Heterozygosity shows a relatively smooth distribution with only a slight uptick for homozygous blocks (0 or 1 heterozygous sites), indicating little evidence for inbreeding driving homozygosity.

Figure 6.

Distinguishing AT from A. coluzzii using amplicon genotyping. A pool of five primer pairs were selected that amplify five polymorphisms diagnostic for AT, and jointly amplified PCR products were sequenced (10 AT, 10 A. coluzzii, and one negative control). For each marker, the vast majority of read pairs are consistent with the known taxonomic category as inferred from whole genome sequencing, allowing for unambiguous identification. A small number of read pairs were erroneously assigned to the negative control (“neg”), indicating that a conservative test should exclude any individual with unusually low coverage and/or intermediate frequencies of both alleles.