Genetic variation at the Cyp6m2 putative insecticide resistance locus in Anopheles gambiae and Anopheles coluzzii

Background The emergence of insecticide resistance is a major threat to malaria control programmes in Africa, with many different factors contributing to insecticide resistance in its vectors, Anopheles mosquitoes. CYP6M2 has previously been recognized as an important candidate in cytochrome P450-mediated detoxification in Anopheles. As it has been implicated in resistance against pyrethroids, organochlorines and carbamates, its broad metabolic activity makes it a potential agent in insecticide cross-resistance. Currently, allelic variation within the Cyp6m2 gene remains unknown. Methods Here, Illumina whole-genome sequence data from Phase 2 of the Anopheles gambiae 1000 Genomes Project (Ag1000G) was used to examine genetic variation in the Cyp6m2 gene across 16 populations in 13 countries comprising Anopheles gambiae and Anopheles coluzzii mosquitoes. To identify whether these alleles show evidence of selection either through potentially modified enzymatic function or by being linked to variants that change the transcriptional profile of the gene, hierarchical clustering of haplotypes, linkage disequilibrium, median joining networks and extended haplotype homozygosity analyses were performed. Results Fifteen missense biallelic substitutions at high frequency (defined as > 5% frequency in one or more populations) are found, which fall into five distinct haplotype groups that carry the main high frequency variants: A13T, D65A, E328Q, Y347F, I359V and A468S. Despite consistent reports of Cyp6m2 upregulation and metabolic activity in insecticide resistant Anophelines, no evidence of directional selection is found occurring on these variants or on the haplotype clusters in which they are found. Conclusion These results imply that emerging resistance associated with Cyp6m2 is potentially driven by distant regulatory loci such as transcriptional factors rather than by its missense variants, or that other genes are playing a more significant role in conferring metabolic resistance. Supplementary Information The online version contains supplementary material available at 10.1186/s12936-021-03757-4.

progress has stagnated and case numbers are stable or on the rise in many countries in Africa [1][2][3]. This is due to multiple factors, including the emergence of insecticide resistance, which threaten the effectiveness of vector control interventions [4].
The mechanisms of insecticide resistance are classified into four main functional categories: target-site insensitivity, penetration resistance, behavioural resistance and metabolic sequestration and detoxification; the latter being the primary focus. Mutations that diminish the binding capacity of insecticides to their targets confer target-site resistance [5], those that lead to a reduction in the insecticide absorption rate confer penetration resistance [6], those that lead to behavioural modification to circumvent the lethal effects of insecticides confer behavioural resistance [7], and those that increase insecticide detoxification and elimination rates confer metabolic resistance [8]. These mechanisms of insecticide resistance may occur concurrently within a single population or even within a single mosquito [9,10]. They have led to increasing resistance to all four common insecticide classes-pyrethroids, organochlorines, carbamates and organophosphates-in all major malaria vectors across Africa [11,12].
Metabolic resistance occurs mainly through the elevated activity of large and functionally diverse multigene enzyme families: glutathione S-transferases (GSTs), carboxylesterases (COEs) and cytochrome P450 monooxygenases (P450s) [13,14]. Metabolic resistance is considered as one of the most serious threats to mosquito control [13], because the only widely accepted occurrence of malaria vector control failure was attributed to the elevated expression of resistance-associated P450s in Anopheles funestus [15][16][17]. A comprehensive understanding of metabolic resistance is therefore vital, and it must involve disambiguating the roles that individual enzymes play and the genetic backgrounds that underlie their significance in vector populations.
The CYP6M2 enzyme exhibits complex insecticide metabolism associated with multiple binding modes for insecticides [18]. Its gene is located within a cluster of 14 Cyp6 P450 genes on chromosome 3R of Anopheles gambiae [19], and is among the 111 known P450 genes across the An. gambiae genome [20,21]. In this genomic region, Cyp6m2 is nested within a sub-cluster of P450s containing Cyp6m3 and Cyp6m4 which have also been associated with xenobiotic detoxification [22].
Cyp6m2 is notably one of the few specific P450s that have shown a consistent association with metabolic resistance [9]. Metabolic resistance is mainly assessed through transcriptional profiling of genes involved in xenobiotic detoxification. Transcriptomic experiments such as quantitative PCR and microarray assays have established a link between Cyp6m2 overexpression and the resistance phenotype in field populations of An. gambiae, Anopheles coluzzii, Anopheles arabiensis and Anopheles sinensis, irrespective of the presence of knockdown resistance (kdr) mutations such as L995F or L995S (corresponding to L1014F or L1014S in Musca domestica) in the voltage gated sodium channel (VGSC) [9,[23][24][25]. In DDT resistant An. gambiae in Ghana, Cyp6m2 has been found to be overexpressed 3.2 to 5.2-fold in combination with the upregulation of additional P450s like Cyp6z2 [22]. In DDT resistant An. coluzzii collected in Benin, Cyp6m2 was also found to be overexpressed 1.2 to 4.6-fold in combination with Gste2 from the GST gene family and in the presence of fixed kdr alleles in the Vgsc gene [26]. In Nigeria, the 2.4 to 2.7-fold upregulation of Cyp6m2 was found to be associated with high levels of permethrin resistance [9] and An. gambiae that exhibited a strong resistance to bendiocarb in In Côte d'Ivoire also had an elevated (up to eightfold) expression of the Cyp6m2 gene [24]. In the same study, transgenic expression of Cyp6m2 in Drosophila melanogaster was shown to produce resistance to both DDT and bendiocarb. In vivo functional analysis of multi-tissue overexpression induced by genetic modification has also shown Cyp6m2 to be sufficient in conferring resistance to permethrin and deltamethrin [27]. However, this overexpression also increased the susceptibility of mosquitoes to the organophosphate malathion (a negative cross-resistance mechanism). Collectively, these studies indicate that Cyp6m2 can confer metabolic resistance against insecticides in 3 of the 4 known classes: both type I and type II pyrethroids [18,22,27], organochlorines [28], and carbamates [24]. It has, therefore, a high potential for cross-resistance, which may make the problem of malaria vector control even more intractable by limiting the options available to malaria control programmes for insecticide rotation or combination. The negative cross-resistance associated with malathion hereby points to potential mitigating strategies [27].
The frequent association of Cyp6m2 with insecticide resistance described above warrants further investigation into whether there is evidence of copy number variation (CNV) or missense mutations at the locus. CNVs have been implicated in augmenting gene dosage leading to increased transcription of metabolic enzymes [29,30]. A genome-wide CNV analysis conducted on the Ag1000G dataset and described in detail elsewhere [29] found CNVs to be significantly enriched in metabolic resistance-associated gene families and to be undergoing positive selection. These CNVs were identified across P450s (such as Cyp9k1 and at both the Cyp6z3-Cyp6z1 and the Cyp6aa1-Cyp6p2 gene clusters) and GSTs (at the Gstu4-Gste3 cluster). However, CNVs across the Cyp6m2 locus were found to be rare, even in populations that are known to exhibit Cyp6m2-mediated resistance [29]. This indicates that CNVs alone are not sufficient to explain the widespread occurrence of the Cyp6m2-associated resistance phenotype: additional factors such as allelic variation might contribute to resistance associated with Cyp6m2 activity.
Allelic variation can play an additional role in P450-mediated resistance by modifying either enzyme catalytic activity or gene expression levels [31]. Allelic variation has been shown to be key in inducing high metabolic efficiency of Cyp6P9b and in conferring metabolic resistance to An. funestus [32]. Allelic variants in metabolic genes have also been identified to reliably and reproducibly associate with resistance, such as in Cyp4J5 and Coeae1d in An. gambiae, and can serve as diagnostic markers of phenotypic resistance [33]. However, there is still a paucity of information about allelic variation associated with metabolic resistance when compared to the well-characterized target-site mutations [33]. Mutations that may modulate metabolic resistance by either altering function or modifying expression in Cyp6m2 are yet to be described.
Following the consistent association of Cyp6m2 with insecticide resistance in many populations, wholegenome Illumina sequence data was examined from phase 2 of the Anopheles gambiae 1000 Genomes Project (Ag1000G) [34], which consists of 1142 wild-caught mosquitoes sequenced to a mean depth above 14× , and report a comprehensive analysis of genetic variation within the Cyp6m2 gene. The wider haplotypes around Cyp6m2 spanning across the Cyp6m sub-cluster and the larger Cyp6 supercluster are also examined for signatures of selection.

Data collection and analysis
This study ollowed the species nomenclature of Coetzee et al. [35], where An. gambiae refers to An. gambiae sensu stricto (s.s.) (S form) and An. coluzzii refers to An. gambiae s.s. (M form). A detailed description of the Ag1000G sample collection, DNA extraction, sequencing, variant calling, quality control and phasing can be found here [36]. Briefly, Anopheline samples were collected from 33 sampling sites across 16 populations in 13 countries in sub-Saharan Africa (Table 1 and Additional file 1). The sampling procedure covered different ecosystems and aimed at collecting a minimum of 30 specimens per country. The specimens consisted of An. gambiae and An. coluzzii: only An. coluzzii were sampled from Angola, both An. gambiae and An. coluzzii were sampled from Burkina Faso, while all other populations consisted of An.
gambiae, except Kenya, The Gambia and Guinea Bissau where the species identity was indeterminate.
Whole genome sequencing of all mosquitoes was performed on the Illumina HiSeq 2000 platform. The generated 100 base paired-end reads were aligned to the An. gambiae AgamP3 reference genome assembly [37] and variants were called using GATK UnifiedGenotyper. Samples with mean coverage ≤ 14× and variants with attributes that correlated with Mendelian error in genetic crosses were removed during quality control.
The SnpEff v4.1b software was used for the functional annotation of Ag1000G variant data [38] using locations from geneset AgamP4.12. All variants in transcript AGAP008212-RA with a SnpEff annotation of "missense" were regarded as nonsynonymous variants. The Cyp6m2 gene has not been shown to exhibit alternative splicing, and no alternative transcripts have been reported.

Haplotype clustering, linkage disequilibrium and mapping of haplotype clusters
To reveal the haplotype structure at Cyp6m2, Cyp6m sub-cluster, Cyp6 supercluster, HAM, ODR-2 and SH2, the Hamming distance is computed between all haplotype pairs and hierarchical clustering of haplotypes is performed. Arbitrary clustering threshold values were worked through to cut the dendrograms at genetic distances that would best highlight the most relevant clusters. Lewontin's D′ [39] was used to compute the linkage disequilibrium (LD) between all pairs of missense Cyp6m2 mutations. Image rendering for the haplotype clustering, linkage disequilibrium and haplotype cluster frequencies map was performed using the matplotlib Python package version 3.1.2 [40]. Geography handling for the haplotype cluster frequencies map was done using cartopy version 0.17.0 [41].

Haplotype networks
Haplotype networks are constructed using the medianjoining algorithm [42] implemented in Python [43]. Haplotypes carrying the main high frequency mutations were analysed with a maximum edge distance of two SNPs. The Graphviz library was used to render the networks and the composite figure was constructed in Inkscape [44].

Extended haplotype homozygosity
The core haplotype is defined on a 1689 base region spanning the Cyp6m2, from chromosome arm 3R, starting at position 6928858 and ending at position 6930547. This region was selected to ensure a 1:1 haplotype correspondence with that used in the hierarchical clustering analysis. Extended haplotype homozygosity (EHH) was computed across all core haplotypes in all populations Table 1 Allele frequencies of non-synonymous variants in Cyp6m2 coding region  [45] using scikit-allel version 1.1.9 [46]. EHH composite plots were made using the matplotlib Python package version 3.1.2 [40].

Cyp6m2 non-synonymous nucleotide variation
Short-read whole-genome sequence data from the Ag1000G phase 2 data resource [47] were used to investigate genetic variation at the Cyp6m2 locus across 16 populations of An. gambiae and An. coluzzii (n = 1142 total individuals) collected between 2000 and 2012 ( Table 1, Additional file 1). The single nucleotide polymorphisms (SNPs) studied here were discovered and QC'd using methods described elsewhere [36]. Focus is placed on SNPs that change the amino acid sequence of the CYP6M2 enzyme as they have a potential to alter gene function in Cyp6m2-associated insecticide resistance (n = 193) (Additional file 2). None of these reported mutations are synonymous or found within cis regulatory elements. As putative resistance variants under selection pressure from insecticides are expected to increase in frequency over time, allele frequencies were subsequently computed for every non-synonymous SNP in each population with reference to species and country of origin. The list is filtered to focus only on those variants that were at high frequency within populations or across populations (defined as > 5% frequency in one or more populations).
In total, this resulted in 15 non-synonymous variants that are explored further ( Table 1). Analysis of the patterns of polymorphism of Cyp6m2 from different populations showed both relative homogeneity within some geographical regions and distinct variants across different regions. The variants with the highest overall frequency were I359V (16%) and D65A (6%) ( Table 1). The most widespread variant was I359V, which was present in West, Central and East African populations of both An. gambiae and An. coluzzii. Populations with the highest frequency of I359V were Gabon (49%) and Ghana (25%) for An. gambiae, and Guinea (37.5%) for An. coluzzii. Another mutation, E328Q, was found across West Africa's An. coluzzii populations in Burkina Faso, Côte d'Ivoire, Ghana, Guinea and The Gambia and ranged in frequency from 6.2 to 13.6%. Several variants were found to exceed the 5% threshold only in one or two populations: A13T and Y347F, in Angola's An. coluzzii (39.7%) and in Kenya (52.1%) respectively and D65A only in Gabon's An. gambiae and in Kenya's populations at 42.8% and 52.1%, respectively (Table 1).

Haplotypic backgrounds of non-synonymous alleles
The Ag1000G data resource contains data that not only spans across exonic regions of any given gene, but also intronic and intergenic regions. This enables a comprehensive analysis of haplotypes that contain putative insecticide resistance alleles, but is constrained by the fact that this resource does not contain samples whose resistance status or Cyp6m2 expression levels are known.
Selection pressure acting upon missense variants or linked cis regulatory variants is likely to affect the haplotype structure of the gene. To study haplotype structure at Cyp6m2, biallelic SNPs across the entire 1689 bp Cyp6m2 gene were extracted to calculate the number of SNP differences between all pairs of 2,284 haplotypes derived from the mosquitoes. A clustering threshold of seven SNPs was manually identified, where the haplotype clusters corresponded to the haplotypes carrying the high frequency alleles (Table 1, Fig. 1). These haplotypes could mostly be grouped into five distinct clusters (labelled C1-C5): C1 contained haplotypes carrying A13T; C2 contained most haplotypes carrying D65A, A468S, and some haplotypes carrying I359V; C3 contained most haplotypes carrying both D65A and Y347F, and C5 contained haplotypes carrying E328Q. C4 contained haplotypes with no signature missense mutation (Fig. 1).
Overall, haplotype cluster distribution resembled the whole genome groupings of individuals described elsewhere using the same dataset [34]: Cluster C5 contained haplotypes from West African An. coluzzii; C4 contained An. gambiae from West, Central and near-East Africa; and the rest of the clusters contained haplotypes from samples from a single country and species (Fig. 2). The variation across the haplotypes largely showed no strict or systematic difference between the two species or across broad geographic regions, which is in line with recent whole genome sequencing reports [36].
Patterns of association among these non-synonymous variants were investigated by computing the normalized coefficient of linkage disequilibrium (D') using haplotypes from the Ag1000G phase 2 resource. Of the two highest frequency variants, I359V was found to be in perfect linkage with A468S but this was driven only by one population (Gabon) with most backgrounds carrying I359V not showing linkage with any other missense mutations ( Fig. 1 and Additional file 3: Fig. S1). D65A was in perfect linkage with A468S and Y347F, showing that D65A was almost only ever found on haplotypes carrying either A468S or Y347F. I359V and D65A, the highest frequency mutations across all populations, were found to be only in moderate linkage disequilibrium (0.36) (Additional file 3: Fig. S1). Other variants were found to be in weak linkage disequilibrium with the six main high frequency alleles and segregated independently within their own populations. While some strong associations were observed through linkage disequilibrium analysis across all populations, a deeper investigation revealed that these associations were driven by population specific dynamics in populations (such as Kenya) where bottlenecking has been an issue [36]. It is, therefore, unlikely that the identified variants are conferring some selective advantage against existing insecticide pressures.
The surrounding genomic region was explored to identify whether it showed a similar hierarchical clustering pattern to Cyp6m2, which might be indicative of either dominant demographic effects or selection acting at other linked loci that is having a major impact on variation within Cyp6m2. The downstream genes selected coded for proteins that were 1-to-1 orthologs with Drosophila melanogaster genes. Three genes were selected: ODR2 [48], HAM [49] and SH2 [50], which were 81280 bases, 457164 bases and 1198636 bases downstream of the Cyp6m2 gene, respectively. The distinctive haplotype clustering pattern observed for Cyp6m2 in the Kenya, Angola and Gabon populations persisted across these genes, indicating that in these populations, the diversity reduction in and downstream of Cyp6m2 is more likely driven by demography rather than by a selective sweep (Additional files 4, 5, 6: Fig. S2-S4). Biallelic SNPs were also extracted across the Cyp6m sub cluster of 3 genes (Cyp6m2, Cyp6m3 and Cyp6m4) and across the Cyp6 supercluster of 14 genes within which the Cyp6m sub cluster is located (Cyp6s2, Cyp6s1, Cyp6r1, Cyp6n2, Cyp6y2, Cyp6y1, Cyp6m1, Cyp6n1, Cyp6m2, Cyp6m3, Cyp6m4, Cyp6z3, Cyp6z2 and Cyp6z1), and performed hierarchical clustering across these regions as described The genetic backgrounds carrying these alleles were examined further by constructing median joining networks (MJNs) [42] using the Ag1000G Phase 2 haplotype data. This enabled the resolving of the radiation of DNA substitutions arising on haplotypes carrying the identified variants. It also allowed the reconstruction and positioning of intermediate haplotypes while revealing the non-hierarchical relationships between haplotypes that could not be resolved by hierarchical clustering alone. The MJNs were constructed with reference to a maximum edge distance of two SNPs. This ensured that the connected components captured only closely related haplotypes. The resulting MJNs had a close correspondence with the hierarchical clustering output in assignment of haplotypes to clusters (88% overall concordance across all clusters).
The median joining networks showed more clearly the distinctive demographic stratification of the high frequency variants that was highlighted by the hierarchical clustering networks (Fig. 3). Most nodes containing secondary variants arising from the main nodes were small, which is inconsistent with directional selection where larger nodes are expected. Only one of the I359V nodes contained haplotypes from mosquitoes of both species, however the secondary nodes did not contain haplotypes from more than one species. This indicates that although I359V is shared by both An. gambiae and An. coluzzii, it is unlikely that this is because of an introgression event across the Cyp6m2 gene.

Positive selection of non-synonymous alleles
Extended Haplotype Heterozygosity (EHH) decay [45] was calculated to explore evidence for directional selection on the haplotypes carrying high frequency nonsynonymous variants. It is expected that the presence of ongoing or recent directional selection pressure would lead to the increase in frequency of haplotypes, which on average will have longer regions of haplotype homozygosity relative to haplotypes that are not under selection. This diversity reduction would produce signatures of selection that would be conspicuous across a large genomic region, including across nearby cis-regulatory elements. EHH analysis would therefore be able to detect diversity reduction caused by ongoing directional selection being driven either by amino acid substitutions identified within the gene or by mutations within cis-acting elements next to the gene that may be under selection.
To perform the EHH decay analysis, a core region of 1689 bases that spans across the entire gene was defined. This was identical to what was used to differentiate the Haplotypes in group C5 carry the E328Q allele. Haplotypes in group C4 had no defining non-synonymous variant, and wild type (wt) haplotypes were all those that did not fall within the C1-C5 clusters identified haplotype groups though hierarchical clustering. This region contained multiple distinct haplotypes above 1% frequency within the cohort, including haplotypes corresponding to the C1-C5 haplotype clusters. All haplotypes that did not correspond to C1-C5 were considered to be wild type (wt). Although there were several different haplotypes in each population that fit this description, they are not distinguished in this paper and all these are called wild type, as Cyp6m2 has no known resistance alleles and a true wild type remains to be discovered. EHH decay was then computed for each core haplotype up to 200 kilobases upstream and downstream (Additional file 9: Fig. S7): beyond 200 kb, the EHH had decayed to zero.
It is noted that haplotype clusters containing high frequency variants (C1-C5) did not exhibit a significantly slower EHH decay relative to the wild types, showing no evidence of positive selection. However, one Kenyan wild type haplotype group had a dramatically slower EHH decay relative to wild type haplotypes from other populations. In order to account for this difference within wild type groups across multiple populations and to reveal potential signs of selection that would be obscured by a collective analysis across all populations, the haplotypes were separated by population and species and EHH decay was recomputed for each core haplotype as above.
Kenyan mosquito populations are known to have an extreme demographic history, as they have experienced a severe recent bottleneck, and the Angola and Gabon populations are known to be geographically unique populations which are strongly differentiated from all other populations [36]. Hence, their haplotypes exhibited a considerably slower decay than West African haplotypes (first three panels: Fig. 4). However, the putative resistance haplotypes C1-C5 did not experience a slower EHH decay relative to their wild type haplotypes, showing no evidence of positive selection acting upon those haplotypes in those populations. As expected, the West African An. coluzzii haplotypes exhibited a much faster decay of EHH than specimens from Kenya, Angola, or Gabon, highlighting the demographic differences previously observed for these collections [36] (last three panels, Fig. 4). The C5 haplotype was a promising candidate for potential selection as it occurred within a more diverse population, and it was interesting to note that some wild type haplotypes in Côte d'Ivoire's An. coluzzii had a slightly slower decay than others within West Africa (fifth panel, Fig. 4). However, these haplotypes were not part of the C5 cluster, and did not carry the widespread E328Q mutation. The C5 haplotype did not exhibit a dramatically slower decay of EHH than wild type haplotypes in the populations in which it was found, suggesting that it is not under positive selection.

Discussion
Cyp6m2 has been implicated in many Anopheles populations as a key P450 that contributes to the insecticide resistance phenotype [9,18,24,28]. It has been reported that allelic variants across some P450s can affect enzyme conformational dynamics and substrate binding affinity [32], offering potential mechanisms that may modulate enzyme activity and efficiency, and thus account for additional Cyp6m2 resistance where CNVs alone may not suffice. However, little is also known about Cyp6m2 allelic variation across Africa.
In this study, a comprehensive account of the distribution of amino acid substitutions occurring within the Cyp6m2 gene is reported. The haplotype structure of the gene is also examined to probe for selective sweeps by performing hierarchical clustering of haplotypes. The genetic background upon which the missense variants are found is also examined by plotting both median joining networks and decay of extended haplotype homozygosity, which are useful for revealing signatures of selection. It is noted that the distinct haplotype groups therein are stratified demographically and largely correspond to signature missense variants found in specific populations. This is in contrast to the strong signals of recent positive selection at other cytochrome P450 gene loci such as at Cyp6p3 [36] which is often upregulated in tandem with Cyp6m2 in multiple pyrethroid resistant populations [9,51,52].
It is still unclear how the identified non-synonymous variants may modulate Cyp6m2 binding activity, in either the presence or absence of multiple competitive substrates and metabolites. The two aromatic residues (Phe 108 and Phe 121) that have been previously identified to be vital in deltamethrin orientation in the Cyp6m2 active site [18] were not found to contain high frequency variants in our dataset. Candidate-gene based genome-wide association studies can be used to supplement this exploration of natural population variation in the identification of variants that reproducibly associate with resistance [33]. Computational techniques can also provide a quick and cheaper alternative to experimental approaches to explore and predict the effects of these identified variants on enzyme-substrate interactions [53]. However, subsequent in vivo functional analysis will be necessary to validate any causative links between variants and resistance phenotypes. This may involve transgenic expression of mutant enzymes, either in cell lines or in mosquito None of the haplotype groups identified that carried missense variants were found to be under directional selection. This is despite the existence of a widespread variant (E328Q) in both An. gambiae and An. coluzzii linked to a geographic region (West Africa) where Cyp6m2 upregulation has been associated with emerging metabolic resistance [24,51]. In An. coluzzii originating from both Côte d'Ivoire and Ghana, the C5 haplotype that carried E328Q was shown to have an even faster decay of EHH than the wild type haplotypes, further indicating an absence of directional selection. The stratification of other main haplotype clusters from Angola (C1), Gabon (C2) and Kenya (C3) was also consistent with the strong demographic differentiation and overall reduced heterozygosity of these populations described elsewhere [36].
While the genomic data quality across the Cyp6m2 gene and its putative promoter region was satisfactory, there was a ~ 10,000 base region of inaccessibility upstream of Cyp6m2 that cut across the intergenic region into Cyp6n1 [54]. A similar inaccessible region was also present 1 kb downstream of the gene in the intergenic region between Cyp6m2 and Cyp6m3, which is likely caused by the presence of repeats that inhibit read mapping. Although it is possible that the upstream region of inaccessibility could contain a regulatory variant that is susceptible to selection, it is unlikely to obscure signatures of selection.
It has been shown in multiple studies that target-site resistance (i.e. VGSC-kdr) provides a strong persistent baseline of resistance as it rises towards fixation within populations [55]. In the presence of insecticide selection pressure, target-site mutations and metabolic resistance have also been shown to act synergistically to confer a stronger resistance phenotype to pyrethroids [33,56]. While signatures of selection have previously been identified in some metabolic gene clusters within populations that have a high kdr frequency [36], further studies need to examine whether directional selection occurring on one locus can obscure selection on another locus. To resolve this conundrum, genomic analysis must be performed on populations sampled across generations and whose transcriptomic and phenotypic characteristics are known, in order to tease out the individual contributions of specific sources of resistance.
Independent studies employing different experimental designs have also shown that metabolic resistance manifests as a cascade of multiple upregulated genes [57]. These genes, like Cyp6m2, are part of the normal cellular mechanism for xenobiotic detoxification that involves a linked, coordinated response of large multi-gene enzyme families in complicated pathways. Therefore, it is likely that identifying signatures of selection due to insecticide pressure will involve thorough analysis across this vast network. The Cap 'n' Collar isoform-C (CncC) transcription factor sub-family has been shown to work in tandem with other transcription factors to regulate the transcription of phase I, II and III detoxification loci of multiple insects such as Culex quinquefasciatus and D. melanogaster [58,59]. CncC knockdown or upregulation has been shown to directly affect phenotypic resistance in An. gambiae as well, modulating the expression of key P450s enzymes such as Cyp6z2, Cyp6z3 and Cyp6m2 that are located in the same genomic region [58,59]. Given that no evidence of selection on amino acid variants were detected in the Cyp6m2 gene, it is possible that the emergence of Cyp6m2 associated resistance is being driven by selection pressures acting upon genes coding for distant regulatory proteins, such as transcription factors. These transcription factors can regulate downstream gene expression across large genomic distances. These transcription factors have also been implicated in the differential expression of other detoxification enzyme families also associated with insecticide resistance (GSTs, COEs, UDP-glucuronosyltransferases (UGTs) and ABC transporters). It is, therefore, likely that the centre of selection leading to the Cyp6m2 associated resistance phenotype will be identified through whole genome selection scans of susceptible and resistant populations rather than by single loci analysis. Further research on Anopheline epigenomics, transcriptomics, proteomics and systems biology will also be game changers in mapping the complex regulatory network of insecticide resistance, aiding the identification of critical targets and the development of new strategies to control the spread of metabolic insecticide resistance.

Conclusion
The scale-up of insecticide-based interventions has caused increased selection pressure and higher levels of insecticide resistance across Africa. While the CYP6M2 enzyme has been associated with emerging metabolic resistance in Africa, our data indicates that allelic variation within the Cyp6m2 gene itself or across its Cyp6 supercluster has not been subject to recent positive selection in any of the populations sampled. This is in contrast to other Cytochrome P450 genes where CNV alleles are clearly under strong selection. Our results do not rule out a role for Cyp6m2 in insecticide resistance in natural populations, but highlight the need for a deeper understanding of the regulatory networks affecting Cytochrome P450 gene expression in malaria vectors. This will require large-scale, holistic experimental work