Xavier Grau-Bové
email: xavier.graubove@gmail.com

Eric Lucas

Dimitra Pipini

Emily Rippon

Arjèn Van't Hof

Edi Constant

Samuel Dadzie

Alexander Egyir-Yawson

John Essandoh

Joseph Chabi

Luc Djogbénou

Nicholas J Harding

Alistair Miles

Dominic Kwiatkowski

Martin J Donnelly

David Weetman
email: david.weetman@lstmed.ac.uk

The Anopheles Gambiae

Genomes Consortium

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

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Pirimiphos-methyl is an organophosphate insecticide that is widely used in control interventions against populations of the malaria vector Anopheles, especially in Africa [START_REF] Oxborough | Trends in US President's Malaria Initiative-funded indoor residual spray coverage and insecticide choice in sub-Saharan Africa (2008-2015): urgent need for affordable, long-lasting insecticides[END_REF][START_REF] Dengela | Multi-country assessment of residual bio-efficacy of insecticides used for indoor residual spraying in malaria control on different surface types: results from program monitoring in 17 PMI/USAIDsupported IRS countries[END_REF]. Since 2013, the World Health Organization (WHO) has recommended the use of pirimiphos-methyl for indoor residual spraying (IRS) interventions, the major anti-vector strategy in malaria control after treated bednet distribution [START_REF] Oxborough | Trends in US President's Malaria Initiative-funded indoor residual spray coverage and insecticide choice in sub-Saharan Africa (2008-2015): urgent need for affordable, long-lasting insecticides[END_REF]World Health Organization 2013). Strategic approaches to vector control often rely on the use of multiple insecticides to avoid or overcome the recurrent emergence of resistance in natural populations (World Health Organization 2018a). In that regard, various insecticide classes have been used in IRS, with pyrethroids-which target the voltage-gated sodium channel-being the dominant choice until recently [START_REF] Sherrard-Smith | Systematic review of indoor residual spray efficacy and effectiveness against Plasmodium falciparum in Africa[END_REF][START_REF] Van Den Berg | Global Trends in the Use of Insecticides to Control Vector-Borne Diseases[END_REF]. However, the increase in pyrethroid resistance in Anopheles populations [START_REF] Ranson | Pyrethroid resistance in African anopheline mosquitoes: what are the implications for malaria control?[END_REF]) has led to a progressive replacement with acetylcholinesterase-targeting insecticide classes, first the carbamate bendiocarb, and, latterly, the organophosphate pirimiphos-methyl [START_REF] Oxborough | Trends in US President's Malaria Initiative-funded indoor residual spray coverage and insecticide choice in sub-Saharan Africa (2008-2015): urgent need for affordable, long-lasting insecticides[END_REF]. Pirimiphos-methyl is the active ingredient in the most widely used insecticide for IRS in Africa, the spray formulation Actellic, which is highly effective and has strong residual performance [START_REF] Dengela | Multi-country assessment of residual bio-efficacy of insecticides used for indoor residual spraying in malaria control on different surface types: results from program monitoring in 17 PMI/USAIDsupported IRS countries[END_REF][START_REF] Sherrard-Smith | Systematic review of indoor residual spray efficacy and effectiveness against Plasmodium falciparum in Africa[END_REF][START_REF] Wagman | An observational analysis of the impact of indoor residual spraying with non-pyrethroid insecticides on the incidence of malaria in Ségou Region[END_REF][START_REF] Abong'o B | Impact of indoor residual spraying with pirimiphos-methyl (Actellic 300CS) on entomological indicators of transmission and malaria case burden in Migori County, western Kenya[END_REF]). However, resistance has recently been reported in several populations of African Anopheles s.l. [START_REF] Kisinza | Multiple insecticide resistance in Anopheles gambiae from Tanzania: A major concern for malaria vector control[END_REF][START_REF] Chukwuekezie | Resistance status of Anopheles gambiae (s.l.) to four commonly used insecticides for malaria vector control in South-East Nigeria[END_REF], and though control failures have yet to be reported, it represents a clear threat to the efficacy of IRS strategies.

Mechanisms of resistance to pirimiphos-methyl are poorly understood, but organophosphates, as well as carbamates, all block the action of the acetylcholinesterase enzyme (ACE1, encoded by the Ace1 gene in mosquitoes) via competitive binding to its active site [START_REF] Oakeshott | Comparing the organophosphorus and carbamate insecticide resistance mutations in cholin-and carboxyl-esterases[END_REF]. Studies in culicine and anopheline mosquitoes have found that non-synonymous mutations in Ace1 can result in resistance to carbamates and organophosphates other than pirimiphos-methyl [START_REF] Weill | Insecticide resistance in mosquito vectors[END_REF][START_REF] Weill | The unique mutation in ace-1 giving high insecticide resistance is easily detectable in mosquito vectors[END_REF][START_REF] Feyereisen | Genotype to phenotype, the molecular and physiological dimensions of resistance in arthropods[END_REF]. The most common mutation in Anopheles is a glycine to serine mutation in codon 280 (G280S, also known as G119S after the codon numbering of a partial crystal structure from the electric ray Torpedo californica [START_REF] Weill | Insecticide resistance in mosquito vectors[END_REF][START_REF] Greenblatt | The complex of a bivalent derivative of galanthamine with torpedo acetylcholinesterase displays drastic deformation of the active-site gorge: implications for structure-based drug design[END_REF][START_REF] Feyereisen | Genotype to phenotype, the molecular and physiological dimensions of resistance in arthropods[END_REF])), which is located near the active site gorge of ACE1 [START_REF] Cheung | Structure of the G119S Mutant Acetylcholinesterase of the Malaria Vector Anopheles gambiae Reveals Basis of Insecticide Resistance[END_REF]) and decreases its sensitivity to organophosphates and carbamates.

Results

Conservation and distribution of Ace1 resistance mutations in Anopheles

We examined the frequency and distribution of the two Ace1 mutations that have been associated with organophosphate and carbamate resistance in A. gambiae and A. coluzzii: the G280S non-synonymous single nucleotide polymorphism (SNP), and copy number variation (CNV) polymorphisms of Ace1 and the surrounding genomic region. G280S is sometimes known as G119S [START_REF] Feyereisen | Genotype to phenotype, the molecular and physiological dimensions of resistance in arthropods[END_REF] based on its position in the truncated crystal structure of its homolog in the electric ray Torpedo californica, where ACE1 protein structure was first elucidated [START_REF] Greenblatt | The complex of a bivalent derivative of galanthamine with torpedo acetylcholinesterase displays drastic deformation of the active-site gorge: implications for structure-based drug design[END_REF]). Due to a culicine-specific N-terminal insertion in ACE1, the exact position of this conserved codon differs among animal orthologs (Supplementary Material SM1). We provide a list of homologous codon positions for Ace1 orthologs from selected animal species, including common insect vectors (Supplementary Material SM2). Henceforth, we will use A. gambiae s.l.-based codon coordinates and refer to this SNP as G280S (wild-type allele: 280G; resistance allele: 280S; gene accession number: AGAP001356-RA in AgamP4.12).

In the Anopheles 1000 Genomes cohort (Phase 2, n = 1142 genomes; Figure 1) (The Anopheles gambiae 1000 Genomes Consortium 2019), the 280S resistance allele is present across West African populations of A. coluzzii (CoB te d'Ivoire, Burkina Faso, Ghana) and A. gambiae (Burkina Faso, Ghana, Guinea), with the highest frequencies being attained in Ghanaian A. gambiae (67%) and Ivorian A. coluzzii (44%; Figure 1A), which is consistent with previous results [START_REF] Djogbénou | Evidence of Introgression of the ace-1R Mutation and of the ace-1 Duplication in West African Anopheles gambiae s[END_REF][START_REF] Clarkson | The genetic architecture of target-site resistance to pyrethroid insecticides in the African malaria vectors Anopheles gambiae and Anopheles coluzzii[END_REF][START_REF] Alou | Distribution of ace-1R and resistance to carbamates and organophosphates in Anopheles gambiae s.s. populations from CoB te d'Ivoire[END_REF][START_REF] Essandoh | Acetylcholinesterase (Ace-1) target site mutation 119S is strongly diagnostic of carbamate and organophosphate resistance in Anopheles gambiae s.s. and Anopheles coluzzii across southern Ghana[END_REF]. Ace1 CNVs have a similar distribution to the SNP in West African populations, as they overwhelmingly overlap with specimens also carrying 280S alleles (Figure 1B).

The combination of CNVs and 280S alleles results in multiple genotypes being observed in the Ace1 locus, defined by a variable number of 280S copies. We used the fraction of sequencing reads supporting the wt and 280S alleles to estimate the number of 280S alleles in each sample (Figure 2A and 2B), which revealed that most specimens with duplications carry both wt and 280S alleles at various frequencies (n = 113, orange box in Figure 2A), whereas the vast majority of non-duplicated specimens only have wt alleles (n = 1021, purple box in Figure 2A). In addition, we identify six high-copy, 280S homozygous specimens, present only in the Ghanaian A. gambiae population in this dataset (Figure 1B). Virtually all Ace1 CNVs identified in the 1000 Genomes cohort share the same duplication breakpoints [START_REF] Lucas | Whole-genome sequencing reveals high complexity of copy number variation at insecticide resistance loci in malaria mosquitoes[END_REF], spanning a region ca. 200 kbp that includes a total of 11 genes (Supplementary Material SM5). The only exception is a single A. gambiae specimen from Guinea that carries a unique wt-homozygous duplication (asterisk in Figure 2B). This Guinea-specific CNV is shorter than the Ace1 duplication found in all other samples and has different breakpoints (ca. 70 kbp, including Ace1 and one other gene; Supplementary Material SM4), implying an independent origin. Henceforth, all mentions of Ace1 CNVs will refer to the major duplication. The absence of specimens carrying the major duplication and lacking 280S alleles strongly implies that this CNV contains both 280S and wt alleles on the same chromosome, and thus results in permanent heterozygosity. Figure 2C summarises the four model haplotypes that can result in the genotypes observed in the 1000 Genomes dataset: (i) non-duplicated wt; (ii) heterozygous and (iii) 280Shomozygous major duplications; and (iv) the wt minor duplication from Guinea.

In addition to G280S, we found seven non-synonymous SNPs in Ace1 with at least 1% frequency in at least one population (Figure 1A). None were in linkage disequilibrium with G280S (Supplementary Material SM3), nor have any of them been previously associated with resistance [START_REF] Feyereisen | Genotype to phenotype, the molecular and physiological dimensions of resistance in arthropods[END_REF]. The most salient of these additional mutations was a serine-to-alanine SNP in codon 65 (S65A), with frequencies of 30-50% in multiple West African populations. Codon 65 maps to the anopheline-specific N-terminal insertion of ACE1 proteins (Supplementary Material SM1), which lacks predicted secondary structure and is far away from the active-site gorge of the enzyme [START_REF] Cheung | Structure of the G119S Mutant Acetylcholinesterase of the Malaria Vector Anopheles gambiae Reveals Basis of Insecticide Resistance[END_REF]). This position is not conserved across ACE1 orthologs, being variously encoded as alanine, serine or threonine in different Anopheles species (Supplementary Material SM1). in 2012, that had been tested for resistance prior to genome sequencing. The G280S SNP and CNV co-occur at high frequencies in this population, with 87.3% of specimens being duplicated heterozygotes with a one or more 280S alleles (Figure 3A). Specimens with at least one 280S allele had a survival rate of 50%, as opposed to 0% in the wt homozygotes (Figure 3B), suggesting that 280S is linked to pirimiphos-methyl resistance in this population (p = 3.7 × 10 -3 in Fisher's exact test, Woolf-corrected odds ratio = 19.0 [95% confidence interval = 1.1 -340.6]) but it does not fully explain the resistance phenotype. The only other non-synonymous SNP present in Ivorian A. coluzzii, S65A, was not associated with pirimiphos-methyl resistance (p = 0.605 Fisher's exact test; Supplementary Material SM6A, SM6B).

F 3 5 Y ( 1 0 4 T > A ) A 3 6 E ( 1 0 7 C > A ) E 4 4 D ( 1 3 2 G > T ) G 6 1 D ( 1 8 2 G > A ) S 6 5 A ( 1 9 3 T > G ) Angola A. coluzzii

Burkina Faso A. coluzzii

Burkina

Pirimiphos-methyl resistance in Ivorian

We found a strong bias towards higher Ace1 copy numbers in resistant specimens, all of which had at least 3 copies of the Ace1 locus and most of them had 4 or more (Figure 3B; odds ratio from a binomial generalised linear model [GLM] = 16.9 [5.1 -56.6]; p = 6.4 × 10 -10 in a χ 2 test compared to a null model).

Concordantly, the estimated number of 280S alleles was also associated with resistance (GLM odds ratio = 10. 7 [3.6 -31.7]; p = 8.9 × 10 -10 in χ 2 test): survival rates increased with higher number of 280S copies (Figure 3B), with a threshold apparent between individuals with zero or one copies of 280S (11.1% survival) and those with two or more (77% survival and two wt alleles (n = 36), 280S provided no survival benefit compared to non-duplicated wt specimens (Fisher's exact test p = 0.553; Woolf-corrected odds ratio = 3.6 [0.2 -74.3]).

We also examined the combined effects of all the mutations detected in the Ace1 locus (280S and S65A presence, CNVs and number of 280S alleles) in Ivorian A. coluzzii. We used a binomial GLM and a step-wise procedure to remove non-informative variants according to the Bayesian Information Criterion (see Methods), finding that the number of Ace1 copies was the minimal sufficient marker to predict pirimiphosmethyl resistance in this population (odds ratio = 16.9 [5.1 -56.9]; p = 6.43 × 10 -10 in a χ 2 test; details in Supplementary Material SM6B).

Pirimiphos-methyl resistance diagnostics are repeatable across A. gambiae and A. coluzzii populations We investigated the extent to which the role of Ace1 and 280S copy numbers in pirimiphos-methyl resistance could be generalised to other populations beyond the A. coluzzii from Tiassalé. Specifically, we surveyed Ace1 mutations and phenotypic resistance in eight West African populations of A. gambiae and A. coluzzii from additional locations in CoB te d'Ivoire, Togo, and Ghana (Table 1 andTable 2; list of specimens in Supplementary Material SM6C and SM6D).

Pirimiphos-methyl resistance was associated with G280S mutations in seven out of eight populations (p < 0.01 in a χ 2 test comparing GLM of each genotype to a null model; Table 1). However, whilst 280S homozygotes were strongly associated with resistance in populations where they were present in significant numbers (e.g. 95.8% survival rate in A. coluzzii from Korle-Bu), heterozygotes often exhibited lower, intermediate survival rates (e.g. 49.6% in A. coluzzii from Korle-Bu; six out of eight populations had survival rates <60%; Table 1) similar to that of the A. coluzzii population from Tiassalé (CoB te d'Ivoire) mentioned above (Figure 3A).

Following our findings for the Tiassalé population (Figure 3), we further investigated the phenotypic variation in duplicated heterozygotes using two CNV-related variables: (i) the ratio of 280S to 280G alleles (measured as the ratio of FAM-to-HEX dye signal in a TaqMan qPCR assay) [START_REF] Bass | Development of high-throughput real-time PCR assays for the identification of insensitive acetylcholinesterase (ace-1R) in Anopheles gambiae[END_REF]; and (ii) the estimated number of Ace1 copies relative to a non-duplicated and wt Kisumu A. gambiae colony (Edi et al. 2014a) (see Methods). Pirimiphos-methyl resistance was associated with CNV-related measures in heterozygotes from A.

gambiae and A. coluzzii populations, according to generalised linear models that included 280S allelic ratios and/or Ace1 copy number as potential informative markers (Table 2). In A. coluzzii from Korle-Bu and Weija, as well as A. gambiae from Madina, both the 280S allelic ratio and the Ace1 copy number were included in the minimal model (at p < 0.01 in a χ 2 comparison with a null model); whereas in A. coluzzii from Aboisso the allelic ratio was the only variable associated with resistance (p = 1.3 × 10 -3 in a χ 2 ). In the Obuasi A. gambiae heterozygotes, the minimal model only included Ace1 copy number (at p = 0.038) but neither measure was strongly predictive, as was the case in A. gambiae from Aboisso and Baguida. Overall, these results indicate that the combination of 280S allelic ratios and CNVs provide a similar pirimiphosmethyl resistance mechanism in heterozygotes from both species. However, diagnostic capacity appears to vary among populations, suggesting the possible existence of additional resistance mechanisms.

Table 2. Association of Ace1 duplications with pirimiphos-methyl resistance in heterozygous specimens from seven West African populations of A. coluzzii and A. gambiae. For each population, we report number of sampled specimens (n), and the results of three binomial generalised linear models (GLM): (i) using the number of Ace1 copies as the predictor variable, (ii) using the 280S-to-280G allelic ratio, and (iii) a minimal model obtained with step-wise reduction, according to the Bayesian Information Criterion (BIC). For each model, we report the p-value from a comparison with a null model (ANOVA χ 2 test) and odds ratios (OR) with 95% confidence intervals for the variables included in the model.

) 1.3 × 10 -3 - 52.9 (0.01 - 2.2 × 10 5 ) 1.3 × 10 -3 coluzzii Korle-Bu (GH) 46 1.8 (1.2 -2.6) 1.0 × 10 -4 2.8 (1.6 -5.1) 4.1 × 10 -7 1.5 (0.9 -2.3) 2.3 (1.3 -4.0) 3.0 × 10 -7 coluzzii Weija (GH) 21 0.5 (0.3 -1.0) 0.013 3.4 (0.4 -28.6) 0.021 0.4 (0.1 -1.1) 50 (0.3 -9505.0) 2.9 × 10 -3 gambiae Aboisso (CI) 31 1.0 (0.6 -1.5) 0.933 2.6 (0.3 -24.9) 0.382 - - - gambiae Madina (GH) 19 5.2 (1.1 -23.9) 7.7 × 10 -3 6.5 × 10 85 (0 -Inf) 2.1 × 10 -6 8.7 × 10 -20 (0 -Inf) 5.5 × 10 45 (0 -Inf) 2.6 × 10 -9 gambiae Obuasi (GH) 26 1.5 (0.9 -2.3) 0.038 3.1 (0.85 -10.8) 0.038 1.5 (0.9 -2.3) - 0.038 gambiae Weija (GH) 0 - - - - - - -
gambiae Baguida (TG) 13 1.4 (0.7 -2.6) 0.350 2.9 (0.1 -69.4) 0.498 ---

Genome-wide identification of pirimiphos-methyl resistance variants in Ivorian A. coluzzii

We examined the genomes of the 71 A. coluzzii specimens from CoB te d'Ivoire to identify additional genetic variants-other than Ace1 G280S and duplications-that could be linked to pirimiphos-methyl resistance. All Ivorian specimens had been collected in the same location and time (Tiassalé, 2012), before the widespread adoption of pirimiphos-methyl in IRS interventions [START_REF] Oxborough | Trends in US President's Malaria Initiative-funded indoor residual spray coverage and insecticide choice in sub-Saharan Africa (2008-2015): urgent need for affordable, long-lasting insecticides[END_REF]World Health Organization 2013). The absence of a long period of pirimiphos-methyl use prior to collection implies that adaptation to this insecticide would presumably be caused by cross-resistance with previously employed insecticides and/or standing variation, rather than novel selective sweeps. In accordance with this expectation of low differentiation, a principal component analysis of genetically unlinked variants (see Methods) did not reveal population stratification between resistant and susceptible mosquitoes (Figure 4), and average Hudson's F ST between them was low along all chromosomal arms (F ST ≈ 0%).

We thus aimed to identify signals of selection in resistant A. coluzzii using the population branch statistic (PBS, Figure 5A), which identifies regions with an excess of genetic differentiation in a focal population (here, resistant Ivorian mosquitoes) relative to a basal level of differentiation between a closely related population (susceptible Ivorian mosquitoes) and a more distant control (a population of A. coluzzii from Angola; see Methods). PBS is a powerful test to detect recent selection acting on incomplete sweeps and standing variation [START_REF] Yi | Sequencing of 50 Human Exomes Reveals Adaptation to High Altitude[END_REF][START_REF] Malaspinas | A genomic history of Aboriginal Australia[END_REF], and is therefore well-suited to investigate a scenario-like ours-in which the population of interest has not yet diverged.

We estimated PBS along genomic windows to identify regions in which resistant Ivorian specimens had an excess of genetic differentiation (see Methods). In total, we identified 43 genes within six regions with excess differentiation in resistant A. coluzzii (PBS > 0, FDR < 0.001; see Methods; Figure 5A). This set of candidates contained eleven mostly uncharacterised genes located downstream of Ace1 and within the duplicated region, two ionotropic receptors (GLURIIb and GLURIIa), and several protesases, among others (Supplementary Material SM7). However, genetic differentiation between resistant and susceptible specimens in each of these genes was low (F ST < 6% in all cases; Supplementary Material SM7), which suggests that they are unlikely to be directly associated with the resistance phenotype. The Ace1 gene itself was on the verge of the significance threshold (FDR = 3.2 × 10 -3 , PBS = 0.031) and exhibited low differentiation (F ST = 1.2%). This apparent contradiction regarding the Ace1 phenotypic association can be explained by the fact that pirimiphos-methyl resistance is caused by a combination of mutations-a SNP within a heterozygous duplication-that is not captured by the diploid allelic frequencies used in PBS calculations [START_REF] Yi | Sequencing of 50 Human Exomes Reveals Adaptation to High Altitude[END_REF]). To overcome this limitation, we performed a k-mer association study comparing the resistant and susceptible Ivorian A. coluzzii. k-mer association studies test for differential frequencies of tracts of nucleotides of length k in different groups of genomes, and are able to identify genetic variation patterns linked to both SNPs and structural variants, such as CNVs, without requiring a reference annotation [START_REF] Rahman | Association mapping from sequencing reads using k-mers[END_REF]. We identified ca. 767 million different k-mers of length k = 31 bp across all 71 samples (see Methods). Among these, 9,603 k-mers were significantly enriched in resistant specimens (Spearman correlation test, FDR < 0.001). The 9,603 significant k-mers were assembled into 446 unique sequences composed of more than one k-mer (median length = 54 bp), which we then aligned to the A. gambiae reference genome. We retained sequences that could be mapped to chromosomes 2, 3 or X (409 out of 446)

12

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • •• • •• • • • 0 100
for further analysis (listed in Supplementary Material SM8). Among these, the vast majority (n = 389, 95.1%) aligned to the Ace1 duplicated region, while the remaining significant sequences aligned in scattered regions across the rest of the genome (n = 20, 4.9%; Figure 5B). Yet, 19 out of 20 k-mers that did not map to Ace1 had a very similar frequency profile in resistant and susceptible samples than the 389

Ace1-linked k-mers, being absent and present in the same genomes (Figure 5C). In fact, the k-mer frequencies of these 19 sequences was strongly correlated with Ace1 copy number (Pearson correlation r > 0 and p < 1 × 10 -4 ; Figure 5D). This suggests that these 19 k-mers represent a non-independent signal that also reflects the association of Ace1 mutations with resistance, and can be parsimoniously explained by non-specific sequence alignments (see Discussion). Only one remaining sequence did not correlate with Ace1 copy number (Pearson correlation r < 0.0, maroon arrow in Figure 5C and5D), indicating that it may represent another variant involved in resistance. The primary alignment of this sequence mapped to a non-coding region of chromosomal arm 2L (from 2L:8,662,023), with no proximity to any gene of known function. Altogether, these results support the conclusion that Ace1 mutations are the primary driver of pirimiphos-methyl resistance in this A. coluzzii population.

Selection and introgression of Ace1 duplications in A. gambiae and A. coluzzii

A high degree of geographical and phylogenetic overlap is evident between Ace1 duplications and 280S

alleles across four countries (CoB te d'Ivoire, Ghana, Burkina Faso and Guinea) and two different species (A.

coluzzii and A. gambiae; Figure 2A). Previous studies provide two key insights to understand this pattern.

First, Ace1 duplications across West African populations have concordant breakpoints [START_REF] Assogba | The ace-1 Locus Is Amplified in All Resistant Anopheles gambiae Mosquitoes: Fitness Consequences of Homogeneous and Heterogeneous Duplications[END_REF][START_REF] Lucas | Whole-genome sequencing reveals high complexity of copy number variation at insecticide resistance loci in malaria mosquitoes[END_REF], which suggests they share a common evolutionary origin despite their multi-species distribution. Second, partial sequencing of Ace1 has shown that 280S alleles share highly similar haplotypic backgrounds in both West African A. gambiae and A. coluzzii [START_REF] Djogbénou | Evidence of Introgression of the ace-1R Mutation and of the ace-1 Duplication in West African Anopheles gambiae s[END_REF][START_REF] Essandoh | Acetylcholinesterase (Ace-1) target site mutation 119S is strongly diagnostic of carbamate and organophosphate resistance in Anopheles gambiae s.s. and Anopheles coluzzii across southern Ghana[END_REF]Weetman et al. 2015), which is suggestive of inter-specific introgression and a selective sweep. The most parsimonious hypothesis linking our results and these observations would posit that (i) the high similarity of 280S-carrying haplotypes across the A. gambiae -A. coluzzii species boundary is shared along the entire duplicated region (ca. 200 kbp); and (ii) this similarity is due to a 280S-linked selective sweep around Ace1, followed by duplications and the joint introgression of the duplication-280S haplotype between A. gambiae and A. coluzzii.

To test these hypotheses, we investigated genetic variation in the Ace1 duplication using networks of haplotype similarity (Figure 6). The analysis of Ace1 haplotypes is limited by the low density of phased variants in the vicinity of this gene (Supplementary Material SM9), caused by a combination of high sequence conservation and the difficulties in phasing haplotypes within a heterozygous duplication. To circumvent these limitations, we built four different networks around tagging variants that were in strong linkage disequilibrium with G280S (Huff and Rogers' r > 0.95) and located within the duplication (located within -26 kbp and +12 kbp from G280S; Supplementary Material SM9), and used these to define four sets of 280S-linked haplotypes. Then, we assumed that concordant observations among all sets of 280S-linked haplotypes would reflect genetic signatures shared with Ace1.

The minimum spanning tree built from haplotypes located around the first tagging variant (Figure 6A) identified a cluster of 124 identical haplotypes carrying the 280S-linked allele (Figure 6A, labelled with green text) and two larger haplotype clusters linked to the wt 280G allele (Figure 6A). identical haplotypes originating in the same A. gambiae and A. coluzzii specimens from West Africa (Figure 6B, Supplementary Material SM10). This result indicates that haplotype similarity between A. gambiae and A. coluzzii extends beyond the Ace1 gene.

We also investigated possible signals of positive selection around the Ace1 duplication breakpoints using Garud's H statistics, haplotypic diversity, and extended haplotype homozygosity (Figure 7, Supplementary Material SM11 and SM12). The profile of Garud's H 12 statistic in 280S-linked haplotypes showed peaks both upstream and downstream of the Ace1 duplication (Figure 7A), which coincided with a low H 2 /H 1 ratio (Figure 7B) and low haplotypic diversity (Figure 7C), and is thus indicative of a hard selective sweep in this region [START_REF] Garud | Recent Selective Sweeps in North American Drosophila melanogaster Show Signatures of Soft Sweeps[END_REF][START_REF] Messer | Population genomics of rapid adaptation by soft selective sweeps[END_REF]. Notably, the region of extended haplotype homozygosity associated with this sweep was still apparent upstream and downstream of the duplication breakpoints, i.e., far away from the focal tagging variant used to define 280S-linked haplotypes (Figure 7D). Statistics are calculated in blocks of 500 variants with 20% block overlap. Includes averages of each statistic upstream and downstream of the breakpoints (with standard errors from jackknife haplotype resampling). D) Extended haplotype homozygosity at the duplication breakpoints for haplotypes carrying the 280S-or wttagging variants at the 2R:3481632 locus. Additional plots for all tagging variants are available in Supplementary Material SM11 and SM12.

Next, we tested whether 280S alleles and the duplication spread jointly between A. coluzzii and A. gambiae via introgression (Figure 8). We examined the chromosome-wide profile of introgression between A. gambiae and A. coluzzii populations using Patterson's D statistic [START_REF] Durand | Testing for Ancient Admixture between Closely Related Populations[END_REF][START_REF] Patterson | Ancient Admixture in Human History[END_REF], which compares allele frequencies between three putatively admixing populations (A, B and C) and one outgroup (O); and can identify introgression between populations A and C (D > 0) or B and C (D < 0). Specifically, we tested whether duplicated A. coluzzii from West African populations (population A) introgressed with A. gambiae populations from West and Central Africa (C), using wt Angolan A. coluzzii as a control (B) and A. arabiensis as an outgroup (O).

We evaluated whether the Ace1 duplication introgressed between A. gambiae and A. coluzzii in two scenarios: (i) in specimens with CNVs from populations with Ace1 resistance mutations (Figure 8A), and (ii) in non-duplicated specimens from the same populations (Figure 8B). We only find evidence of introgression in the first case, i.e. between West African A. gambiae and A. coluzzii specimens carrying Ace1 duplications (Figure 8A). For example, we found evidence of introgression between duplicated Ghanaian A. gambiae and duplicated Ivorian A. coluzzii (D = 0.352 +/-0.040, p = 3.1 × 10 -18 from a standard distribution of block-jackknifed estimates; green lines in Figure 8A); but not with non-duplicated Ivorian A. coluzzii (D = -0.026 +/-0.043, p = 0.54; green lines in Figure 8B). In all comparisons where introgression was apparent, the genomic region of increased D values extended along the entire duplicated region. On the other hand, there was no evidence of introgression in any comparison involving non-duplicated A. gambiae (D ≈ 0; purple lines in Figure 8A and 8B).

To establish the direction of introgression, we performed a phylogenomic analysis of haplotypes from the duplicated region (Figure 9). Firstly, haplotypes with Ace1 duplications formed a single clade containing A. gambiae and A. coluzzii sequences to the exclusion of all non-duplicated sequences from both species. Secondly, we found that duplicated sequences were phylogenetically closer to the wt A. gambiae clade than to wt A. coluzzii (Figure 9A). Thirdly, an analysis of differences in allelic frequencies between wt A. gambiae, wt A. coluzzii and duplicated haplotypes also indicated that the duplicated haplotypes are more similar to wt A. gambiae than to wt A. coluzzii: the branch leading to wt A. coluzzii since the divergence from the duplicated specimens was longer (L = 0.0273 ± 0.0015 standard error) than the branch leading to wt A.

gambiae (L = 0.0102 ± 0.0011). This topology indicates that duplicated haplotypes emerged from a A. gambiae wt background and later introgressed into A. coluzzii. In contrast, phylogenies built from non-introgressed regions upstream and downstream of the duplication exhibited a topology with species-specific clades that included both duplicated and non-duplicated specimens (Figure 9B,C).

Altogether, these results show that the 280S resistance alleles present across West Africa appear in similar genetic backgrounds irrespectively of the sampled population and species, and that this similarity extends beyond the Ace1 gene to encompass the entire duplicated region (Figure 6). This low haplotypic diversity is due to a hard selective sweep, detectable at the duplication breakpoints (Figure 7). This resistance haplotype, which spans ca. 200 kbp and includes Ace1 and ten other genes (Supplementary Material SM3) emerged in A. gambiae and later spread to A. coluzzii (Figures 8 and9).

Discussion

Evolutionary history of G280S and Ace1 duplications

Our results show that pirimiphos-methyl resistance is associated with a combination of two mutations in Ace1: the G280S SNP and CNVs. In West Africa, virtually all CNVs are found as duplications of a wide region (ca. 200 kbp) that includes Ace1 and 10 other genes [START_REF] Assogba | The ace-1 Locus Is Amplified in All Resistant Anopheles gambiae Mosquitoes: Fitness Consequences of Homogeneous and Heterogeneous Duplications[END_REF][START_REF] Lucas | Whole-genome sequencing reveals high complexity of copy number variation at insecticide resistance loci in malaria mosquitoes[END_REF]. This duplication has a unique evolutionary origin in the A. gambiae/A. coluzzii species pair, as its breakpoints are consistent in all populations from both species studied so far [START_REF] Assogba | The ace-1 Locus Is Amplified in All Resistant Anopheles gambiae Mosquitoes: Fitness Consequences of Homogeneous and Heterogeneous Duplications[END_REF][START_REF] Lucas | Whole-genome sequencing reveals high complexity of copy number variation at insecticide resistance loci in malaria mosquitoes[END_REF]. The distribution of duplication-280S resistance haplotypes observed in the 1000 Genomes cohort can explained by three evolutionary events: the G280S mutation, an internally heterozygous duplication, and interspecific introgression. Furthermore, recent sampling efforts have identified internal deletions within the duplication in both A. gambiae and A. coluzzii [START_REF] Assogba | Adaptive deletion in resistance gene duplications in the malaria vector Anopheles gambiae[END_REF]. Homogeneous duplications are rare in the 1000 Genomes dataset but appear to be quite common in some A. gambiae populations, though less so in A. coluzzii (Table 2; (Weetman et al. 2015;[START_REF] Assogba | Adaptive deletion in resistance gene duplications in the malaria vector Anopheles gambiae[END_REF].

At present, we cannot establish the relative order of occurrence of the SNP and CNV in Ace1 because these mutations are tightly linked in the 1000 Genomes dataset (Figure 2A). Nevertheless, the detection of wt homozygous duplications in samples collected in 2002 in southern Ghana (Weetman et al. 2015)-a country where 280S alleles have more recently been found at high frequencies (Figure 1)-raises the possibility that the SNP might have occurred on an already duplicated background. This order of events would have initially resulted in permanent heterozygosity, providing resistance and reducing the fitness costs associated with impaired acetylcholinesterase activity in the absence of insecticide [START_REF] Labbé | Independent Duplications of the Acetylcholinesterase Gene Conferring Insecticide Resistance in the Mosquito Culex pipiens[END_REF][START_REF] Assogba | An ace-1 gene duplication resorbs the fitness cost associated with resistance in Anopheles gambiae, the main malaria mosquito[END_REF]. If that were the case, the emergence of 280S-homozygous duplications would require either a secondary loss of a wt copy or a parallel G280S mutation in the initially heterozygous duplication.

After these two mutations, the joint duplication-280S resistance haplotype introgressed from A. gambiae into A. coluzzii (Figures 8 and9). This A. gambiae origin is consistent with earlier studies of Ace1 variation in West African locations, which generally reported higher frequencies of 280S [START_REF] Djogbénou | Evidence of Introgression of the ace-1R Mutation and of the ace-1 Duplication in West African Anopheles gambiae s[END_REF][START_REF] Clarkson | The genetic architecture of target-site resistance to pyrethroid insecticides in the African malaria vectors Anopheles gambiae and Anopheles coluzzii[END_REF][START_REF] Essandoh | Acetylcholinesterase (Ace-1) target site mutation 119S is strongly diagnostic of carbamate and organophosphate resistance in Anopheles gambiae s.s. and Anopheles coluzzii across southern Ghana[END_REF] and CNVs [START_REF] Assogba | Adaptive deletion in resistance gene duplications in the malaria vector Anopheles gambiae[END_REF] in A. gambiae than in A. coluzzii, as one would expect if they had a longer evolutionary history in the former. Whilst previous studies had suggested that the similarity of 280S haplotypes in A. coluzzii and A. gambiae was due to inter-specific introgression [START_REF] Djogbénou | Evidence of Introgression of the ace-1R Mutation and of the ace-1 Duplication in West African Anopheles gambiae s[END_REF][START_REF] Essandoh | Acetylcholinesterase (Ace-1) target site mutation 119S is strongly diagnostic of carbamate and organophosphate resistance in Anopheles gambiae s.s. and Anopheles coluzzii across southern Ghana[END_REF]Weetman et al. 2015), they were focused on partial sequencing of the Ace1 gene and could not assess the relationship between introgression and the wider duplicated region (ca. 30 times longer than the Ace1 coding region). Furthermore, analyses focused on the Ace1 gene body can be confounded by the lack of sequence variation in this region (Supplementary SM9). By leveraging variation in linkage disequilibrium beyond the gene, we are able identify clear signatures of introgression along the entire Ace1 duplication region (Figure 8).

Genetic basis of pirimiphos-methyl resistance in A. coluzzii and A. gambiae

We have uncovered a strong association between resistance to pirimiphos-methyl and the possession of two or more copies of the Ace1 280S allele in A. coluzzii from CoB te d'Ivoire (Figure 3). While the G280S mutation alone has been previously linked to resistance to carbamates and organophosphates in multiple insects [START_REF] Weill | Insecticide resistance in mosquito vectors[END_REF][START_REF] Weill | The unique mutation in ace-1 giving high insecticide resistance is easily detectable in mosquito vectors[END_REF][START_REF] Labbé | Independent Duplications of the Acetylcholinesterase Gene Conferring Insecticide Resistance in the Mosquito Culex pipiens[END_REF][START_REF] Alout | Different Amino-Acid Substitutions Confer Insecticide Resistance Through Acetylcholinesterase 1 Insensitivity in Culex vishnui and Culex tritaeniorhynchus (Diptera: Culicidae) from China[END_REF][START_REF] Djogbénou | Evidence of Introgression of the ace-1R Mutation and of the ace-1 Duplication in West African Anopheles gambiae s[END_REF][START_REF] Dabiré | Distribution and Frequency of kdr Mutations within Anopheles gambiae s.l. Populations and First Report of the Ace.1G119S Mutation in Anopheles arabiensis from Burkina Faso (West Africa)[END_REF][START_REF] Liebman | Novel mutations on the ace-1 gene of the malaria vector Anopheles albimanus provide evidence for balancing selection in an area of high insecticide resistance in Peru[END_REF][START_REF] Feyereisen | Genotype to phenotype, the molecular and physiological dimensions of resistance in arthropods[END_REF], our results indicate that resistance to pirimiphos-methyl specifically relies on the presence of multiple 280S alleles.

Across the 1000 Genomes cohort, specimens with 280S alleles can be grouped into two categories according to the number of Ace1 copies (Figure 2C): (i) heterozygous duplications with multiple 280S alleles and/or multiple wt, and (ii) high-copy, 280S-homozygous specimens. Heterozygous duplications are the most frequent combination of Ace1 mutations (113 out of 119 samples, Figure 2B) and the only one identified in the A. coluzzii population from CoB te d'Ivoire (Tiassalé). Therefore, we decided to examine their contribution to pirimiphos-methyl resistance in a wider array of A. gambiae and A. coluzzii populations from Ghana, CoB te d'Ivoire, and Togo (Tables 1 and2). We identified a common theme in most populations surveyed here: (i) 280S heterozygotes were the most common resistance genotype in six out of eight populations analysed, including all A. coluzzii sampling locations and three out of five in A. gambiae (Table 1); and (ii) among the heterozygotes, we found that the ratio of 280S alleles and Ace1 copy number were significantly associated with resistance in five out of seven tested populations (Table 2). This provides a wider demonstration that the combined effect of G280S and Ace1 duplications is significantly associated with pirimiphos-methyl resistance across multiple West African populations of both species; consistent with a requirement for multiple copies of 280S for resistance.

We also identified six high-copy, 280S-homozygous genomes sampled from the Ghanaian A. gambiae population in the 1000 Genomes dataset (Figures 1 and2). These samples had not been assayed for pirimiphos-methyl resistance prior to sequencing. However, our extended genotype-phenotype analysis in eight West African populations (Table 1) showed that (i) 280S homozygotes exhibited increased resistance to pirimiphos-methyl, and that, in fact, they were the most abundant resistance allele in two A.

gambiae populations (Baguida and Obuasi). It is worth mentioning that survival rates among 280S homozygotes were noticeably lower these two locations (56.8% -58.9%) than in populations where heterozygotes were more common (> 95%; Table 1). Our analyses cannot fully capture the resistance landscape in these populations, because we have focused on the role of CNVs in heterozygous specimens (as they were the most abundant resistant type in the 1000 Genomes cohort; Figure 1). Thus, we cannot presently explain these differences, which may reflect an effect of different phenotyping strategies employed in each population (see Methods), or additional resistance mechanisms. A detailed examination of 280S-homozygous duplications-e.g. effects on gene dosage or fitness-could shed light onto the phenotypic variability observed in populations with high penetrance of Ace1 mutations.

Finally, we have performed a genome-wide scan in a population of A. coluzzii from CoB te d'Ivoire (Tiassalé) to identify loci associated with pirimiphos-methyl resistance (Figure 5). Our investigation of the association between k-mer frequencies and resistance confirmed the overwhelming concentration of phenotypic association in the Ace1 duplication (Figure 5B), supporting the conclusion that it is the primary driver of resistance in this population. The correlations between Ace1 duplication and all but one of the kmers significantly associated with resistance, but mapping elsewhere in the genome strongly suggests non-independence of the signals, but the exact cause of this is unclear. Whilst we cannot entirely discount linkage disequilibrium arising from epistasis, a simpler and parsimonious explanation would be that this result is caused by non-specific or incorrect sequence alignments, possibly linked to currently unrecognised variation in the Ace1 duplication.

The lack of clear signals of pirimiphos-methyl adaptation other than Ace1 in this population may reflect its collection at an early stage of the development of pirimiphos-methyl resistance in West Africa (Edi et al. 2014b). Therefore, our genome-wide analyses would only capture pre-existing variants that were enriched in the resistant specimens. Ace1 mutations and duplications, which cause cross-resistance with previously employed insecticides [START_REF] Djogbénou | Estimation of allele-specific Ace-1 duplication in insecticide-resistant Anopheles mosquitoes from West Africa[END_REF][START_REF] Djogbénou | Evidence of Introgression of the ace-1R Mutation and of the ace-1 Duplication in West African Anopheles gambiae s[END_REF][START_REF] Essandoh | Acetylcholinesterase (Ace-1) target site mutation 119S is strongly diagnostic of carbamate and organophosphate resistance in Anopheles gambiae s.s. and Anopheles coluzzii across southern Ghana[END_REF]Edi et al. 2014a;[START_REF] Assogba | An ace-1 gene duplication resorbs the fitness cost associated with resistance in Anopheles gambiae, the main malaria mosquito[END_REF][START_REF] Assogba | The ace-1 Locus Is Amplified in All Resistant Anopheles gambiae Mosquitoes: Fitness Consequences of Homogeneous and Heterogeneous Duplications[END_REF][START_REF] Assogba | Adaptive deletion in resistance gene duplications in the malaria vector Anopheles gambiae[END_REF]Weetman et al. 2015), fit this definition and are sufficient to sustain high levels of resistance. However, future analyses on recent collections from populations subjected to regular pirimiphos-methyl-based IRS may reveal additional mechanisms more specifically selected by this insecticide.

Implications for insecticide intervention programmes

Our results demonstrate that the duplication-280S haplotype represents a powerful marker to diagnose resistance to pirimiphos-methyl in A. coluzzii and A. gambiae, permitting the spread of resistance to be tracked from any preserved sample from which DNA can be extracted. Our initial survey of an A. coluzzii population from Tiassalé (CoB te d'Ivoire) shows that discrimination between susceptible and resistant specimens can be attained with high accuracy in this population according to the number of 280S alleles.

Specifically, following the resistance thresholds apparent in Figure 3B, we can (i) classify wt homozygous samples as susceptible (100% predictive value) and (ii) separate heterozygous individuals between those with more wt than 280S alleles (susceptible) and those with equal or more 280S than wt alleles (resistant).

This would yield a positive predictive value of 77% for resistance, and a negative predictive value of 85% for susceptibility.

In practice, precise quantification of the number of resistance 280S alleles will often be difficult.

Nevertheless, we show that key variation in Ace1 can be captured using two qPCR genotyping assays that predict pirimiphos-methyl resistance in populations of A. gambiae and A. coluzzii (Table 2): (i) a measurement of the ratio of 280S to 280G alleles [START_REF] Bass | Development of high-throughput real-time PCR assays for the identification of insensitive acetylcholinesterase (ace-1R) in Anopheles gambiae[END_REF], and (ii) the number of Ace1 copies (Edi et al. 2014a). Crucially, measurements of the 280S allelic ratio can be obtained from the same assay that is currently used for standard genotyping of the G280S mutation [START_REF] Bass | Development of high-throughput real-time PCR assays for the identification of insensitive acetylcholinesterase (ace-1R) in Anopheles gambiae[END_REF], which thus has good diagnostic capacity in homozygotic and heterozygotic specimens. The phenotypic assays we employed were intended to provide more accurate separation than standard dose approaches [START_REF] Weetman | Evolution of insecticide resistance diagnostics in malaria vectors[END_REF], but the strong effect of Ace1 CNVs and G280S mutations on pirimiphos-methyl resistance is sufficiently robust for resistance diagnosis in single time-dose discriminating assays as well -as used in the Ivorian A. coluzzii specimens from the 1000 Genomes cohort.

Furthermore, monitoring studies must note that this resistance haplotype is still evolving, especially in the light of the fitness trade-offs associated with Ace1 resistance alleles -namely, that heterozygous duplications offset the fitness cost of 280S alleles in the absence of insecticide exposure [START_REF] Labbé | Independent Duplications of the Acetylcholinesterase Gene Conferring Insecticide Resistance in the Mosquito Culex pipiens[END_REF][START_REF] Assogba | An ace-1 gene duplication resorbs the fitness cost associated with resistance in Anopheles gambiae, the main malaria mosquito[END_REF], and that variably copy number can result in gene dosage disturbances [START_REF] Assogba | The ace-1 Locus Is Amplified in All Resistant Anopheles gambiae Mosquitoes: Fitness Consequences of Homogeneous and Heterogeneous Duplications[END_REF]. A recent survey has shown that internal deletions within the duplicated region, downstream to Ace1, are spreading in both A. gambiae and A. coluzzii [START_REF] Assogba | Adaptive deletion in resistance gene duplications in the malaria vector Anopheles gambiae[END_REF]). The authors have proposed that these deletions reduce the fitness costs of duplications by ameliorating deleterious changes in gene expression [START_REF] Assogba | Adaptive deletion in resistance gene duplications in the malaria vector Anopheles gambiae[END_REF]. These small deletions are absent from the Ivorian A. coluzzii genomes analysed here [START_REF] Lucas | Whole-genome sequencing reveals high complexity of copy number variation at insecticide resistance loci in malaria mosquitoes[END_REF]), but future genome monitoring studies could investigate their spread and identify signals of selection associated with their proposed selective advantage. Likewise, the independently evolved Ace1 CNV that we have identified in one A. gambiae specimen from Guinea (Figure 2B, Supplementary Material SM4) might merit further attention in the future. Finally, as noted above, we cannot rule out the possibility that, following several years of widespread use of pirimiphos-methyl for IRS, additional non-target site-based resistance mechanisms might have emerged.

The duplication-280S haplotype pre-dates the widespread adoption of pirimiphos-methyl in 2013 [START_REF] Oxborough | Trends in US President's Malaria Initiative-funded indoor residual spray coverage and insecticide choice in sub-Saharan Africa (2008-2015): urgent need for affordable, long-lasting insecticides[END_REF], which indicates that it likely spread due to its adaptive value with respect to previously employed insecticides used in disease control or agriculture, such as carbamates (bendiocarb (Edi et al. 2014a;[START_REF] Assogba | The ace-1 Locus Is Amplified in All Resistant Anopheles gambiae Mosquitoes: Fitness Consequences of Homogeneous and Heterogeneous Duplications[END_REF])) or other organophosphates (fenitrothion, chlorpyrifos-methyl [START_REF] Essandoh | Acetylcholinesterase (Ace-1) target site mutation 119S is strongly diagnostic of carbamate and organophosphate resistance in Anopheles gambiae s.s. and Anopheles coluzzii across southern Ghana[END_REF][START_REF] Assogba | The ace-1 Locus Is Amplified in All Resistant Anopheles gambiae Mosquitoes: Fitness Consequences of Homogeneous and Heterogeneous Duplications[END_REF]). In that regard, we show that the evolution of the pirimiphos-methyl resistance haplotype is deeply intertwined with that of previous insecticide adaptations, and reveal how continued genomic monitoring studies can help us understand the influence of previous interventions on future population control efforts. Furthermore, whilst current efficacy of pirimiphos-methyl in Actellic IRS formulations appears to be retained [START_REF] Sherrard-Smith | Systematic review of indoor residual spray efficacy and effectiveness against Plasmodium falciparum in Africa[END_REF], its potential could be limited in the Anopheles populations where resistance to these insecticides is already common, most notably in West Africa (Figures 1 and2). This emphasises the crucial importance of now-commencing resistance management strategies using additional IRS insecticides with alternative modes of action [START_REF] Oxborough | Susceptibility testing of Anopheles malaria vectors with the neonicotinoid insecticide clothianidin; Results from 16 African countries, in preparation for indoor residual spraying with new insecticide formulations[END_REF].

Materials and methods

Data collection

We used genomic variation data from individual mosquitoes, obtained from the Phase 2 (AR1) of the Anopheles gambiae 1000 Genomes project, as described in (Miles et al. 2017 [START_REF] Delaneau | Haplotype estimation using sequencing reads[END_REF], and variant effects were predicted using SnpEff 4.1b [START_REF] Cingolani | A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3[END_REF]. We downloaded the genotype calls, SNP effect predictions, and haplotype phasing as available in Anopheles gambiae 1000 Genomes Phase 2 online archive (The Anopheles gambiae 1000 Genomes Consortium 2017).

The catalog of Ace1 CNVs in the 1000 Genomes dataset (Phase 2) was obtained from [START_REF] Lucas | Whole-genome sequencing reveals high complexity of copy number variation at insecticide resistance loci in malaria mosquitoes[END_REF].

There, the authors used calculated the coverage (sequencing depth) of each whole-genome sequenced sample in non-overlapping windows along the genome (300 bp) and normalised this value so as to obtain an expected average value of coverage = 2 in non-duplicated, diploid regions. A Gaussian HMM procedure was then applied to the normalised windowed coverage data so as to call windows with heightened normalised coverage (>2). A detailed account of these methods can be found in the original publication [START_REF] Lucas | Whole-genome sequencing reveals high complexity of copy number variation at insecticide resistance loci in malaria mosquitoes[END_REF], and the results for the Ace1 region are available in Supplementary Material SM4.

We retrieved the reference genome information for A. gambiae from VectorBase (Giraldo-Calderón et al. coluzzii based on PCR assay [START_REF] Fanello | Simultaneous identification of species and molecular forms of the Anopheles gambiae complex by PCR-RFLP[END_REF]. The 71 larvae were tested for resistance to pirimiphosmethyl resistance using a WHO tube assay with 0.25% impregnated papers (World Health Organization 2018b), which led to the identification of 31 resistant (live) and 40 susceptible (dead) specimens.

We tested the association of these resistance phenotypes with various genetic variants present in Ace1 in this population. This includes the non-synonymous mutations G280S and S65A, the number of Ace1 copies, and the number of 280S alleles in each individual sample (listed in Supplementary Material SM6A).

For the non-synonymous mutations, we assessed genotype-phenotype associations for each variant independently using Fisher's exact test (gmodels R library (Warnes et al. 2018)) and estimated odds ratio and 95% confidence intervals using the Woolf method with Haldane-Anscombe correction (Prop.or function in the pairwiseCI R library [START_REF] Schaarschmidt | pairwiseCI: Confidence Intervals for Two Sample Comparisons[END_REF]). For the Ace1 and 280S copy numbers, we used used binomial generalised linear models (glm function in R stats library, family binomial) to obtain odds ratio estimates for each of the four variables.

We also built a binomial generalised linear models with all four variables (G280S and A65S genotypes;

Ace1 and 280S copy number) and a step-wise procedure to remove non-informative variants from the model (step function in R stats using k = log(n) as the degree of freedom, as required by the Bayesian Information Criterion; n represents the number of samples in our dataset). The number of Ace1 copies and the number of mutated alleles were encoded as continuous variables, and the genotypes of G280S and S65A

were encoded as categorical variables. The performance of all generalised linear models was assessed relative to a null model (intercept as the only variable) using a χ 2 test in an analysis of variance (anova function in the R stats library). A detailed statistical analysis of all comparisons mentioned above can be found in Supplementary Material SM6A and SM6B.

The number of Ace1 copies in each genome of Anopheles gambiae 1000 Genomes Phase 2 was retrieved from [START_REF] Lucas | Whole-genome sequencing reveals high complexity of copy number variation at insecticide resistance loci in malaria mosquitoes[END_REF]. In that study, copy number states of multiple CNVs were inferred using a HMMbased predictive model that took as input the normalised sequencing depth calculated along nonoverlapping 300 bp genomic windows.

The estimated number of 280S alleles in each genome was estimated in the following manner: (i) we retrieved the copy numbers of Ace1 in each genome (C; see above); (ii) we calculated the fraction of reads supporting the reference and alternate alleles; and (iii) assigned the number of 280S alleles (S) as the value that minimised the difference between the fraction of alternate alleles and S/C, for all discrete S values between 0 and C. For example, a genome with three Ace1 copies (C = 3) and 30% of reads supporting 280S would carry one 280S allele (S = 1), given that S/C = 1/3 ≈ 30%.

Genotype-phenotype association in additional West African populations

We collected 1080 female specimens of two species (A. coluzzii and A. gambiae) from six different locations across West Africa (CoB te d'Ivoire: Aboisso; Ghana: Madina, Korle-Bu, Weija and Obuasi; Togo: Baguida). Species identity was determined using two methods designed to discriminate between A.

gambiae, A. coluzzii and A. arabiensis: a PCR of species-specific SINE insertion polymorphisms as described in [START_REF] Santolamazza | Insertion polymorphisms of SINE200 retrotransposons within speciation islands of Anopheles gambiae molecular forms[END_REF], and a high-resolution melt curve analysis [START_REF] Chabi | Rapid high throughput SYBR green assay for identifying the malaria vectors Anopheles arabiensis, Anopheles coluzzii and Anopheles gambiae s.s. Giles[END_REF] For each sample, we determined the G280S genotype using a qPCR TaqMan assay as described by [START_REF] Bass | Development of high-throughput real-time PCR assays for the identification of insensitive acetylcholinesterase (ace-1R) in Anopheles gambiae[END_REF]) (list of samples in Supplementary Material SM6C), in which 280S-and 280G-specific primers where labeled with different fluorescent dyes (FAM and HEX, respectively). We also calculated the ratio of FAM-to-HEX fluorescent dye signal in hetrerozygotes, which label 280S and 280G alleles respectively, as an index of the fraction of 280S allele copies (Edi et al. 2014a;[START_REF] Djogbénou | Estimation of allele-specific Ace-1 duplication in insecticide-resistant Anopheles mosquitoes from West Africa[END_REF].

Detailed methods for these two genotyping assays are available in Supplementary Material SM6G.

Introgression analysis

We used Patterson's D statistic [START_REF] Durand | Testing for Ancient Admixture between Closely Related Populations[END_REF][START_REF] Patterson | Ancient Admixture in Human History[END_REF]) to detect introgression between A. coluzzii and A. gambiae. We retrieved the allele frequencies of all biallelic, non-singleton variants from chromosome arm 2R that were segregating in West African populations where we had identified Ace1 duplications (namely: 75 A. coluzzii from Burkina Faso, 55 from Ghana, and 71 from CoB te d'Ivoire; and 92 A. gambiae from Burkina Faso, 12 from Ghana and 40 from Guinea; total n = 345 genomes) [START_REF] Lucas | Whole-genome sequencing reveals high complexity of copy number variation at insecticide resistance loci in malaria mosquitoes[END_REF], as well as Central African populations that we used as non-admixed controls (78 A. For each combination of (((A,B),C),O) populations, we calculated the average D statistic within the duplicated (position 3,436,800 to 3,639,600, spanning 202.8 kbp around Ace1), and estimated its deviation from the null expectation (no introgression: D = 0) with a block-jackknife procedure (block length = 100 variants; average_patterson_d function in scikit-allel), which we used to estimate the standard error, Z-score and the corresponding p value from the two-sided normal distribution (Supplementary Material SM13).

We assessed the degree of genetic differentiation between resistant and susceptible A. coluzzii from CoB te d'Ivoire. We focused on non-singleton variants that were segregating in this population. In total, 8,431,869 out of 57,837,885 (14.57%) variants from chromosomes 2, 3 and X were retained for further analysis.

For each chromosome arm (2R, 2L, 3R and X), we calculated genotype counts in each subpopulation (resistant and susceptible), and calculated their genetic differentiation using the Hudson's F ST statistic [START_REF] Hudson | Estimation of levels of gene flow from DNA sequence data[END_REF][START_REF] Bhatia | Estimating and interpreting FST: The impact of rare variants[END_REF]) along non-overlapping windows of 1,000 variants (moving_hudson_fst function in the scikit-allel v1.1.10 library (Miles and Harding 2017) from Python 3.4). We also calculated the average F ST in each chromosomal arm, with standard errors obtained from a block-jackknife procedure (using non-overlapping windows of 1,000 variants; average_hudson_fst function).

We calculated the normalised population branching statistic (PBS) [START_REF] Yi | Sequencing of 50 Human Exomes Reveals Adaptation to High Altitude[END_REF][START_REF] Malaspinas | A genomic history of Aboriginal Australia[END_REF] in non-overlapping windows of 1,000 variants along each chromosomal arm (pbs function in scikit-allel).

using resistant and susceptible Ivorian A. coluzzii as test populations, and A. coluzzii from Angola as an outgroup. Angolan A. coluzzii (n = 78) where selected as outgroup due to their relative isolation relative to West African A. coluzzii populations (The Anopheles gambiae 1000 Genomes Consortium 2019; Miles et al. 2017) and their putatively naive profile of organophosphate resistance [START_REF] Oxborough | Trends in US President's Malaria Initiative-funded indoor residual spray coverage and insecticide choice in sub-Saharan Africa (2008-2015): urgent need for affordable, long-lasting insecticides[END_REF]World Health Organization 2018a). The distribution of PBS estimates along each chromosomal arm was standardised to unit variance (standardize function in scikit-allel), and the resulting distribution of Z-scores was used to derive a two-sided p-value that reflected extreme values of PBS (highly positive or negative). We corrected for multiple testing using local estimation of false discovery rates (fdrtool function in the fdrtool 1.2.15 R library [START_REF] Klaus | fdrtool: Estimation of (Local) False Discovery Rates and Higher Criticism[END_REF]). Finally, we selected genomic windows with high differentiation in resistant specimens (standardised PBS > 0, significance threshold FDR < 0.001).

The estimates of PBS and Hudson's F ST in each genomic window are available as Supplementary Material SM7.

We performed a principal component analysis of resistant and susceptible Ivorian A. coluzzii. We first obtained a set of genetically unlinked variants from chromosomal arms 3R and 3L (so as to avoid the confounding effects of chromosomal inversions in arms 2R, 2L and X). We discarded linked variants within 500 bp consecutive windows (using a 200 bp step), using a Rogers and Huff's r threshold value of 0.1

  6

Figure 1 .

 1 Figure 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 peptideand 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.

Figure 2 .

 2 Figure2. 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 homozygotes (purple), duplicated heterozygotes (orange) and duplicated 280S homozygotes (green). B) Breakdown of of possible 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) Summary of haplotypes observed in the Ace1 locus.

Figure 4 .

 4 Figure 4. Principal component analysis of CoB te d'Ivoire A. coluzzii. PCA constructed using genotypes of 791 unlinked variants from chromosome 3.

Figure 5 .

 5 Figure5. 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.

Figure 6 .

 6 Figure6. 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 Supplementary Material SM10). 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 four tagging variants (Supplementary Material SM10).

Figure 7 .

 7 Figure 7. Positive selection around the Ace1 duplication. A-C) Profile of Garud's H 12 , H 2 /H 1 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 upstream and downstream of the breakpoints (with standard errors from jackknife haplotype resampling). D) Extended haplotype homozygosity at the duplication breakpoints for haplotypes carrying the 280S-or wttagging variants at the 2R:3481632 locus. Additional plots for all tagging variants are available in Supplementary Material SM11 and SM12.

Ace1B)Figure 8 .

 8 Figure 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 CoB 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 CoB 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 CoB te d'Ivoire as population A. Detailed lists of all statistical tests and replicate analyses with additional populations are available in Supplementary SM13. Species abbreviations: gam, A. gambiae; col, A. coluzzii; ara, A. arabiensis. Country abbreviations: AO, Angola; CI, CoB te d'Ivoire; GH, Ghana.

  2015), including the genome assembly (PEST, AgamP4 version), gene annotation coordinates in GFF format (AgamP4.9) and gene names and descriptions.Genotype-phenotype association in CoB te d'Ivoire A. coluzzii genomesAs part of Anopheles gambiae 1000 Genomes Phase 2, 71 A. coluzzii specimens from CoB te d'Ivoire were collected and phenotyped for pirimiphos-methyl resistance before whole-genome sequencing (The Anopheles gambiae 1000 Genomes Consortium 2019). These samples were collected as larvae in rice fields near 5.898) between May and September 2012, and were identified as A.

  coluzzii from Angola; 69 A. gambiae from Gabon and 297 from Cameroon). In addition, we retrieved the genotypes at the same variant positions for the populations of A. arabiensis (12 genomes), A. quadriannulatus (10), A. merus (10) and A. melas (4) populations analysed in[START_REF] Fontaine | Extensive introgression in a malaria vector species complex revealed by phylogenomics[END_REF], which we used as outgroups in the calculation of Patterson's D. In total, we retained ca. 12% of the 47,817,813 variants in 2R in each analysis.To determine whether duplicated and non-duplicated populations had different introgression patterns in Ace1, we calculated windowed means of Patterson's D (length 5,000 variants, 20% step; moving_patterson_d in scikit-allel) using various combinations of populations with and without the Ace1 duplication (Figure5A). This test requires genotype frequencies in four populations (two ingroups A and B; one candidate donor/receptor C with which A or B might have introgressed; and one unadmixed outgroup O) branching in a (((A,B),C),O) topology (Figure5A). We observed the following convention: (i) we used West African A. coluzzii populations as A, discriminating between duplicated and non-duplicated subpopulations; (ii) we used non-duplicated A. coluzzii from Angola as control B; (iii) we used various A. gambiae populations as C, again analysing duplicated and non-duplicated specimens within each population separately; and (iv) A. arabiensis, A. quadriannulatus, A. merus and A. melas as outgroup O.

  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 Supplementary Material SM6.

				). Crucially, within the group of duplicated specimens with one 280S
	A) Ace1 & 280S copies		B) Phenotypic association
		4	susceptible resistant		no yes 280S presence
						resistant		31
		3				susceptible	9	31
	280S copies	2				resistant survival:	2 0%	5 3 Ace1 copies 22 4 50%	4 5
		1				susceptible	9	24	7
						survival:	13.2%	78.8%
		0					0	280S copies 1 2 3	4
			2	Ace1 copies 3 4	5	resistant susceptible	9	4 23	18 7	6 1	3
						survival:	11.1%	77.1%
	Figure 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
					8

Table 1 .

 1 Association of resistance genotypes with pirimiphos-methyl resistance in nine West African populations of A. coluzzii and A. gambiae. For each population, we report number of sampled specimens (n), the number of wt homozygous (wt), heterozygous (het) and 280S homozygous specimens (hom), survival frequency among mutated groups, and the odds ratios of survival in heterozygous and 280S homozygous specimens (if available) in a binomial generalised linear model, and the p-value of this model in a χ 2 ANOVA comparison with a null model. Odds ratios (OR) are reported with 95% confidence intervals. Dagger signs ( †) denote survival rates in genotypes with less than 5 specimens. Samples listed in Supplementary Material SM6C; statistical analysis in Supplementary Material SM6D.

	Country abbreviations: CI, CoB te d'Ivoire; GH, Ghana; TG, Togo.
	Population	n	Frequencies wt / het / hom	% survival het	OR het	% survival hom	OR hom	p
	coluzzii Aboisso (CI)	55	42 / 13 / 0	46.1%	35.1 (3.6 -338.0)	-	-	1.3 × 10 -4
	coluzzii Korle-Bu (GH)	214	25 / 165 / 24	49.6%	4.2 × 10 7 (0 -Inf)	95.8%	9.8 × 10 8 (0 -Inf)	1.2 × 10 -13
	coluzzii Weija (GH)	131	60 / 70 / 1	54.1%	7.3 (3.0 -17.6)	100% ( †)	3.7 × 10 7 (0 -Inf)	2.4 × 10 -6
	gambiae Aboisso (CI)	82	3 / 45 / 34	77.8%	7.0 (0.6 -85.4)	97%	66.0 (2.9 -1491.2)	3.7 × 10 -3
	gambiae Madina (GH)	172	19 / 112 / 41	58.9%	12.2 (2.7 -55.4)	100%	9.8 × 10 8 (0 -Inf)	4.4 × 10 -14
	gambiae Obuasi (GH)	140	0 / 28 / 112	17.8%	-	58.9%	6.6 (2.3 -18.6)	6.0 × 10 -5
	gambiae Weija (GH)	13	7 / 4 / 2	75% ( †)	7.3 × 10 17 (0 -Inf)	100% ( †)	2.6 × 10 9 (0 -Inf)	1.6 × 10 -3
	gambiae Baguida (TG)	102	0 / 14 / 88	57.1%	-	56.8%	0.9 (0.3 -3.1)	0.982

  Sample phenotypes and genotypes listed in Supplementary Material SM6E; statistical analysis in Supplementary Material SM6F. Country abbreviations: CI, CoB te d'Ivoire; GH, Ghana; TG, Togo.

	Population	n	GLM Ace1 copies (single variable) OR p	GLM 280S/280G allelic ratio (single variable) OR p	OR Ace1 copies	Minimal GLM (BIC) OR allelic ratio	p
	coluzzii Aboisso (CI)	11	1.3 (0.7 -2.3)	0.361	52.9 (0.01 -2.2 × 10 5

coluzzii dup gambiae no dup gambiae dup coluzzii no dup

	A) Ace1 duplication	B) Upstream to duplication
	coluzzii dup gambiae dup	coluzzii dup coluzzii no dup
	dup cluster: 97%
	dup within gam: 80%
		gambiae dup gambiae no dup
		0.05
		C) Downstream to duplication
	no dup gambiae	coluzzii no dup
	0.1
	A. coluzzii, no duplication A. coluzzii, duplication A. gambiae, no duplication A. gambiae, duplication
		0.05
	Figure 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-C) Id., using

variants upstream and downstream to the Ace1 duplication

(3,935 and 3,302 variants)

. Tips are color-coded according to species and duplication presence. Source alignments and complete phylogenies with supports in all nodes are available in Supplementary Material SM14 and SM15.

  . Details for each of these methods, including primer sequences, are available in Supplementary Material SM6G. The complete list of specimens, sampling times and locations, and species identification are available in Supplementary Weija and Madina samples, resistant specimens were determined as surviving a threshold concentration of pirimiphos-methyl after exposure during the larval stage, and susceptible specimens were identified as dead after a lower threshold dose (constant exposure time of 60 minutes). The exact concentration thresholds for the dose-response assay were determined on a per-site basis, and are listed in Supplementary Material SM6C. For example, A. gambiae from Aboisso were deemed resistant if they survived a 1x concentration of pirimiphos-methyl (0.25%) for 60 minutes, and susceptible if they died at

	Material SM6C and SM6E.
	All 1080 specimens were phenotyped for pirimiphos-methyl resistance using a custom dose-response
	assay with WHO standard tubes (World Health Organization 2018b). For the Aboisso (CoB te d'Ivoire), Korle-
	Bu,

0.5x concentration (0.125%). In the Obuasi (Ghana) and Baguida (Togo) populations, phenotypes were determined according to survival after different exposure times (susceptible: dead after less than 45 minutes; resistant: alive after 45 minutes or more) using a single concentration of pirimiphos-methyl (0.5x, 0.125%). The exact exposure time for each sample in intervals of 15 minutes are listed in Supplementary Material SM6C.

Out of 1080 specimens, 909 could be assigned to one of the following populations, defined as combinations of species and collection locations with determined resistance phenotypes: A. coluzzii from Aboisso, CoB te d'Ivoire (n = 55); A. coluzzii from Korle-Bu (n = 214), and Weija (n = 131), Ghana; A. gambiae from Aboisso, CoB te d'Ivoire (n = 82); A. gambiae from Baguida, Togo (n = 102); A. gambiae from Madina (n = 172), Obuasi (n = 140) and Weija (n = 13), Ghana.

Acknowledgements

We thank Sean Tomlinson (LSTM) and Chris Clarkson (Wellcome Sanger Institute) for fruitful discussions on the analyses. Funding This work was supported by the National Institute of Allergy and Infectious Diseases (R01-AI116811), the Wellcome Trust (090770/Z/09/Z; 090532/Z/09/Z; 098051), the Medical Research Council and the Department for International Development (MR/M006212/1), and the Medical Research Council (MR/P02520X/1 and MR/T001070). This UK funded award is part of the EDCTP2 programme supported by the European Union. The content of this manuscript is solely the responsibility of the authors and does not necessarily represent the official views of the National Institute of Allergy and Infectious Diseases, or the National Institutes of Health.

Availability of data and materials

Python 3.4.7 and R 3.6.2 scripts to reproduce the main analyses from the paper are available in the following Github repository: github.com/xgrau/ace1-anopheles-report.

All genome variation data has been obtained from the publicly available repositories of the Anopheles gambiae 1000 Genomes project (Phase 2 AR1) (The Anopheles gambiae 1000 Genomes Consortium 2017): https://www.malariagen.net/data/ag1000g-phase-2-ar1 (download instructions can be found in the Github repository). All quantitative PCR data are provided as Supplementary Material SM6.

We assessed the association of G280S mutations with resistance for each of the populations listed above using generalised linear models (glm function in R stats library, binomial family). The results from these tests are available in Supplementary Material SM6D.

A sub-set of these 909 specimens (n = 167; listed in Supplementary Material SM6E) was also genotyped for CNV polymorphisms in the Ace1 locus using a qPCR assay and a combination of primers for Ace1 and two non-duplicated control genes (the CYP4G16 AGAP001076; and the elongation factor AGAP005128).

Detailed methods for CNV genotyping assay are available in Supplementary Material SM6G and (Edi et al. 2014a). These samples were selected from the subset of heterozygotes in order to investigate the effect of both CNVs and fraction of 280S alleles on resistance.

We used generalised linear models (glm function in R stats library, binomial family) to assess the effect of these two variables within each of the seven populations with available CNV data (three from A. coluzzii, four from A. gambiae; listed in Table 2 and Supplementary Material SM6F). Then, we obtained the minimal significant model for each population using a stepwise reduction procedure and the BIC criterion (step function in the R stats library), and assessed their fit against a null model (intercept as the only variable) with an analysis of variance and a χ 2 test (anova function in R stats). The results from these tests are available in Supplementary Material SM6F.

Haplotype clustering and analysis of selection signals in Ace1

Haplotypes along chromosome arm 2R were reconstructed using the dataset of phased variants from the Anopheles gambiae 1000 Genomes Phase 2 data. Specifically, we retained phased variants that were biallelic, non-singletons and segregating in at least one population.

Given the paucity of phased variants within the Ace1 gene body (Supplementary SM9A, B) and the associated difficulties in phasing the G280S mutation itself, we analysed the haplotypes around nearby genetically linked variants instead. To do so, we calculated the linkage disequilibrium between position G280S and all variants located between positions 3416800 and 3659600 along 2R (i.e., 20kbp upstream and downstream of the duplication breakpoints ). Linkage disequilibrium was estimated from genotype counts in each sample using Rogers' and Huff r correlation statistic [START_REF] Rogers | Linkage disequilibrium between loci with unknown phase[END_REF] as implemented in scikit-allel v1.1.10 library (rogers_huff_r) (Miles and Harding 2017). Then, we retrieved all variants with r > 0.95 that were also present in the subset of phased variants (four in total). These four linked variants had minor allele frequencies similar to that of 280S in CoB te d'Ivoire (43%) and were located in intergenic regions -26 kbp and +12 kbp of G280S (listed in Supplementary Material SM9C).

We retrieved haplotypes surrounding each of the four tagging variants and clustered them by similarity using minimum spanning trees. Specifically, we retrieved phased variants ± 300 bp of each tagging variant (retaining variants that were biallelic, non-singletons and segregating in at least one population; these regions contained between 67 and 164 phased variants depending on the analysis; Supplementary Material SM10). We used allele presence/absence data from each haplotype to build minimum spanning trees of haplotypes based on Hamming distances (breaking edges for distances >1). This distance matrix was then used to build medium spanning trees (minimum_spanning_tree function in SciPy 1.1.0 Python library [START_REF] Jones | SciPy: Open source scientific tools for Python[END_REF], from the sparse.csgraph submodule). Tree visualizations were produced using the graphviz 2.38.0 Python library (Ellson et al.) and the hapclust library (Clarkson et al. 2018;Clarkson and Miles 2018), and clusters of haplotypes were color-coded according to the population of origin or the presence/absence of the alternate allele (Supplementary Material SM10).

We used these trees to identify groups of highly similar haplotypes linked to the 280S or wt alleles in Ace1 (Supplementary Material SM11). We calculated the profile of Garud's H statistics (moving_garud_h function, scikit-allel) [START_REF] Garud | Recent Selective Sweeps in North American Drosophila melanogaster Show Signatures of Soft Sweeps[END_REF][START_REF] Messer | Population genomics of rapid adaptation by soft selective sweeps[END_REF] and haplotype diversity (moving_haplotype_diversity in scikit-allel) around the Ace1 duplication region in windows of consecutive variants (500 variants with 20% overlap). These calculations were performed separately for the main groups of identical haplotypes (280S-linked or wt-linked) as identified around each of the tagging variants (minimum size = 100 haplotypes). In addition, we calculated the averages of these same statistics in the regions immediately upstream and downstream of the duplication breakpoints (2R:3436800 -5000 bp and 2R:3639600 + 5000 bp, respectively), and we calculated standard errors of these estimates using sample jack-knifing.

Finally, we calculated the rates of extended haplotype homozygosity decay for each cluster (Supplementary Material SM12), using 10,000 phased variants upstream and downstream of the duplication breakpoints (ehh_decay function in scikit-allel).

In parallel, we performed a phylogenomic analysis of the reconstructed haplotypes of the Ace1 CNV region (position 3436800 to 3639600) in the West African populations where duplications had been identified (75 A. coluzzii from Burkina Faso, 55 from Ghana, and 71 from CoB te d'Ivoire; and 92 A. gambiae from Burkina Faso, 12 from Ghana and 40 from Guinea; total n = 345 genomes, 690 haplotypes). Specifically, we built an alignment of phased variants along the duplication (2,787 positions) and computed a Maximum-Likelihood phylogenetic tree using IQ-TREE 1.6.10 [START_REF] Nguyen | IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies[END_REF]. We used the GTR nucleotide substitution model with ascertainment bias correction, empirical state frequencies observed from each alignment, and four gamma distribution (Γ) categories (GTR+F+ASC+G4 model in IQ-TREE). This model was selected by the IQ-TREE implementation of ModelFinder [START_REF] Kalyaanamoorthy | ModelFinder: Fast model selection for accurate phylogenetic estimates[END_REF]) (TEST mode) as the best-fitting model according to the BIC criterion. The best-fitting tree was found after 22,000 iterations (correlation coefficient threshold = 0.99). We estimated node statistical supports from 1,000 UF bootstraps [START_REF] Hoang | UFBoot2: Improving the Ultrafast Bootstrap Approximation[END_REF][START_REF] Minh | Ultrafast approximation for phylogenetic bootstrap[END_REF]. Using the same approach, we built two additional phylogenies from genomic regions located upstream and downstream of the duplication breakpoints (3,935 and 3,302 variants extracted from 50 kbp segments located -1Mb and +1Mb), finding the best-fitting trees after 7,659 and 27,400 iterations, respectively (correlation coefficient threshold = 0.99). The resulting phylogenetic trees were plotted, unrooted, using the phytools 0.6-60 [START_REF] Revell | phytools: an R package for phylogenetic comparative biology (and other things)[END_REF]) and ape 5.3 libraries (plot.phylo) [START_REF] Paradis | ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R ed. R. Schwartz[END_REF]. Phylogenies and source alignments are available as Supplementary Material SM14 and SM15.

We also used allelic frequencies within the inversion to estimate the divergence times between specimens carrying the duplication (both A. gambiae and A. coluzzii), wt A. coluzzii and wt A. gambiae. This three-way calculation of divergence times reflects the amount of allele frequency change between one of the three groups relative to the other two, and is therefore analogous to the population branch statistic [START_REF] Cavalli-Sforza | Human Diversity[END_REF][START_REF] Yi | Sequencing of 50 Human Exomes Reveals Adaptation to High Altitude[END_REF]. Following this logic, we estimated separately (i) the branch length of wt A. coluzzii relative to the separation of wt A. gambiae and the duplicated sequences, and (ii) the branch length of wt A. gambiae relative to the separation of wt A. coluzzii and the duplicated sequences. We performed these calculations in non-overlapping windows of 100 variants with the pbs function in scikitallel, and estimated standard errors using a block-jackknifing procedure.

Genetic differentiation in A. coluzzii from CoB te d'Ivoire (locate_unlinked function in scikit-allel), and repeated this process for ten iterations for each chromosome arm. This filtering procedure resulted in 791 genetically unlinked variants from both chromosomal arms. We used the alternate allele count in each variant to construct a PCA using singular value decomposition (pca function from scikit-allel), and scaling the resulting coordinates using Patterson's procedure [START_REF] Patterson | Population Structure and Eigenanalysis[END_REF].

k-mer enrichment analysis in A. coluzzii from CoB te d'Ivoire

We obtained k-mer counts for each of the 71 A. coluzzii samples using the count function in jellyfish v.

2.2.10 (Marçais and Kingsford 2011), using a k = 31 bp (parameters: -C -m 31 --out-counter-len 2 -s 300M --bf-size 10G). To reduce the computer memory footprint of the k-mer count tables, the k-mer strings were recoded as integers, split into separate lexicographical groups according to the leading nucleotides (i.e. k-mers beginning with AAA, AAC, etc. were assigned to different groups), and, within each set, sorted lexicographically again. In total, we recorded the frequency of 1,734,834,987 k-mers across 71 samples. The resulting count tables were filtered to retain k-mers showing variation in the population: we discarded kmers present in fewer than 3 samples, or absent in fewer than 3 samples. These filters removed 967,274,879 k-mers, leaving 767,560,108 k-mers for further analysis.

To test whether the frequencies of any k-mers were associated with pirimiphos-methyl resistance, normalised k-mer frequencies were obtained by dividing the k-mer counts by the total number of variant k-mers in each sample. We calculated the Spearman's rank correlation of each k-mer frequency with the resistance phenotype for each of the 767,560,108 k-mers (cor.test function in the R stats library (R Core Team 2017)), correcting for multiple testing using local estimation of false discovery rates (fdrtool function in fdrtool library [START_REF] Klaus | fdrtool: Estimation of (Local) False Discovery Rates and Higher Criticism[END_REF])), and using a significance threshold of FDR < 0.001.

Given that multiple k-mers can overlap a single mutation, it was likely that many of the k-mers identified as significant were overlapping. We therefore took the k-mers that showed a significant association with pirimiphos-methyl resistance and assembled them by joining any k-mers that overlapped perfectly over at least 10 bp. The resulting assembled k-mers were aligned against the A. gambiae reference genome using bwa mem version 0.7.12-r1034 (parameters: -T 0) [START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF]. The mapping coordinates and sequences of the significant assembled k-mers are available in Supplementary Material SM8.

A minority of assembled k-mers aligned in regions other than the Ace1 duplication. To determine whether this was due to mis-alignment, we correlated the frequency of these assembled k-mers in each sample against the Ace1 copy number. Frequency for each assembled k-mer was calculated as the mean normalised frequency of all the k-mers that were used in its assembly. Because these k-mers and the Ace1 copy number share the property of being correlated with pirimiphos-methyl resistance phenotype, they are statistically likely to have correlated frequencies even if physically unlinked. To control for this, the residuals of k-mer frequency and Ace1 copy number were calculated within each of the resistant and susceptible groups of samples. Pearson's correlation was then calculated on these residuals for each of the assembled k-mers (available in Supplementary Material SM8C).

Scripts to reproduce the k-mer analysis steps described above (counting, filtering, significance testing, assembly, mapping, and Ace1 correlation analysis) are available on Github (see Availability of data and materials).

Alignment and phylogenetic analysis of ACE proteins in animals

To obtain a candidate list of homologs of the Ace1 gene, we retrieved all genes belonging to the orthogroup 339133at33208 from OrthoDB [START_REF] Kriventseva | OrthoDB v10: sampling the diversity of animal, plant, fungal, protist, bacterial and viral genomes for evolutionary and functional annotations of orthologs[END_REF]) (682 genes in total). We aligned these candidates to a curated database of predicted proteomes from 89 complete animal genomes (listed in Supplementary Material SM2, including data sources) using Diamond 0.9.22.123 [START_REF] Buchfink | Fast and sensitive protein alignment using DIAMOND[END_REF], and retained all alignments with an identity >95% to any of the candidate queries (130 genes). Next, we performed a multi-sequence alignment with MAFFT 7.310 (1,000 rounds of iterative refinement, L-INS-i algorithm) (Katoh and Standley 2013) of these 130 genes, together with the truncated sequence of the Pacific electric ray Torpedo californica protein used for its crystal structure analysis (PDB accession number: 1W75) and the full-sequence of the ortholog from its close relative T. marmorata (Uniprot accession number: P07692). This alignment (n=132 genes) was trimmed position-wise using trimAL 1.4 (Capella-Gutiérrez et al. 2009) (automated1 mode), and 277 conserved columns were retained. A maximumlikelihood phylogenetic analysis of the trimmed alignment was then performed using IQ-TREE 1.6.10 [START_REF] Nguyen | IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies[END_REF], using a LG substitution matrix [START_REF] Le | An improved general amino acid replacement matrix[END_REF] with four Γ categories and accounting for invariant sites (LG+I+G4 model). IQ-TREE was run for 469 iterations until convergence was attained (correlation coefficient ≥ 0.99). Node statistical supports were calculated using the UF bootstrap procedure (1,000 replicates) [START_REF] Hoang | UFBoot2: Improving the Ultrafast Bootstrap Approximation[END_REF][START_REF] Minh | Ultrafast approximation for phylogenetic bootstrap[END_REF]. Alignment visualizations (Supplementary Material SM1) were obtained from Geneious 11.1.4 (Geneious 2019). Complete alignments are available as Supplementary Material SM2.

Competing interests

The authors declare no conflicts of interest. Supplementary Material SM7. Genetic differentiation in Ivorian A. coluzzii. A) Genetic differentiation statistics (Hudson's F ST and PBS) between resistant and susceptible A. coluzzii from CoB te d'Ivoire, in windows of 10,000 variants along the genome. PBS values are calculated using Angolan A. coluzzii as outgroup, and we also report a Z-score and p-value derived from a two-sided normal distribution. B) Genes overlapping the regions of high PBS (FDR < 0.001 and standardised PBS > 0). Supplementary Material SM8. k-mer analysis in Ivorian A. coluzzii. A) Alignment coordinates of k-mers that are significantly associated with pirimiphos-methyl resistance in Ivorian A. coluzzii samples from the Anopheles gambiae 1000 Genomes Phase 2. For each k-mer, we report the alignment sequence and coordinates, and whether it overlaps the Ace1 gene or the Ace1 duplication. This table includes both primary and secondary alignments as reported by bwa mem, as well as non-aligned k-mers (chromosome "NA"). B) Correlation of alignment frequencies and Ace1 copy number in each sample, for each of the 409 kmers that mapped to chromosomal arms 2R, 2L, 3R, 3L or X. We report Pearson's correlation coefficients

36

Supplementary Material legends