Resistance to pirimiphos-methyl in West African Anopheles is spreading via duplication and introgression of the Ace1 locus

Abstract

Vector population control using insecticides is a key element of current strategies to prevent malaria transmission in Africa. The introduction of effective insecticides, such as the organophosphate pirimiphos-methyl, is essential to overcome the recurrent emergence of resistance driven by the highly diverse Anopheles genomes. Here, we use a population genomic approach to investigate the basis of pirimiphos-methyl resistance in the major malaria vectors Anopheles gambiae and A. coluzzii. A combination of copy number variation and a single non-synonymous substitution in the acetylcholinesterase gene, Ace1, provides the key resistance diagnostic in an A. coluzzii population from Côte d'Ivoire that we used for sequence-based association mapping, with replication in other West African populations. The Ace1 substitution and duplications occur on a unique resistance haplotype that evolved in A. gambiae and introgressed into A. coluzzii, and is now common in West Africa primarily due to selection imposed by other organophosphate or carbamate insecticides. Our findings highlight the predictive value of this complex resistance haplotype for phenotypic resistance and clarify its evolutionary history, providing tools to for molecular surveillance of the current and future effectiveness of pirimiphos-methyl based interventions.

Fig 1. Ace1 mutations in African populations.

A) Frequency of non-synonymous SNPs in the Ace1 gene in African A. gambiae and A. coluzzii populations from Anopheles gambiae 1000 Genomes, Phase 2. For each SNP, we indicate peptide- and transcript-level coordinates and substitutions. B) Ace1 CNVs across African populations, including the frequency of specimens with >2 copies in each population. A diploid genome without duplications would have 2 copies. Populations where G280S and duplications are present are highlighted in bold text. Note: populations denoted with an asterisk (The Gambia, Guinea-Bissau and Kenya) have high frequency of hybridisation and/or unclear species identification.

Fig 2. Combinations of Ace1 G280S and CNV genotypes.

A) Fraction of reads supporting 280S alleles and number of Ace1 copies (1000 Genomes dataset, n = 1142). Boxes highlight groups of haplotypes: non-duplicated wt (purple), duplications with both wt and 280S alleles (orange) and duplications with only 280S alleles (green). B) Breakdown of observed genotypes at the Ace1 locus according to Ace1 and 280S copy number. The asterisk (*) denotes a wt-homozygous specimen from Guinea that carries an independently evolved Ace1 duplication. C) Diagrammatic summary of Ace1 haplotypes and their surrounding flanks; light grey rectangles show duplicated regions.

Fig 3. Genotype-phenotype association in Ivorian A. coluzzii.

A) Number of Ace1 copies compared to the estimated number of 280S copies in Ivorian A. coluzzii (n = 71), color-coded according to resistance phenotypes. Random jitter has been added for clarity. B) Cross-tabulation of pirimiphos-methyl resistance and three Ace1 mutations in Ivorian A. coluzzii: 280S allele presence, number of Ace1 copies, and number of 280S alleles. Orange lines denote ad hoc groups of genotypes where we identify changes in survival rates (included at the bottom of each table). Source data available in S6.

Fig 4. Principal component analysis of Côte d’Ivoire A. coluzzii.

PCA constructed using genotypes of 791 unlinked variants from chromosome 3.

Fig 5. Genome-wide scan of variants associated with pirimiphos-methyl resistance in Ivorian A. coluzzii.

A) Profile of population branching statistics along all chromosomal arms, calculated in consecutive blocks of 1000 segregating variants, using resistant and susceptible Ivorian A. coluzzii as populations A and B, and Angolan A. coluzzii as outgroup. Orange triangles indicate windows with extreme PBS values (p-values derived from a standardised distribution of PBS along each chromosomal arm, and FDR < 0.001), and the number of genes therein. B) Mapping coordinates of k-mers significantly associated with pirimiphos-methyl (n = 439). Most k-mers map to the Ace1 duplication region (n = 414) or, despite mapping elsewhere in the genome (n = 24), are correlated with Ace1 copy number (orange triangles). Only one k-mer mapping outside of the Ace1 duplication is not correlated with Ace1 copy number (purple triangle). C) Normalised frequency of each significant k-mer (n = 439, horizontal axis) in each genome (n = 71, vertical axis). k-mers are sorted according to their mapping location (in Ace1 or elsewhere), and genomes are sorted according to their phenotype (resistant/susceptible). D) Pearson’s correlation coefficients (r) between k-mer frequency and number of Ace1 copies in each genome (n = 439 significant k-mers). k-mers are coloured according to their mapping location (in Ace1 or elsewhere) and sorted by the values of r.

Fig 6. Haplotype clustering network around Ace1.

A) Minimum spanning tree of haplotypes around the 280S-linked variant 2R:3481632 (± 300 bp, n = 104 phased variants). Node size reflects number of haplotypes belonging to each cluster, which are color-coded according to species and geographical origin. Edges link haplotype clusters separated by one substitution. Singleton clusters have been removed from this view (see S10 Data). A cluster of identical haplotypes linked to the 280S-tagging variant is highlighted in green. B) Venn diagram representing the overlap of samples belonging to the 280S-linked haplotype clusters identified around each of the three tagging variants (S10 Data).

Fig 7. Positive selection around the Ace1 duplication.

A-C) Profile of Garud’s H₁₂, H₂/H₁ and haplotypic diversity around the Ace1 duplication, for haplotypes carrying the 280S- or wt-tagging variants at 2R:3481632. Statistics are calculated in blocks of 500 variants with 20% block overlap. Includes averages of each statistic outside the duplication, upstream and downstream of the breakpoints (with standard errors from jackknife haplotype resampling). D) Extended haplotype homozygosity (EHH) at the duplication breakpoints for haplotypes carrying the 280S- or wt-tagging variants at the 2R:3481632 locus. Additional plots for all tagging variants are available in S11 and S12 Data. Note that Garud’s H statistics and EHH assume a diploid genome, and estimates from within the duplication can thus be biased.

Fig 8. Introgression of the Ace1 duplication.

A) Profile of Patterson’s D statistic around the Ace1 duplication, testing evidence of introgression between duplicated A. coluzzii from Côte d’Ivoire (population A) and various A. gambiae populations (population C) with or without duplications (green and purple lines, respectively). We used non-duplicated Angolan A. coluzzii as a contrast (population B) and A. arabiensis as outgroup (O). D is calculated in windows of 5,000 variants with 20% overlap, and the average value of D along the duplication region is shown for the Côte d’Ivoire A. coluzzii / Ghana A. gambiae comparison, with standard errors derived from block-jackknife. B) Id., but using non-duplicated A. coluzzii from Côte d’Ivoire as population A. Detailed lists of all statistical tests and replicate analyses with additional populations are available in S13 Data. Species abbreviations: gam, A. gambiae; col, A. coluzzii; ara, A. arabiensis. Country abbreviations: AO, Angola; CI, Côte d’Ivoire; GH, Ghana.

Fig 9. Phylogenetic analysis of introgression in Ace1.

A) Maximum-Likelihood phylogenetic analysis of 690 Western African haplotypes from the Ace1 duplicated region (2,787 phased variants), using a GTR model with ascertainment bias correction, empirical state frequencies and four Γ rate categories. B-D) Id., using variants beyond the downstream duplication breakpoint (512 variants) and 1Mb upstream and downstream of the duplication (3,935 and 3,302). Tips are color-coded according to species and duplication presence. Source alignments and complete phylogenies with supports in all nodes are available in S14 and S15 Data. E) Distance in allelic frequencies between groups of duplicated specimens and wt A. coluzzii and A. gambiae, calculated using the three-population branch statistic in windows of 5,000 variants along 2R (see Methods). Includes estimated branch lengths (L) from within the duplication region.