Genomic Footprints of Selective Sweeps from Metabolic Resistance to Pyrethroids in African Malaria Vectors Are Driven by Scale up of Insecticide-Based Vector Control

Abstract

Insecticide resistance in mosquito populations threatens recent successes in malaria prevention. Elucidating patterns of genetic structure in malaria vectors to predict the speed and direction of the spread of resistance is essential to get ahead of the 'resistance curve' and to avert a public health catastrophe. Here, applying a combination of microsatellite analysis, whole genome sequencing and targeted sequencing of a resistance locus, we elucidated the continent-wide population structure of a major African malaria vector, Anopheles funestus. We identified a major selective sweep in a genomic region controlling cytochrome P450-based metabolic resistance conferring high resistance to pyrethroids. This selective sweep occurred since 2002, likely as a direct consequence of scaled up vector control as revealed by whole genome and fine-scale sequencing of pre- and post-intervention populations. Fine-scaled analysis of the pyrethroid resistance locus revealed that a resistance-associated allele of the cytochrome P450 monooxygenase CYP6P9a has swept through southern Africa to near fixation, in contrast to high polymorphism levels before interventions, conferring high levels of pyrethroid resistance linked to control failure. Population structure analysis revealed a barrier to gene flow between southern Africa and other areas, which may prevent or slow the spread of the southern mechanism of pyrethroid resistance to other regions. By identifying a genetic signature of pyrethroid-based interventions, we have demonstrated the intense selective pressure that control interventions exert on mosquito populations. If this level of selection and spread of resistance continues unabated, our ability to control malaria with current interventions will be compromised.

Fig 1. Africa-wide population structure analysis.

A) Gene diversity of 11 microsatellites from throughout the genome, showing a loss of diversity at the telomeric end of chromosome 2R. B) Gene diversity of 8 microsatellites on chromosome 2R (including AFUB6, FunR and FunO) shows that the loss of diversity is restricted to AFUB6 and FunR, the microsatellites located near the pyrethroid resistance QTL rp1. C) Neighbor-Joining tree based on pairwise F_st among population samples, estimated using 9 neutrally evolving microsatellites from throughout the genome (excluding AFUB6 and FunR on chromosome 2R, which may be evolving under positive selection). D) Barplot of assignment probabilities of individual genotypes to two ancestral clusters (the most likely number estimated from the data) estimated using Bayesian population structure analysis under an admixture model. Each bar is an individual genotype for 9 neutrally evolving microsatellites from throughout the genome (excluding AFUB6 and FunR). In panel A, points from left to right represent the following microsatellites: FunQ (on chromosome X); AFUB6, FunR, FunO (on 2R); AFUB11, FunL, AFUB10 (on 2L); AFND7, AFND19 (on 3R); FunF and AFUB12 (on 3L). In panel B, points from left to right represent AFUB3, AFND40, AFUB6, FunR, AFND6, AFND30, AFND32 and FunO (all on 2R).

Fig 2. Africa-wide analysis of genetic diversity across the rp1 pyrethroid resistance locus.

A) Genetic diversity (pi) across the 120kb rp1 in Cameroon (orange), Malawi (blue), and Mozambique (green). B) Phylogeny of the ‘BAC25’ fragment (located 25kb along the BAC sequence and 9kb upstream of the CYP6P9a gene), which shows the most extreme difference in diversity between Cameroon (orange) and southern populations Malawi (blue) and Mozambique (green), which form a single clade. C) Phylogeny of the CYP6P9a gene sampled from throughout Africa, showing clear geographical divergence between southern Africa (Malawi and Mozambique; 100% bootstrap support), East Africa (Uganda; 87% bootstrap support) and West/Central Africa (Ghana, Benin, Cameroon; 78% bootstrap support). D) Haplotype network for Non-synonymous nucleotide variants in CYP6P9a. Light blue = haplotype shared between Malawi and Mozambique (MAL/MOZ); pink = haplotype shared between Benin, Cameroon and Ghana (BN/CAM/GH); teal = Benin (BN); red = Ghana (GH); orange = Cameroon (CAM); green = Mozambique (MOZ); blue = Malawi (MAL); purple = Uganda (UG). The size of the shapes is proportional to the frequency of the haplotype and numbers on each branch show the mutational steps separating haplotypes.

Fig 3. The selective sweep at rp1 in southern Africa coincides with the scale-up in pyrethroid use in malaria control.

A) Gene diversity at microsatellites on chromosome 2R in mosquitoes collected before widespread pyrethroid-based intervention in Malawi, 2002 (purple) and Mozambique, 2002 (red), and ‘post-intervention’ in Malawi, 2010 (green) and Mozambique, 2009 (yellow). B) Fine-scale nucleotide sequence analysis of the 120kb rp1 locus in Cameroon (orange) or pre-intervention samples (Malawi (CKW2002) = brown and Mozambique (MOZ2009) = pink) compared to Post intervention samples from Chikwawa (green), Salima (light blue), Nkhotakota (blue) Malawi, Chokwe (purple) Mozambique, and Kaoma (dark blue) Zambia. C) TCS haplotype network at ‘BAC25’. Haplotypes including sequences from more than one location are denoted by an H and include the frequency, haplotypes from one location Malawi (blue) and Mozambique (yellow) are only present in pre-intervention samples. D) Change in average pairwise diversity (Pi) in the CYP6P9a gene pre- vs. post-intervention in Malawi and Mozambique. E) Phylogeny of CYP6P9a sequences including both Malawi (Ckw) and Mozambique (Moz) shows that post-intervention samples (red) are almost completely homogenous while pre-intervention samples (black) are diverse.

Fig 4. Genome-wide analysis of selection on the An. funestus genome.

A) Sites with a significantly different allele frequency between 2002 and 2014 per genomic scaffold plotted against scaffold length identifies scaffold KB669169 as an outlier in which significantly skewed sites are over-represented. B) Plot of P-values of the difference in allele frequency for sites on scaffold KB669169. Grey dotted lines indicate the 120kb region represented by a sequenced BAC clone containing the rp1 locus. C) Mean frequency of non-reference alleles (for 101 adjacent variant sites) on scaffold KB669169. D) Allele frequencies as in C, showing a magnified region spanning the rp1 locus. E) Mean frequency of non-reference alleles (for 51 adjacent variant sites) on the BAC clone containing the rp1 locus. Positions of genes in the P450 cluster are indicated and CYP6P9a and CYP6P9b are labelled.