Genomic diversity of the African malaria vector Anopheles funestus

Abstract

Anopheles funestus s.s. is a formidable human malaria vector across sub-Saharan Africa. To understand how the species is evolving, especially in response to malaria vector control, we sequenced 656 modern specimens (collected 2014-2018) and 45 historic specimens (collected 1927-1967) from 16 African countries. We find high levels of genetic variation with clear and stable continental patterns. Six segregating inversions might be involved in adaptation of local ecotypes. Strong recent signals of selection centred on canonical insecticide resistance genes are shared by multiple populations. A promising gene drive target in An. gambiae is highly conserved in An. funestus. This work represents a significant advance in our understanding of the genetic diversity and population structure of An. funestus and will enable smarter targeted malaria control.

Fig. 1.. Population structure of 656 Anopheles funestus specimens collected across Africa.

(a) Map showing sequenced individuals from sub-Saharan Africa grouped into 17 geographic cohorts based on their collection location (filled circular and triangular markers), and six PCA cohorts consisting of one or more geographic cohorts (open elliptical markers), with the grey shaded area showing the An. funestus range across the continent (adapted from Sinka et al. (1)). (*) One sample collected in GH-N clusters with Equatorial on the PCA and has therefore been removed in analyses using PCA-cohorts. (b) Number of heterozygous sites (in millions) per individual, with each cross representing one individual in subset_2 (females only, for subset_1 see fig. S2a) and ordered along the horizontal axis by geographic cohort (colour, same order as the legend) and number of heterozygous sites. (c) Projection along the first two principal components (PCs) computed on chromosome arm 2L, with the percentage of explained variance in the axis labels, and ticks denoting the 0 position along each axis. The PCA cohorts are indicated with open elliptical markers. The inset PC plot was computed only on samples from the Equatorial PCA cohort. (d) Sliding window PCA along the genome, each line corresponds to one individual, individuals from GH-N are excluded. The genome is split in overlapping windows of approximately 1 Mbp and 200 kbp steps, a PCA is performed for each window, PC1 values are plotted along the y-axis, and the window centres in genomic coordinates are plotted along the x-axis. (The corresponding plot displaying the PC2 values for the same windows is shown in fig. S5a). Inversion regions are indicated by horizontal grey bars, centromeres by vertical grey bars, lined ellipses indicate regions of divergence across several cohorts (Fig. 3, figs. S6–8), while the dotted ellipse showcases divergence specific to the South Benin PCA cohort (fig. S3c,d).

Fig. 2.. Segregating inversions on chromosome 3.

Sliding window PCAs (left panels) were computed on all samples except North Ghana. For visualisation purposes, only a subset of samples is shown in the sliding window PCAs. (a) (left) Sliding window PCA showing CM (red) and MW (blue) samples on the 3Ra inversion region; window centres in Mbp along the x-axis, PC1 values for each individual along the y-axis. The different karyotypes are visible as three different horizontal trajectories, from top to bottom 3R+/+, 3Ra/+, and 3Ra/a. All three karyotypes are present in both cohorts, and because PC1 captures a signal that is a combination of the inversion and geographic population structure, the upper and middle trajectories (3R+/+ and 3Ra/+, respectively) do not overlap for the different cohorts, see fig. S5d. (middle) Projection along the first two principal components computed on the entire inversion region using all samples; samples are coloured by inversion karyotype. (right) Map of karyotype frequencies per geographic cohort. (b) (left) Sliding window PCA on the 3Rb inversion region displaying GA (orange) and MZ-M (cyan). Several samples change trajectory for part of the inversion region (indicated by thicker, darker lines); these samples appear to be double recombinants that locally exhibit a different karyotype than they have for the rest of the inversion (see fig. S5e, Supplementary Information). (middle and right) as in (a). (c) Sliding window PCA displaying MW on the 3La inversion region shows a strong separation of karyotypes near the breakpoints, but the signal decays towards the inversion centre (see fig. S5b). (middle and right) as in (a).

Fig. 3.. Selection sweeps and known insecticide resistance variants in Gste2 and Gaba.

(a) Genome-wide H12 scans for signals of recent selection. H12 values range from 0 to 1, with higher values indicating excessive haplotype sharing, which is a signature of recent selection. North Ghana was excluded due to elevated runs of homozygosity, which adds noise and makes H12 values unreliable for detecting genuine selection signals. The y-axis runs from 0 to 1 for each cohort, the x-axis shows positions along the genome. Peaks of H12 values ≥0.4 are highlighted with a grey vertical bar; for peaks where only one cohort reaches this threshold, the bar is restricted to this cohort. (b) Sliding window PCA computed on modern samples and historic samples combined, excluding North Ghana (fig. S9d). For visualisation purposes, only present-day Equatorial and South Benin cohorts and the historic Equatorial cohort are shown. Here displaying a 20 Mbp region around the Gste2 gene. Present-day samples show a peak in the Gste2 region, which is likely a signal of selection, while historic samples (shown in black) do not exhibit this signal. (c) Sliding window PCA computed and samples subsetted as in (b), here showing a 20 Mbp region around the Gaba gene. Several historic samples follow the peak seen in present-day samples at the Gaba locus. (d) Haplotype clustering within a region containing seven Gste genes. The dendrogram is obtained by hierarchical clustering of phased haplotypes, and used to define haplotype clusters as groups of haplotypes with SNP divergence below 0.0005 (cutoff indicated as dashed horizontal line on the dendrogram). The first bar below the dendrogram shows the population of origin for each haplotype, the second bar shows the genotype for the known Gste2 L119F mutation (note that this mutation was filtered out before haplotype phasing, so each haplotype is coloured by the genotype of the individual it belongs to; fig. S6a, table S3). (e) Haplotype clustering within the Gaba gene. Same structure as in panel (d) with the red bar showing the genotype for the rdl A296S variant (fig. S6b, table S3). (f) Maps showcasing the frequency of the L119F mutation in present-day (above) and historic (below) sample sets. (g) Maps showcasing the frequency of the rdl A296S mutation in present-day (above) and historic (below) sample sets.

Fig. 4.. Challenges for vector control.

(a) Cumulative number of gene drive targets (left) and number of genes containing gene drive targets (right) per geographic cohort as well as all subset_2 individuals (funestus) and An. gambiae and An. coluzzii individuals from Ag1000G phase 1 (16) (gambiae). The insets are zoomed out versions of the main plot, showing the results for the fully explored datasets; the areas of the main plots are indicated as dashed boxes in the insets. (b) Sequence variation in the dsx gene drive target on chromosome arm 2R. The dsx gene drive target is located at the boundary of intron 4 and exon 5 of the female-specific isoform. At the top, the seven exons of the female-specific isoform are depicted as boxes on the AGamP4 reference genome; coding sequence is shown in black, untranslated exonic regions (UTR) in grey. The target sequence is shown below, nucleotides in white boxes are located in intron 4, nucleotides in grey boxes are located in exon 5; the last three nucleotides constitute the PAM. The AGamP4 (An. gambiae) and AFunGA1 (An. funestus) reference genomes differ at the fourth base of the target sequence (highlighted in pink) in the intronic region. Variant 1 is found at very low frequencies in two cohorts from this study, Variant 2 is found at low to intermediate frequencies in four cohorts of An. gambiae and An. coluzzii and Variant 3 is found at very low frequency in one cohort of An. gambiae (see Supplementary Information). (d) Map showing fixation indices (F_ST) for neighbouring populations. F_ST is calculated on all accessible sites on the 2L chromosome arm. F_ST > 0.01 lines are dotted and become thicker the higher the F_ST. The background of the map is coloured by the Köppen-Geiger climate classification, adapted from Beck et al. (51).