Identification of a rapidly‐spreading triple mutant for high‐level metabolic insecticide resistance in Anopheles gambiae provides a real‐time molecular diagnostic for antimalarial intervention deployment

Abstract Studies of insecticide resistance provide insights into the capacity of populations to show rapid evolutionary responses to contemporary selection. Malaria control remains heavily dependent on pyrethroid insecticides, primarily in long lasting insecticidal nets (LLINs). Resistance in the major malaria vectors has increased in concert with the expansion of LLIN distributions. Identifying genetic mechanisms underlying high‐level resistance is crucial for the development and deployment of resistance‐breaking tools. Using the Anopheles gambiae 1000 genomes (Ag1000g) data we identified a very recent selective sweep in mosquitoes from Uganda which localized to a cluster of cytochrome P450 genes. Further interrogation revealed a haplotype involving a trio of mutations, a nonsynonymous point mutation in Cyp6p4 (I236M), an upstream insertion of a partial Zanzibar‐like transposable element (TE) and a duplication of the Cyp6aa1 gene. The mutations appear to have originated recently in An. gambiae from the Kenya‐Uganda border, with stepwise replacement of the double‐mutant (Zanzibar‐like TE and Cyp6p4‐236 M) with the triple‐mutant haplotype (including Cyp6aa1 duplication), which has spread into the Democratic Republic of Congo and Tanzania. The triple‐mutant haplotype is strongly associated with increased expression of genes able to metabolize pyrethroids and is strongly predictive of resistance to pyrethroids most notably deltamethrin. Importantly, there was increased mortality in mosquitoes carrying the triple‐mutation when exposed to nets cotreated with the synergist piperonyl butoxide (PBO). Frequencies of the triple‐mutant haplotype remain spatially variable within countries, suggesting an effective marker system to guide deployment decisions for limited supplies of PBO‐pyrethroid cotreated LLINs across African countries.


| INTRODUC TI ON
Insecticide resistance in disease vectors has become an influential model for understanding rapid contemporary evolution but, more importantly, identifying how resistance arises and spreads is crucial for disease control. Resistance to pyrethroid insecticides in African malaria vector mosquitoes has spread to near ubiquity (Ranson & Lissenden, 2016;WHO, 2019) and, though it is often difficult to demonstrate its impact on malaria infections , in some cases it has reached levels that threaten the effectiveness of vector control programmes (Staedke et al., 2020;Killeen & Ranson, 2018). A better understanding of resistance distribution and mechanisms will permit a more informed selection and deployment of insecticides to combat evolving mosquito populations. Whilst our understanding of the genetic basis of insecticide resistance in mosquitoes has advanced substantially (The Anopheles gambiae 1000 Genomes Consortium, 2020), especially in the important vector Anopheles funestus , molecular diagnostics for the major vectors in the An. gambiae complex remain limited to a handful of mutations  which explain a relatively small fraction of the variance in phenotype Mitchell et al., 2014) or which are now at such high frequency as to provide limited diagnostic resolution .
Long-lasting insecticidal nets (LLINs) are the principal tool for vector control to combat malaria, especially in sub-Saharan Africa (Bhatt et al., 2015). The majority of LLINs are treated only with pyrethroid insecticides, to which resistance is now widespread (Hancock et al., 2020). Though behavioural variation and physical or physiological modifications affecting insecticide uptake may sometimes play a role, pyrethroid resistance is caused predominantly by two distinct mechanisms. The first is resistance via point mutations in the target-site of the insecticide, for pyrethroids the Voltage-gated sodium channel (Vgsc), which results in decreased sensitivity to the insecticide (Clarkson et al., 2021); the second is metabolic resistance due to over-expression or altered activity of detoxification enzymes, of which the cytochrome P450 family is commonly considered most important (Müller et al., 2008Weedall et al., 2019. Cytochrome P450 activity is inhibited by the synergist piperonyl-butoxide (PBO), and bed nets incorporating PBO are effective against P450mediated resistance, as demonstrated by large-scale field trials (Staedke et al., 2020;Protopopoff et al., 2018). Given the continued operational use of pyrethroid-containing nets, it is vital that we understand the genetic mechanisms that may impact their efficacy, to optimize bednet deployment, preferably using information from rapidly-applied DNA markers. In advance of a randomized control trial of PBO-LLINs (Staedke et al., 2020), we sought to characterize pyrethroid resistance mechanisms in the primary malaria vector An. gambiae s.s.  in Uganda and Kenya.
The recent development of the Anopheles gambiae 1000 genomes project (Ag1000g) has led to a step change in our ability to identify DNA variation driven by selection pressure. We have been able to perform genome-wide searches for regions under recent natural selection in insecticide resistant populations across Africa, and work has shown that the strongest selective sweeps in the genome are all found around genes known to be important for resistance (The Anopheles gambiae 1000 Genomes Consortium, 2020, 2017). Furthermore, a whole-genome scan of copy number variants (CNVs) in the Ag1000g data revealed that increases in gene copy number were highly enriched in clusters of detoxification genes, pointing to a potentially widespread mechanism for increased gene expression (Lucas et al., 2019), which, in some cases, may elevate resistance. A number of gene duplications were observed around the Cyp6aa/Cyp6p gene family cluster on chromosome 2R, and the majority of these duplications included the gene Cyp6aa1 ( Figure 1). Cyp6aa1 has been found to be overexpressed in pyrethroid resistant populations in congeneric species (Ibrahim et al., 2018;Zhou et al., 2015;Kwiatkowska et al., 2013), but it has received very little attention compared to known insecticide-metabolizing genes such as Cyp6m2 (Mitchell et al., 2012;Edi et al., 2014), Cyp6p3 (Müller et al., 2008) and Cyp9k1 (Vontas et al., 2018) and its importance in resistance in An. gambiae remains unknown. from the Kenya-Uganda border, with stepwise replacement of the double-mutant  with the triple-mutant haplotype (including Cyp6aa1 duplication), which has spread into the Democratic Republic of Congo and Tanzania. The triple-mutant haplotype is strongly associated with increased expression of genes able to metabolize pyrethroids and is strongly predictive of resistance to pyrethroids most notably deltamethrin. Importantly, there was increased mortality in mosquitoes carrying the triple-mutation when exposed to nets cotreated with the synergist piperonyl butoxide (PBO). Frequencies of the triple-mutant haplotype remain spatially variable within countries, suggesting an effective marker system to guide deployment decisions for limited supplies of PBO-pyrethroid cotreated LLINs across African countries.

K E Y W O R D S
adaptation, contemporary evolution, disease biology, ecological genetics, insects In this study, we examined a strong selective sweep detected in the Cyp6aa/Cyp6p genomic region in samples of An. gambiae s.s. from Uganda and Western Kenya. We found that the sweep was closely associated with three mutations (a SNP in Cyp6p4, a duplication of Cyp6aa1 and a partial transposable element insertion termed ZZB-TE) in tight physical and statistical linkage. The triple-mutant haplotype was associated with a high-level of pyrethroid resistance, most notably to deltamethrin. The three mutations appeared sequentially, leading to successive selective sweeps, with the triple-mutant haplotype replacing earlier variants and then apparently spreading rapidly across East and Central Africa. We show that this haplotype is under positive selection and is associated with increased expression of key cytochrome P450s and through recombinant protein expression using both an E.
coli and an Sf9-baculovirus system we show that both CYP6AA1 and CYP6P4 are capable of metabolizing pyrethroid insecticides.
We also used the haplotype groupings identified above to calculate the Garud H statistics (Garud et al., 2015) and the haplotypic diversity at the Cyp6aa/Cyp6p cluster locus (coordinates 2R:28,480,576-2R:28,505,816). Specifically, we used the moving_ garud_h and moving_haplotype_diversity functions in scikit-allel to obtain a series of estimates for each statistic in blocks of 100 variants located within the cluster, and used a block-jackknife procedure to calculate averages and standard errors of each estimate (jackknife function in scikit-allel misc module).

| Molecular screening of colony and wildcaught mosquitoes
Locked-nucleic acid (LNA) probe-based PCR diagnostics were designed for all three mutations (Appendix S1).

| Mosquito colonies
Genotype: phenotype association testing was performed using two colonies of An. gambiae s.s., BusiaUG (resistant) and Mbita (susceptible). The BusiaUG strain was established in the lab in November 2018 from Busia, eastern Uganda and exhibits high resistance to pyrethroids and organochlorines (<10% mortality following WHO exposure test; fixed for Vgsc-995S resistance mutation), but full susceptibility to organophosphates and carbamates. The Mbita strain was first colonized from Mbita Point, Kenya in 1999, and is fully susceptible to pyrethroids. Colonies were reared in insectaries targeted to 25-27°C and 70%-80% relative humidity.
Freshly emerged females from the Mbita line were mated with 3-5 day-old BusiaUG males and then blood fed. The offspring from this cross were then crossed back to the parental BusiaUG line.
This design was chosen as resistance variants are often recessive.
The resultant backcrossed 3-5-day-old females were exposed for 1 h to deltamethrin, permethrin or α-cypermethrin (the three insecticides most commonly used on LLINs) or DDT (a nonpyrethroid sodium channel antagonist), following WHO standard procedures (WHO, 2016). Mosquitoes were maintained on a 10% sugar solution after exposure and mortality was recorded 24 h post exposure.

| Democratic Republic of Congo
Anopheles gambiae s.s. were obtained from the President's Malaria Initiative supported entomological surveillance project (Wat'senga et al., 2020) and from collections conducted by Lynd et al. (2018).
Mosquitoes were collected by human landing catch or pyrethrum spray collection from 15 locations between 2013 and 2018.
Resistance-phenotyped individuals were also obtained from Pwamba, Bassa and Fiwa in Nord Ubangi Province in 2016. These mosquitoes were assessed for susceptibility to deltamethrin or permethrin using a standard WHO tube assay or a cone bioassay where Permanet 3.0 (deltamethrin plus PBO) or Olyset Plus (permethrin plus PBO) were the test nets (Kwiatkowska et al., 2013).

| Kenya and Uganda
Collection details for specimens from contiguous areas of western Kenya and eastern Uganda have been published previously (Lucas et al., 2019;Mitchell et al., 2012;Edi et al., 2014).

| Tanzania
Mosquito collections were conducted in Geita, Bagamoyo and

Muleba districts of Tanzania in 2018.
Mosquitoes from all collections were genotyped at the ZZB-TE, Cyp6p4-236M and Cyp6aap-Dup1 loci and genotype: phenotype association testing was performed by Fisher's exact tests.

| Analysis of temporal change
The most complete time series of ZZB-TE, Cyp6p4-236 M and Cyp6aap-Dup1 allele frequencies were available from Kabondo, DRC; western Kenya and eastern Uganda. To model allele frequency changes over time we estimated the three parameters of the standard recursive population genetic model (allele frequency at time zero, selection coefficient and dominance coefficient) using a maximum likelihood approach assuming a binomial distribution: an approach previously applied to insecticide target-site resistance mutations (Lynd et al., 2010). The analysis was performed in R (http:// www.r-proje ct.org). An estimated generation time of one calendar month was used as in previous studies (Lynd et al., 2010).

| Recombinant protein and insecticide metabolism
The CYP6P4-236M variant, and its redox partner cytochrome P450 reductase gene (CPR), were expressed in E. coli as per standard protocols (Yunta et al., 2019) (Appendix S2). Initial efforts to generate recombinant CYP6AA1 in an E. coli system with optimized codon usage failed, and we therefore used an Sf9-baculovirus-based expression system. Since P450 catalytic activity is dependent on electrons supplied by NADPH via CPR, insecticide metabolism was assayed with cell pellets of CYP6P4/CPR or CYP6AA1/CPR in the presence or absence of NADPH. The depletion of the substrate and the appearance of metabolites were monitored by reverse-phase HPLC (Appendix S3).

| Estimation of gene expression of key P450s in triple mutant haplotype
To determine whether the presence of the triple mutant haplotype type was associated with differential expression of genes, we examined individual females from the BusiaUg colony (Triple mutant Freq. 0.297; 95% CI 0.233-0.370). Two legs were removed from individual mosquitoes for DNA analysis with the remaining mosquito kept for RNA analysis. DNA was extracted from legs by boiling in STE buffer at 95°C for 90 min and individuals were genotyped using the LNA qPCR assays. RNA was then extracted individually from eight mosquitoes in each genotypic group -homozygotes for the triple mutant haplotype, wild-type homozygotes, and heterozygotes, using the Arcturus Picopure RNA isolation kit (Thermofisher). We then performed SYBR green based qPCR to measure the expression of Cyp6aa1 and Cyp6p4 together with the known resistance-linked variant Cyp6p3 using the housekeeping genes 40s ribosomal protein S7 (AGAP010592) and elongation factor Tu (AGAP005128) for normalization. The ΔΔCT values were tested for normality and homogeneity of variances using the Shapiro-Wilks test, and the Bartlett test, respectively. A significant difference in gene expression between the genotypic groups was determined by a two-tailed two-sample Student's t-test on ΔΔCT values, with a threshold of p = .05.

| RE SULTS
Hierarchical clustering of 206 Ag1000g haplotypes (n = 103 individual females) from An. gambiae s.s. from Uganda resulted in identification of a common haplotype (n = 122 out of 206) around the Cyp6aa/Cyp6p gene cluster, putatively representing a swept haplotype in this region. To characterize the signatures of selection in Ugandan haplotypes around this cluster, we examined the profile of extended haplotype homozygosity around the position of the Cyp6p4-236 M SNP and around the CNV in Cyp6aa1. In both cases, we found that the putative swept haplotype had longer stretches of homozygosity than wild-type haplotypes ( Figure 2).
In addition, we found that An. gambiae s.s. from Uganda had reduced haplotypic diversity along the entire Cyp6aa/Cyp6p gene cluster (h = 0.339746 ± 0.005664 standard error) and a combination of Garud's H statistics that was indicative of a hard selective sweep in this region (high H 12 = 0.821867 ± 0.006308 SE; low H 2 / H 1 = 0.016779 ± 0.000228 SE) (Garud et al., 2015). These results strongly suggest that the haplotypes we have identified have undergone a selective sweep. We then used iterative read mapping of individuals homozygous for the sweep to search for additional mutations that might be distinctive of the haplotype. This revealed that a partial copy of a Ty3/Gypsy Zanzibar transposon insertion F I G U R E 2 Selective sweep around the Cyp6aa/Cyp6p cluster in Anopheles gambiae from Uganda. Upper panel Showing the relationship between the selective sweep ("Sweep")and the three mutations Cyp6p4-236M SNP ("P4"), the ZZB-TE insertion ("ZZB") and Cyp6aap-Dup1 ("Dup1"). Each vertical bar represents a single haplotype, colour-coded to show whether it is a copy of the swept haplotype (blue) and whether it carries the SNP (purple), the TE insertion (green) and the duplication (orange). The swept haplotype and the Cyp6p4-236M SNP overlap almost completely, while the duplication is found on a subset of these haplotypes. Lower panel extended haplotype homozygosity (EHH) plots around the CYP6P4-236M and CYP6AA1-Dup1 variants show slower loss of homozygosity in the swept haplotypes than in the wild-type western Kenya and eastern Uganda, the 95% CIs for each observed data point were calculated according to Newcombe (1998). The model data, illustrated with blue or red curves, were generated from simultaneous maximum likelihood estimates of initial frequency and selection and dominance coefficients which were then used to parameterize standard recursive allele frequency change equation. The purple dotted lines indicate approximate timings of LLIN mass coverage campaigns in the sampled locations. The lower panel shows data from DRC tracking the emergence and spread of the triple mutant haplotype. Triangles indicate collection locations within each province and the scale bar indicates triple mutant frequency  heterogeneity, with very low frequencies in the more southerly provinces (Figure 3). wild-type genotypes, respectively. Figure 4). As a control we examined a neighbouring, very commonly resistance-associated gene,   Table 2). No association was found in samples exposed to permethrin (24 h WHO tube assay) or permethrin-treated Olyset Plus nets (3-min WHO cone assay) ( Table 2). Complete linkage of the three mutants in the BusiaUG colony and the DRC wild-caught collections precludes determination of the relative contribution of each of the three mutations to the resistance phenotype but taken together these results demonstrate a strong impact of the triple mutant on the efficacy of pyrethroid resistance.

| DISCUSS ION
We have identified a series of adaptive mutations in An. gambiae culminating in a triple-mutant haplotype with a large effect on pyrethroid resistance and which is spreading rapidly across East and Central Africa. The mutation that is probably the oldest in this in Vgsc of single origins and spread over many thousands of kilometres (Clarkson et al., 2021).
Second generation nets treated with the synergist PBO were shown to be much more effective than conventional nets in killing mosquitoes in populations where the triple mutant haplotype is present ( Figure 5). Therefore the use of PBO bednets should be prioritized in the regions where this mutation is present. A strong corollary of this finding comes from the cluster randomized control trial conducted in Uganda, where the mutation is at high frequency, which demonstrated that malaria parasite prevalence in children <10 years old (12% vs. 14%; Prevalence ratio = 0.84, 95% CI: 0.72-0.98; p = .029) and mean number of mosquitoes (density ratio = 0.25, 95% CI: 0.18-0.3; p < .0001) per house were significantly lower in villages that had received PBO LLINs relative to standard LLINs (Staedke et al., 2020). There is some evidence that the haplotype may be less strongly associated with resistance to permethrin than deltamethrin (and perhaps also alphacypermethrin), although both pyrethroids were metabolized by Cyp6aa1 and Cyp6p4.

F I G U R E 5
Patterns of insecticide resistance in Anopheles gambiae following exposure to permethrin and deltamethrin treated nets. Female Anopheles gambiae s.s. from the Democratic Republic of Congo (DRC) and Uganda exhibited very high resistance (low % mortality) in WHO cone assays to permethrin (Olyset) and deltamethrin (PermaNet 2.0) LLINs. There was an increase in mortality following exposure to an increased concentration of deltamethrin (PermaNet 3.0 side) and very high mortality following exposure to pyrethroid plus the P450 inhibitor piperonyl butoxide (PBO) Olyset+ and PermaNet 3.0 top. 95% confidence intervals are shown Our results highlight the importance of gene duplications for the evolution of insecticide resistance. In An. gambiae, duplications have recently been shown to be concentrated in regions associated with metabolic resistance, and over 40 such duplications have been described across the genome (Lucas et al., 2019). Thirteen different duplications have so far been described that encompass Cyp6aa1 (Figure 1), both in West and East Africa and in the two sister-species An. gambiae ss and An. coluzzii (Lucas et al., 2019). It seems likely that these other Cyp6aa1 duplications are also associated with pyrethroid resistance. For example, in An. coluzzii sampled from a highly insecticide-resistant population from Cote d'Ivoire, greater than 95%  malaria vector, An. funestus, from west and central Africa (Ibrahim et al., 2018). The Cyp6aa1 orthologue in An. funestus, which shares an 87% identity with An. gambiae, was also observed to metabolize permethrin and deltamethrin and, when expressed in transformed Drosophila, was associated with significant increases in resistance to both permethrin and deltamethrin relative to control mosquitoes (Ibrahim et al., 2018). This evolution of multiple Cyp6aa1 duplications suggests this is an important Africa wide resistance mechanism.
There are now several cases of insecticide resistance evolution where an initial mutation in a genomic region is followed by the spread of additional mutations on the resistant haplotype background (Clarkson et al., 2021;Schmidt et al., 2010;Jones et al., 2012;Djogbenou et al., 2009;Grau-Bové et al., 2020 (Schmidt et al., 2010) and are abundant in mosquito genomes (Neafsey et al., 2015;Nene et al., 2007). Interestingly, in the common house mosquito, Culex pipiens, TEs have been found in the flanking regions of the Ester locus, a genomic region in which many independent gene duplications have arisen and spread worldwide in response to selection from organophosphates (Buss & Callaghan, 2004). Clearly, both TEs and gene duplications are an understudied, yet common source of variation that may have important implications for vector control efforts. The appearance and rapid spread of the three mutations described here is broadly coincident with the scale up of LLIN coverage in DRC, Kenya and Uganda (Bertozzi-Villa et al., 2021). The haplotype is a strong predictive marker of high-level resistance to pyrethroids and the Cyp6aap-Dup1mutation which is the key identifier of the triple mutant, is easily screened with a single diagnostic assay.
We suggest that this assay should be used for both insecticide resistance monitoring strategies and for informing LLIN selection (The