Genetic Diversity of Cytochrome P450s CYP6M2 and CYP6P4 Associated with Pyrethroid Resistance in the Major Malaria Vectors Anopheles coluzzii and Anopheles gambiae from Yaoundé, Cameroon

Assessing the genetic diversity of metabolic resistance genes, such as cytochrome P450s, helps to understand the dynamics and evolution of resistance in the field. Here, we analyzed the polymorphisms of CYP6M2 and CYP6P4, associated with pyrethroid resistance in Anopheles coluzzii and Anopheles gambiae, to detect potential resistance markers. Field-caught resistant mosquitos and susceptible lab strains were crossed, and F4 was exposed to permethrin for 15 min and 90 min to discriminate highly susceptible (HS) and highly resistant (HR) mosquitos, respectively. Significant permethrin mortality reduction was observed after pre-exposure to PBO, suggesting the gene involvement of P450s. qPCR analysis revealed significant overexpression of CYP6M2 (FC = 19.57 [95% CI 13.96–25.18] for An. coluzzii; 10.16 [7.86–12.46] for An. gambiae) and CYP6P4 (FC = 6.73 [6.15–7.30] An. coluzzii; 23.62 [26.48–20.76] An. gambiae). Full-gene and ≈1 kb upstream were sequenced. For CYP6M2, the upstream region shows low diversity in HR and HS (overall Hd = 0.49, π = 0.018), whereas the full-gene shows allelic-variation but without evidence of ongoing selection. CYP6P4 upstream region showed a lower diversity in HR (Hd = 0.48) than HS (Hd = 0.86) of An. gambiae. These results highlighted that CYP6P4-associated resistance is potentially driven by modification in upstream region. However, further work is needed to determine the real causative variants that will help design rapid detection tools.


Introduction
Malaria cases significantly decreased from 2000 to 2015, mainly due to the scaling up of vector control interventions. This decrease was attributed largely to the enormous distribution of Insecticidal Treated Nets (ITNs) with the number of households in sub-Saharan Africa having at least one net increasing from 5% in 2000 to 65% in 2020 [1,2]. Unfortunately, the indiscriminate use of pesticides in the agricultural sector combined with indoor residual spraying (IRS) and the massive distribution of pyrethroid-impregnated nets selected resistance in vector populations in endemic areas [3][4][5][6]. Moreover, the COVID-19 pandemic disrupted vector control activities in most countries, leading to a rebound in mortality from 409,000 malaria associated deaths in 2019 to 627,000 in 2020 [2]. Insecticide resistance, which now extends to almost all classes of insecticides used in public health, remains a major concern for the effectiveness of vector control tools [7]. It is therefore imperative to implement appropriate insecticide resistance management strategies as a proactive approach to maintain the effectiveness of these interventions and to reduce the negative impact of resistance.
The development of insecticide resistance is a complex process depending directly on genetic, physiological, behavioral and ecological factors [8]. In Anopheles mosquitos, resistance is mainly conferred by target site and metabolic mechanisms [6,9]. Metabolic resistance involves an increase in production, or more efficient forms, of detoxication enzymes capable of breaking down the insecticide before it reaches its target; the cytochrome P450 family is commonly considered the most important in insecticide detoxification, especially for the metabolism of pyrethroid [9,10]. Cytochrome P450-mediated resistance can often be reversed by the use of synergists, such as piperonyl butoxide (PBO), which inhibit the activity of this enzyme family. Large-scale field trials by Protopopoff et al. [11] and Staedke et al. [12] have shown that nets incorporating this synergist are effective against P450-mediated resistance. Given the continued importance of pyrethroids in ITNs, it is imperative to understand the genetic mechanisms that impact their efficacy to optimize the deployment of PBO-nets and to establish good insecticide resistance management protocols.
Members of the CYP family, including CYP6M2 and CYP6P4, have been functionally associated with metabolic resistance to at least one class of insecticides widely used in vector control, and both are pyrethroid metabolizers [13,14]. CYP6M2 is overexpressed in multiple populations resistant to type I and II pyrethroids [13,15], organochlorines [16] and carbamates [17]. Interrogation of the Anopheles gambiae 100 Genomes project (Ag1000G) identified several amino acid substitutions within CYP6M2, grouping into 5 haplotype groups, but none of these was directly linked with resistance [18]. The CYP6P4 gene has been established as a resistance-associated gene that explains the resistance to permethrin observed in the field in An. arabiensis populations in the Sudano-Sahelian region [14,19]. Recently, Njoroge et al. [20] identified in An. gambiae a selective sweep localized within a cluster of P450 genes including CYP6P4 in mosquitos from Uganda. A haplotype, containing three mutations, I236M in CYP6P4, an upstream transposable element (TE) insertion and a CYP6AA1 duplication is strongly associated with pyrethroid resistance in An. gambiae from Uganda. In particular, the I236M mutation is strongly associated with increased overexpression of the CYP6P4 gene capable of metabolizing pyrethroids, notably deltamethrin. Moreover, the overexpression of a group of genes, including CYP6P4 in Anopheles populations, was associated with the loss of efficacy of some ITNs in Tanzania [21].
In Cameroon, P450s have been found widely implicated in metabolic resistance to pyrethroids, DDT and Bendiocarb in An. gambiae s.l. populations [22][23][24] and are at least partially responsible for the reduced efficacy of bed nets for An. gambiae s.s. [25]. The challenge in detecting simple point-mutation markers for metabolic resistance in contrast to target-site resistance, such as knockdown resistance (kdr), includes the size of the gene families involved in the detoxication process, the redundancy between their members and the multiple mechanisms by which metabolic resistance can arise [26]. A major milestone was recently reached with the development of field-applicable diagnostic tools for P450-mediated resistance in An. funestus, taking advantage of a mutation linked to the overexpression of CYP6P9a and the variants in cis-regulatory elements of CYP6P9b [27][28][29]. Such tools are less readily available for An. gambiae although, as described above, a triple mutant haplotype in CYP6P4 was recently linked to resistance in East Africa [20].
Transcriptomic analysis in field populations of An. gambiae, An. coluzzii and An. arabiensis have revealed several overexpressed detoxication genes, irrespective of the presence of kdrmutations [4,17,22,30]. Similarly, studies conducted by the Ag1000G have identified strong selective sweeps in clusters of CYPs known to be implicated in resistance [31]. The elucidation of the roles played by major resistance genes and the changes in the genetic diversity underlying their expression are essential to design DNA-based diagnostic tools to detect and track metabolic resistance in the field.
In this study, we analyzed the polymorphism of CYP6M2 and CYP6P4 in An. gambiae and An. coluzzii to detect potential resistance markers following the same approach as Weedall et al. [29] based on target DNA-sequencing of highly permethrin-resistant and highly susceptible samples selected from reciprocal crosses. We report a comprehensive analysis of the putative promoter and full-gene region of these genes in An. coluzzii and An. gambiae from Yaoundé-Cameroon.  51.67 E), two neighborhoods of Yaoundé, the capital city of Cameroon. These collection areas have been chosen based on species distribution and the detection of pyrethroid resistance as described in previous studies [22,32,33] predicting the presence of An. coluzzii in Ngousso and An. gambiae in Nkolondom, with detected metabolic resistance. Larvae were sampled by dipping methods [34], in more than 20 breeding sites and pooled areas, then taken to the insectarium of the Centre for Research in Infectious Diseases (CRID) in Yaoundé, where they were reared under standard laboratory conditions to the adult stage.

Insecticide Bioassays and PBO-Based Synergist Test
Bioassays were conducted with female mosquitos obtained from field-collected larvae. Two-to five-day-old, unfed F 0 females of An. coluzzii and An. gambiae were exposed for 1 h in WHO tube tests according to WHO procedures to 0.75% permethrin (1×). After observing resistance (mortality less than 95% [35]), we followed up with higher concentrations (3.5% (5×), then 7.5% (10×)) for the evaluation of resistance intensity [35].
To determine the possible implication of cytochromes P450 in the observed phenotypic resistance to pyrethroids, a test with synergist (PBO) that inhibits the effect of those P450s was done. F 0 females were exposed to 4% PBO [36] for 1 h followed by an immediate second exposure to permethrin (0.75%) for another 1 h. The mortality was determined 24 h post-exposure and compared with mortality obtained for the mosquitos exposed to the 0.75% permethrin only using a Chi-square test of significance.

Mosquito Selection
In the high pyrethroid resistance context of our collection areas, it was difficult to have enough susceptible samples from the field to run a comparative analysis between resistant and susceptible. Therefore, crossings between selected field permethrin-resistant and lab-susceptible strains from the same species were done as follows: An. coluzzii with Ngousso lab strain; An. gambiae with Kisumu lab strain. This follows the process described by Wondji et al. [37] and Cattel et al. [38]. To have distinguished phenotypes, F4 individual progenies obtained from those crosses were then segregated based on their resistance phenotype by exposing 3-to 5-day-old females to 0.75% permethrin at two different exposure times: 15 min and 90 min. After 24 h, dead individuals from the 15 min exposed batch were considered highly susceptible (HS) groups, while those who survived after the 90 min exposure time were considered highly resistant (HR) groups [37]. Samples were then stored on silica gel for the dead and in RNAlater ® for the survivors for subsequent molecular analysis.

Molecular Analysis
Genomic DNA (gDNA) was extracted using the LIVAK method [39] from 30 female mosquitos (F 0 ) morphologically identified as belonging to the An. gambiae complex [40,41] Genes 2023, 14, 52 4 of 19 originating from both collection sites and two sets each of 60 HR and HS mosquitos (F 4 ) for both An. coluzzii and An. gambiae. Molecular identification was achieved through a polymerase chain reaction (PCR) described by Santolamazza et al. [42].

Genotyping of Kdr-Mutation Alleles
To establish the frequency of each of the kdr-mutations (L1014F, L1014S, N175Y), genomic DNA of previously extracted samples were used to run the TaqMan assays described by Bass et al. [43] and Jones et al. [44]. The reaction mixture of 10 µL final volume containing 1 × Sensimix (Bioline), 1 × primer/probe mix and 1 µL template DNA was used for this assay. The probes were labeled with two distinct fluorophores [43]: FAM to detect the resistant alleles (L1014F-kdr and L1014S-kdr) and HEX to detect the susceptible allele (L1014). For N1575Y, the detection was performed by using the primer/probe described by Jones et al. [44]. The assay was performed following cycling conditions of 95 • C for 10 min, 40 cycles at 95 • C for 15 s and 60 • C for 1 min.

Investigation of the Expression Profile of Six Main Candidate-Resistance Genes
The expression profiles of six genes were investigated by quantitative Reverse Transcriptase PCR (qPCR): CYP6M2, CYP6P3, CYP9K1, CYP6Z1, CYP6Z2 and CYP6P4. For three biological replicates of 10 F 0 females that recovered 24 h post-exposure to permethrin and 3 batches of 10 susceptible females of laboratory Ngousso (An. coluzzii) or Kisumu (An. gambiae), total RNA was extracted using the Arcturus PicoPure RNA isolation kit (Life Technologies, Carlsbad, CA, USA), followed by cDNA synthesis using Superscript III (Invitrogen ® ) with oligo-dT 20 and RNase H, according to the manufacturer instructions. A qPCR reaction was performed following the protocol of Riveron et al. [45]. Fold changes and expression levels of each gene in field-caught permethrin-resistant (R) mosquitos and susceptible laboratory strain (S) samples were computed for both An. coluzzii and An. gambiae according to the 2-∆∆CT method [46] following standardization with the housekeeping genes ribosomal protein S7 (RPS7; VectorBase assession number: AGAP010592) and Elongation factor (VectorBase assession number: AGAP005128). The primers used here are those previously described by Mavridis et al. [47].

Amplification and Direct Sequencing of Upstream and Full Gene Regions of CYP6M2 and CYP6P4
Exploiting the An. gambiae reference genome (VectorBase accession number: AGAP008212 and AGAP002867, respectively), both An. coluzzii and An. gambiae primers were designed for the amplification of a 1-kb upstream and full-gene length of CYP6M2 and CYP6P4 using Primer3 online software (v4.0.0; http://bioinfo.ut.ee/primer3/, accessed on 23 August 2019). The upstream region and full-gene length of each species were amplified from the gDNA of 10 samples of susceptible strain and two sets each of 10 HR and 10 HS female  Table S1). Briefly, the conditions are: initial denaturation at 95 • C for 5 min followed by 35 cycles each of 30 s at 94 • C (denaturation), 2 min at primer annealing temperature, 1 min 30 s at 72 • C (elongation) and a final extension step at 72 • C for 10 min. All the PCR products CYP6M2 and CYP6P4 were purified using QIAquick ® PCR Purification Kit (Qiagen, Hilden, Germany) and ExoSAP-ITTM, respectively, according to the manufacturer's protocol, and sent for sequencing at GENEWIZ UK LTD ® using the same PCR primers.

In Silico Prediction of Promoter Region of CYP6M2 and CYP6P4
The CYP6M2 and CYP6P4 1-kb upstream region was analyzed to identify any regulatory features, such as TATA box, CCAAT box and GC box, found in the sequence using GPMiner [48] as done to detect key resistance variants for CYP6P9a/b in An. funestus [28,29]. Additionally, the potential transcription factor binding sites (TFBSs) involved in the regulatory mechanism of resistance were identified in each upstream sequence [49,50], using a prediction software ALGGEN Promo (http://alggen.lsi.upc.es/cgi-bin/promo_v3/, accessed on 20 June 2022) in which positional weight matrices (PWM) are constructed from eukaryotic binding sites extracted from TRANSFAC [51], with a focus on human and mouse.

Polymorphism-Based Analysis of Genomic DNA Sequence
Polymorphism analysis was carried out through the manual examination of the sequence traces using BioEdit version 7.2.3.0 [52] and/or nucleotides/amino acid differences from multiple sequence alignments. Genetic parameters, such as the number of haplotypes (h) and their diversity (Hd); number of polymorphic sites (S); and nucleotide diversity (π) were computed using DnaSP v6 [53]. Population Analysis with Reticulate Trees (PopART) version 1.7 software was used to construct haplotype networks showing the distribution of haplotypes per study site or species. MEGA 6.06 software was used to construct the Maximum Likelihood tree based on a specific parameter distance model according to the specific region with 1000 bootstrap replicates [54].
After polymorphism analysis, simple assays of Restriction fragment length polymorphism (RFLP) or Allele-specific PCR (AS-PCR) based on key mutations consistently found in higher frequencies in HR than in HS from F4 populations and the susceptible lab strain were designed to assess their frequency in field-caught populations.

Detection of Key Polymorphisms in An. gambiae CYP6M2 and CYP6P4
PrimerQuest tools (https://eu.idtdna.com/Primerquest/ accessed on 20 June 2022) were used to design a PCR-RFLP and AS-PCR that could detect the key mutations found in CYP6M2 and CYP6P4 genes. These mutations were selected based on their frequency of occurrence by comparing HR, HS and susceptible lab strain Kisumu. AS-PCR assays were then designed for the genotyping of non-synonymous mutations A392S in CYP6M2 and C168S in CYP6P4, using different sets of primers (Supplementary Table S2). Primers were designed manually to match the mutation, and an additional mismatched nucleotide was added in the 3rd nucleotide from the 3 end of each inner primer to enhance the specificity as done previously by Tchouakui et al. [55]. RFLP PCRs, capable of discrimination between the mutant allele of the upstream region of CYP6M2 and for both upstream and gene regions of CYP6P4, were designed based on specific restriction enzymes: BsrDI for the upstream region of CYP6M2, PvuII and EagI for CYP6P4, respectively. The restriction site was identified in the sequences of HR samples and absent in the sequences of Kisumu. Ten microliters of the digestion mix made of 1 µL of 10× NEBuffer 2.1, 0.2 µL of 2 U of respective restriction enzymes (New England Biolabs, Ipswich, MA, USA), 5 µL of PCR product and 3.8 µL of dH 2 O was incubated at specific temperatures for 2 h, as presented in Supplementary Table S2. According to the specific region, PCR products obtained in all the assays were separated on 2% agarose gel, stained with Midori Green Advance DNA Stain (Nippon Genetics Europe GmbH) and visualized using a gel imaging system to confirm the product sizes (Supplementary Table S2).
To evaluate the assay described above, the F4 progeny from a cross between fieldcaught permethrin-resistant from Nkolondom and susceptible laboratory strain was genotyped for 30-50 mosquitos and correlated with the established resistance phenotype using the odds ratio and Fisher's exact test [56].

Data Analysis
The results of bioassays were interpreted based on percentage mortalities with standard error on the means (SEM) calculated and corrected when needed following WHO protocol [35]. Results of mortalities from synergist-permethrin exposure were compared with values obtained from exposure to pyrethroid alone using a two-tailed Chi-Square test of independence, with a level of significance set as p < 0.05, then implemented in GraphPad Prism 7.02 (GraphPad Inc., La Jolla, CA, USA). To investigate the association between specific mutations and the ability of the mosquitos to survive, Vassarstats was used to estimate the Odds Ratio (OR) based on a fisher exact probability test with a 2 × 2 contingency table.

Species Identification and Susceptibility to Insecticides
The molecular identification of collected mosquitos confirmed the spatial distribution of the two species of An. gambiae complex: on 90 randomly tested samples per area, An. coluzzii was found exclusively in Ngousso, the urban area, whereas An. gambiae was found exclusively in Nkolondom, the peri-urban area.
F 0 populations of An. coluzzii showed resistance to permethrin with a mortality rate of 8.33% ± 1.43% for all the samples tested. A pre-exposure to PBO revealed a 23.67% recovery of susceptibility ( Figure 1a), with mortality rising to 32% ± 2.74% (χ 2 = 64.37, p < 0.0001). The susceptibility of F 0 An. gambiae to permethrin, showed a similar trend to what was observed in An. coluzzii with a mortality of 9.98% ± 3.43%. An increase in mortality with PBO was also observed, up to 25% ± 9.65% (χ 2 = 56.39, p < 0.0001), representing a recovery of 16% ( Figure 1a). These results suggest an involvement of cytochromes P450 in the resistance pattern observed in this population.
tocol [35]. Results of mortalities from synergist-permethrin exposure were compared with values obtained from exposure to pyrethroid alone using a two-tailed Chi-Square test of independence, with a level of significance set as p < 0.05, then implemented in GraphPad Prism 7.02 (GraphPad Inc., La Jolla, CA, USA). To investigate the association between specific mutations and the ability of the mosquitos to survive, Vassarstats was used to estimate the Odds Ratio (OR) based on a fisher exact probability test with a 2 × 2 contingency table.

Species Identification and Susceptibility to Insecticides
The molecular identification of collected mosquitos confirmed the spatial distribution of the two species of An. gambiae complex: on 90 randomly tested samples per area, An. coluzzii was found exclusively in Ngousso, the urban area, whereas An. gambiae was found exclusively in Nkolondom, the peri-urban area.
F0 populations of An. coluzzii showed resistance to permethrin with a mortality rate of 8.33% ± 1.43% for all the samples tested. A pre-exposure to PBO revealed a 23.67% recovery of susceptibility ( Figure 1a), with mortality rising to 32% ± 2.74% (χ 2 = 64.37, p < 0.0001). The susceptibility of F0 An. gambiae to permethrin, showed a similar trend to what was observed in An. coluzzii with a mortality of 9.98% ± 3.43%. An increase in mortality with PBO was also observed, up to 25% ± 9.65% (χ 2 = 56.39, p < 0.0001), representing a recovery of 16% ( Figure 1a). These results suggest an involvement of cytochromes P450 in the resistance pattern observed in this population.

Expression Profile of CYP6M2 and CYP6P4
Resistant Anopheles populations from Nkolondom and Ngousso showed significant variation in metabolic gene expression between the two species for some of the tested genes. Both CYP6P4 and CYP6M2 were overexpressed in An. coluzzii F 0 at a significant level (p < 0.05) with a fold change (FC) of 19 (Figure 1b).

Detection of the Knockdown Resistance Mutations
Among the three kdr mutations tested, only the 1014F kdr-mutations were detected in any of the sites/species. This mutation was found fixed in the An. gambiae population from Nkolondom (1.00) and nearly fixed (0.93) in the An. coluzzii population from Ngousso (Supplementary Table S3). In F4 from crosses for both species, the HR samples showed the highest 1014F allelic frequency, 0.80 for An. coluzzii and 0.84 for An. gambiae (Figure 2), while in the HS sample, we obtained 0.07 for An. coluzzii to 0.10 for An. gambiae on 30 samples per batch.

Expression Profile of CYP6M2 and CYP6P4
Resistant Anopheles populations from Nkolondom and Ngousso showed significant variation in metabolic gene expression between the two species for some of the tested genes. Both CYP6P4 and CYP6M2 were overexpressed in An. coluzzii F0 at a significant level (p < 0.05) with a fold change (FC) of 19.57 [95% CI 13.96-25.18] for CYP6M2 and 6.73 [6.15-7.30] for CYP6P4 compared to the susceptible laboratory colony (Figure 1b). In An. gambiae, those genes were also found to be overexpressed (CYP6M2: FC = 10. 16 (Figure 1b).

Detection of the Knockdown Resistance Mutations
Among the three kdr mutations tested, only the 1014F kdr-mutations were detected in any of the sites/species. This mutation was found fixed in the An. gambiae population from Nkolondom (1.00) and nearly fixed (0.93) in the An. coluzzii population from Ngousso (Supplementary Table S3). In F4 from crosses for both species, the HR samples showed the highest 1014F allelic frequency, 0.80 for An. coluzzii and 0.84 for An. gambiae (Figure 2), while in the HS sample, we obtained 0.07 for An. coluzzii to 0.10 for An. gambiae on 30 samples per batch. Correlation between L1014F kdr and permethrin-resistance phenotype in An. coluzzii and An. gambiae from F4. (a1 and b1) Distribution of kdr L1014F genotypes in dead and alive F4 mosquitos from An. coluzzii and An. gambiae crossings, respectively, after WHO bioassays with 0.75% permethrin. (a2, a3 and b2, b3) Association between frequency of 1014F and ability of F4 mosquitos, respectively, from An. coluzzii and An. gambiae crossings to survive to WHO bioassays with 0.75% permethrin, showing a very strong link between 1014F kdr mutation and permethrin-resistance phenotype. HR, Highly resistant; HS, Highly susceptible; FF, homozygous resistant; FL, heterozygous resistant; LL, homozygous susceptible; F, resistant allele (L1014F); L, susceptible allele (L1014). Significance is shown by *p < 0.05, ****p < 0.0001, as estimated using Fisher's exact test. Correlation between L1014F kdr and permethrin-resistance phenotype in An. coluzzii and An. gambiae from F4. (a1,b1) Distribution of kdr L1014F genotypes in dead and alive F4 mosquitos from An. coluzzii and An. gambiae crossings, respectively, after WHO bioassays with 0.75% permethrin. (a2,a3,b2,b3) Association between frequency of 1014F and ability of F4 mosquitos, respectively, from An. coluzzii and An. gambiae crossings to survive to WHO bioassays with 0.75% permethrin, showing a very strong link between 1014F kdr mutation and permethrin-resistance phenotype. HR, Highly resistant; HS, Highly susceptible; FF, homozygous resistant; FL, heterozygous resistant; LL, homozygous susceptible; F, resistant allele (L1014F); L, susceptible allele (L1014). Significance is shown by * p < 0.05, **** p < 0.0001, as estimated using Fisher's exact test.

Polymorphism Analyses of the Upstream Region of CYP6M2
To have a view of the potential regulatory elements driving the overexpression of CYP6M2, An. gambiae mosquito samples were used as a model. After cleaning and aligning the sequences, an indel of 8 bp (TAGTTACT) was found and is linked to the presence of the G/T mutation (Guanine replaced Thymine at position 316) in the HR (7/8 sequences) and HS samples (5/6). This Indel was absent for all the Kisumu tested as the G/T mutation is here absent. Supplementary Table S4 shows the frequency of this indel across the groups used. The core promoter elements detected by GPMiner include TATA boxes (9 in the insertion carriers and 6 in the others), CCAAT box (1 in the insertion carriers and 2 in deletion carriers) and 1 GC box (Supplementary Figure S1). Furthermore, the CYP6M2 upstream region also contains transcription factors binding sites in all the samples tested, including 2 AhR/ARNT, 4 nrf2/MAF as well as several minority sites for GATA and MYB transcription factors (Supplementary Figure S1).
The 932 bp upstream sequences of CYP6M2 were aligned for 20 individuals of An. coluzzii (5 HR, 8 HS and 7 Ngousso strain) and 17 of An. gambiae (8 HR, 6 HS and 6 Kisumu lab strain) (Supplementary Figure S2). For An. coluzzii, we observed that HR mosquito samples exhibited a low genetic diversity marked by a reduced haplotype diversity with one haplotype and no polymorphic sites detected (Table 1, Supplementary Figure S3), compared to HS where 23 polymorphic sites were found in 4 haplotypes (h) with a high haplotype diversity (Hd = 0.60) and to the Ngousso strain (h = 3; Hd = 0.44). The haplotype network (Figure 3(a1)) and phylogenetic tree (Figure 3(a2)) showed a predominance of Hap_1 (75%) in the HR mosquitos compared to other groups with a marked difference in the haplotype diversity, suggesting that a possible selection could be acting on this gene or in nearby genes. Higher diversity is observed in the HS group with a high number of mutational steps between haplotypes, with the presence of two Indels of 3 bp and 8 bp. A positive Tajima's D was observed in An. Coluzzii, indicating a decrease in population size and/or balancing selection. In An. gambiae, no difference was observed when comparing the polymorphism parameters of the HR and HS groups, with a similar number of polymorphic sites (23 vs. 20) and a low haplotype diversity, 0.23 and 0.30 in HR and HS, respectively. However, a difference was observed between these groups and the susceptible lab Kisumu where the number of haplotypes is h = 4 with an equivalent high haplotype diversity Hd = 0.60 ( Table 1). The haplotype network (Figure 3(b1)) and phylogenetic tree (Figure 3(b2)) showed that Hap_3 was the most represented (70%) and shared by HR and HS. Other haplotypes, such as Hap_1 (12%) and Hap_2 (6%), were found in the HS group, while Hap_5 was only found in the HR group (Figure 3(b1), Supplementary Figure S3). This Hap_3 has the same sequence as the main haplotype (Hap_1) found in An. coluzzii. The negative values of the neutral test Tajima are indicative of purifying selection in HR and HS groups but are not statistically significant.

Polymorphism Analyses of the Full-Gene Region of CYP6M2
The full-gene of CYP6M2 was successfully sequenced in 21 and 28 individuals of both An. coluzzii (10 HR, 9 HS, 9 Ngousso strain) and An. gambiae (7 HR, 7 HS, 7 Kisumu lab strain), respectively (Supplementary Figure S2). The nucleotide sequence analysis of this fragment exhibited a lower number of haplotypes (8) and polymorphic sites (20) in An. coluzzii HR as compared to 14 haplotypes and 27 polymorphic sites in HS individuals (Table 2, Supplementary Figure S4). This is evident in the lower number of synonymous (11 and 13) and non-synonymous (7 and 10, respectively) mutations but not when compared with lab strain Ngousso ( Table 2). The haplotype network (Figure 4(a1)) and phylogenetic tree (Figure 4(a2)) analysis showed a predominance of Hap_3 (21%), which was shared by all mosquito groups tested. Most haplotypes are separated from each other by 1 to 4 mutational steps and have probably derived from common ancestral susceptible strains.
In An. gambiae, HR mosquitos showed 36 polymorphic sites with 14 non-synonymous mutations compared to 44 polymorphic sites giving all 14 non-synonymous mutations in HS ( Table 2). The haplotype network and phylogenetic profile showed two clear clusters among 20 An. gambiae haplotypes, a resistant cluster with 6 haplotypes and a susceptible cluster with 5 haplotypes (Figure 4(b1,b2)). But, in both species, a star-like shape of the haplotype network indicates a high frequency of intermediate variants, which can be an indication of population expansion.

Polymorphism Analyses of the Full-Gene Region of CYP6M2
The full-gene of CYP6M2 was successfully sequenced in 21 and 28 individuals of both An. coluzzii (10 HR, 9 HS, 9 Ngousso strain) and An. gambiae (7 HR, 7 HS, 7 Kisumu lab strain), respectively (Supplementary Figure S2). The nucleotide sequence analysis of this fragment exhibited a lower number of haplotypes (8) and polymorphic sites (20) in An. coluzzii HR as compared to 14 haplotypes and 27 polymorphic sites in HS individuals (Table 2, Supplementary Figure S4). This is evident in the lower number of synonymous (11 and 13) and non-synonymous (7 and 10, respectively) mutations but not when compared with lab strain Ngousso ( Table 2). The haplotype network (Figure 4(a1)) and phylogenetic tree (Figure 4(a2)) analysis showed a predominance of Hap_3 (21%), which was shared by all mosquito groups tested. Most haplotypes are separated from each other by 1 to 4 mutational steps and have probably derived from common ancestral susceptible strains. In An. gambiae, HR mosquitos showed 36 polymorphic sites with 14 non-synonymous mutations compared to 44 polymorphic sites giving all 14 non-synonymous mutations in HS ( Table 2). The haplotype network and phylogenetic profile showed two clear clusters among 20 An. gambiae haplotypes, a resistant cluster with 6 haplotypes and a susceptible cluster with 5 haplotypes (Figure 4(b1, b2)). But, in both species, a star-like shape of the haplotype network indicates a high frequency of intermediate variants, which can be an indication of population expansion.  Moreover, the GCC-> TCC polymorphism in codon 392-CYP6M2, inducing an amino acid change of alanine to serine A392S (Supplementary Table S4 and Figure S5), was observed in An. gambiae samples, with a frequency of 71.43% (5/7) in the HR and 57.14% (4/7) in HS but not observed in Kisumu lab strain (0/7) as shown in Supplementary  Table S4. Samples carrying the 392A allele were more polymorphic (S = 33, Hd = 0.89 and π = 0.0089; Figure 4(b3)) than 392S allele carriers (S = 25, Hd = 0.88 and π = 0.0059; Figure 4(b3)), suggesting a purifying or balancing selection in mosquitos carrying the A 392 S and confirmed by a positive value obtained for the Tajima's D and Fu and Li statistics (Supplementary Table S4).
Moreover, the GCC-> TCC polymorphism in codon 392-CYP6M2, inducing an amino acid change of alanine to serine A392S (Supplementary Table S4 and Figure S5), was observed in An. gambiae samples, with a frequency of 71.43% (5/7) in the HR and 57.14% (4/7) in HS but not observed in Kisumu lab strain (0/7) as shown in Supplementary Table S4. Samples carrying the 392A allele were more polymorphic (S = 33, Hd= 0.89 and π = 0.0089; Figure 4(b3)) than 392S allele carriers (S = 25, Hd = 0.88 and π = 0.0059; Figure 4(b3)), suggesting a purifying or balancing selection in mosquitos carrying the A 392 S and confirmed by a positive value obtained for the Tajima's D and Fu and Li statistics (Supplementary Table S4).

Polymorphism Analyses of the Upstream Region of CYP6P4
The CYP6P4 upstream region (the intergenic region between CYP6P4 and CYP6P5) was amplified and sequenced for individual HR and HS mosquitos. As above, An. gambiae mosquito samples were used as a model to characterize this intergenic region. After aligning sequences, an indel of 7 bp (GGGGTGC) was consistently observed with a deletion associated with an A/T substitution in all the HR (6/6) and less than 50% in HS, and present in susceptible strain Kisumu. Supplementary Table S5 shows the frequency of this indel across the groups used. Overall, the core elements detected by GPMiner include TATA boxes (5) Table S5).
Regarding the 870-bp of the intergenic region between CYP6P4 and CYP6P5, a comparative analysis was done using 18 sequences (6 HR, 4 HS and 8 Ngousso) and 22 (6 HR, 9 HS and 7 Kisumu), respectively, for An. coluzzii and An. gambiae (Supplementary Figure S7). In An. coluzzii, HR and HS mosquitos revealed a high haplotype diversity of 0.88 and 0.71, respectively ( Table 3). The haplotype 2 was predominant (39%) and shared by all groups, while Hap_3 (22%) and Hap_1 (5%) were specific to susceptible Ngousso lab strain; Hap_4 (5%) and Hap_5 (5%) to HS; and the rest to HR (Figure 5(a1,a2), Supplementary Figure S8). In addition, most HR haplotypes are separated from each other and from other haplotypes by a large number of mutational steps confirming the high genetic diversity observed in their sequences. In An. gambiae, the number of haplotypes for this region was low in HR (2 haplotypes) with a low Hd of 0.485 for 2 polymorphic sites, compared to HS (7 haplotypes) with Hd = 0.85 (Table 3, Supplementary Figure S8) and 17 polymorphic sites, suggesting a possible ongoing directional selection in resistant mosquitos. The haplotype Hap_2 (34%) was shared by all groups with Hap_3 (13%), present only in the susceptible groups (HS and susceptible lab strain Kisumu), while Hap_11 (9%) was only found in the HR groups ( Figure 5(b1,b2)). This Hap_2 corresponds to the main haplotype (Hap_2) found in An. coluzzii. This reduced diversity in An. gambiae HR groups supports the hypothesis of ongoing selection, but it is not confirmed by the positive values of both neutrality tests. In addition, based on the allelic profile at position 273 (Supplementary Table S5), the genetic variability parameters showed that a higher number of haplotypes occurs in 273-A allele carriers (h = 7), with an equivalent high haplotype diversity, Hd = 0.84, and nucleotide diversity, π = 0.0096, while the lowest haplotype diversity and nucleotide diversity (h = 5, Hd = 0.58, π = 0.0013) were found in 273-T allele carriers ( Figure 5(b3), Supplementary Table S5).

Polymorphism Analyses of the Full-Gene Region of CYP6P4
Analysis of variability of a 1051 bp-fragment of the CYP6P4-gene (1521 bp) was done on 58 individual sequences from An. coluzzii (9 HR, 10 HS, 10 Ngousso strain) and An. gambiae (10 HR, 9 HS, 10 Kisumu lab strain) (Supplementary Figure S7). In An. coluzzii, HR mosquitos showed a lower number of haplotypes (8) than the HS individuals (9) (Table 4, Supplementary Figure S9). The haplotype network analysis and phylogenetic tree for this fragment showed that the observed haplotypes are genetically diverse and do not cluster according to resistance phenotype (Figure 6(a1, a2)). However, Hap_2 was pre-

Polymorphism Analyses of the Full-Gene Region of CYP6P4
Analysis of variability of a 1051 bp-fragment of the CYP6P4-gene (1521 bp) was done on 58 individual sequences from An. coluzzii (9 HR, 10 HS, 10 Ngousso strain) and An. gambiae (10 HR, 9 HS, 10 Kisumu lab strain) (Supplementary Figure S7). In An. coluzzii, HR mosquitos showed a lower number of haplotypes (8) than the HS individuals (9) ( Table 4, Supplementary Figure S9). The haplotype network analysis and phylogenetic tree for this fragment showed that the observed haplotypes are genetically diverse and do not cluster according to resistance phenotype (Figure 6(a1,a2)). However, Hap_2 was predominant (14%) and shared by all sample groups tested, the others being distributed as shown in Figure 6(a1,a2) below. The absence of predominant haplotypes indicates that this gene is not subject to selection pressure in An. coluzzii specimens from Yaoundé, Cameroon.
In addition, two synonymous and one non-synonymous mutations were identified (Supplementary Figure S10) in the An. gambiae CYP6P4 gene. The synonymous substitution GCG-> GCA in codon 115 was linked to the mutation of TGC-> AGC in codon 168, which led to an amino acid change of cysteine to serine (C168S) substitution; the two mutations were found with a frequency of 70% in the HR, 28% in HS and 15% in Kisumu (Supplementary Table S6). The carriers of 168C allele show the highest polymorphism with 23 haplotypes (Hd = 0.97, π = 0.0005) in 17 polymorphic sites (Supplementary Table S6) compared to samples carrying the 168S allele where 11 haplotypes (Hd = 0.86, π = 0.0003) (Figure 6(b3) and Supplementary Table S6). The non-synonymous GGC-> GGT at codon 144 was found with a frequency of 95% HR compared to HS and Kisumu of 50% and 20% frequency respectively (Supplementary Table S6).

Distribution of Key Mutations Found in CYP6M2 and CYP6P4 Genes among the Populations of An. gambiae
The presence of mutations in the upstream and the full-gene length regions of CYP6M2 and CYP6P4 were in the HS and HR strains, and direct field-collected samples from Nkolondom (Table 5, Supplementary Figures S11 and S12) were determined. No significant link between any of the mutations was identified, and pyrethroid resistance was detected (Table 5).  In addition, two synonymous and one non-synonymous mutations were identified (Supplementary Figure S10) in the An. gambiae CYP6P4 gene. The synonymous substitution GCG-> GCA in codon 115 was linked to the mutation of TGC-> AGC in codon 168, which led to an amino acid change of cysteine to serine (C168S) substitution; the two mutations were found with a frequency of 70% in the HR, 28% in HS and 15% in Kisumu (Supplementary Table S6

Discussion
Insecticide resistance is a major threat to the effectiveness of insecticide-based malaria control tools. Besides the already well-studied target site mutations, metabolic resistance is one of the most important mechanisms. Although the expression of cytochrome P450 genes, such as CYP6M2 and CYP6P4, capable of breaking down the insecticide before it reaches its target is a well-known mechanism of insecticide resistance, the underlying factors modulating the expression of this enzyme family are poorly understood. The present study analyzed the polymorphisms of these two cytochrome-P450s to detect potential mutations that are linked to their overexpression and may be used as resistance markers in two main malaria vectors populations in Yaoundé, Cameroon.
The study shows resistance to permethrin in both species in line with previous reports in different studies in the country [23][24][25]57,58] and highlighted an escalation of resistance over years in the country with observed resistance to higher insecticide concentrations [32]. Additionally, pre-exposure to PBO reveals a partial recovery of susceptibility to permethrin in both collection sites, indicating an involvement of oxidases, such as cytochromes P450 in the observed resistance. However, the recovery of susceptibility after PBO is moderate (<20%) showing that other mechanisms than P450s are involved, including kdr, which were found fixed in An. gambiae. The escalation in resistance could be explained by the insecticide selection pressure due to the scale up of control interventions with the distribution of ITNs taking place in Cameroon since 2011 and also by the presence of pesticide residues or components associated with anthropogenic activities (soapy water, organic waste) in the urban breeding sites. Those residues may activate certain detoxication enzymes in mosquitos by oxidative stress, enzymes being also active against insecticides molecules [59,60].
The consistent over-expression of CYP6M2 and CYP6P4 genes in field permethrinresistant An. gambiae and An. coluzzii, among all the up-regulated genes found in this study, shows their link with resistance patterns in Cameroon. These two genes have been firstly recorded upregulated in Cameroon DDT and pyrethroid resistant An. gambiae and An. coluzzii mosquitos [22], and more recently, in An. gambiae populations from the same peri-urban area [23,24]. Meanwhile, Ibrahim et al. [14] have already established that CYP6P4 is responsible for resistance to permethrin in An. arabiensis populations from Central Africa. Many studies have shown that the up-regulation of certain P450s can be modulated by mutations of cis-regulatory elements in the promoter regions of P450 genes and/or changes in the expression level of transcription factors binding to these cis-regulatory elements [28,61,62]. Knowing that other mechanisms, such as allelic variations of resistance genes, may also be involved, we decided in this study to amplify around 1 kb of the putative promoter and full-gene region of CYP6M2 and CYP6P4 in both An. coluzzii and An. gambiae to detect any variations that might influence their expression in resistant individuals.
Polymorphism analysis of CYP6M2 sequences in An. coluzzii from Ngousso, Yaoundé, revealed that for both the upstream and the full-gene length regions, the HS strain was more polymorphic than the HR, but the predominant haplotype was shared by both HR and HS. In An. gambiae the upstream region of CYP6M2 in HS mosquitos showed a similar polymorphism to that in HR, and this, coupled with the high diversity in HR samples compared to HS, suggests that there is not a strong selective sweep acting upon CYP6M2 in either species in Yaoundé. The absence of polymorphisms associated with resistance in the coding region or in the 1 kb upstream may indicate that the emergence of CYP6M2associated resistance is due to selection pressures acting on genes encoding distantly related regulatory proteins in trans position, as concluded in a similar study comparing CYP6M2 sequences of An. gambiae s.l. sampled from 13 countries [18]. This larger study identified 6 amino acid polymorphisms within CYP6M2, none of which were detected in the current study. The mutation A392S was found to be the most frequent in F 4 HR (50%) of An. gambiae from Nkolondom crossing and was also detected in field F 0 populations (7%).
A detailed analysis of the polymorphisms in both the upstream and full-gene regions of CYP6P4 revealed that the HS was more polymorphic than HR samples of An. gambiae populations. The key amino acid change (C168S) observed in CYP6P4 and substitution mutations in the nrf2/MAF TFBs of their upstream region appear to be subject to positive selection in An. gambiae populations and could signal a directional selection as seen for CYP6P9a/b [45].
In the full-gene region, one major amino-acid change at codon 168 where cysteine (C) is replaced by serine (S) were found in almost all the HR sequences. This amino-acid change was different from what has been seen in Uganda with the I 236 M mutation [20], with only the susceptible allele 236I found in all the An. coluzzii and An. gambiae sequence from Yaoundé, Cameroon. The impact of such allelic variation on the metabolic efficiency of detoxication genes has previously been demonstrated in An. funestus for P450s, such as CYP6P9a/b [27], and is similar to the case of CYP6A2 in Drosophila melanogaster for which 3 amino acid substitutions located close to the active site in the allele predominant in DDT-resistant flies have been shown to confer the increased metabolism of DDT [62]. Further studies on the variants are needed to determine whether they vary in their ability to detoxify pyrethroids.

Conclusions
This study explored the genetic diversity of two cytochrome P450s associated with pyrethroid resistance in An. gambiae s.l. revealing some potentially interesting molecular markers in CYP6P4 that merit further investigation in the resistance context in Cameroon. No molecular markers associated with resistance were found in CYP6M2 in either species agreeing with previous studies. Further studies need to be performed to detect the specific genetic variants driving the overexpression of both CYP6M2 and CYP6P4 in the central African population of An. gambiae s.l., exploring structural variants and both cis-and trans-regulatory loci. Such work will help design robust DNA-based assays to detect and track resistance.  Figure S11: Representative diagram of DNA-based assay to genotype a key mutation in An. gambiae CYP6M2. (a) CYP6M2 upstream region: alignment of sequences showing differences by resistance phenotype including the deletion of 7 bp found in HR and HS groups link to the G/A variant generating a restriction site for the BsrDI restriction enzyme; a schematic representation of the CYP6M2 promoter PCR-RFLP illustration digestion of the PCR amplicon and genotyping results for F4 field-resistant from Nkolondom and Susceptible Kisumu crossing. (b) CYP6M2-gene region: alignment of sequences showing A392S-mutation differences by resistance phenotype; a schematic representation of the CYP6M2-gene AS-PCR illustration digestion of the PCR amplicon and genotyping results for F4 fieldresistant from Nkolondom and Susceptible Kisumu crossing; Figure S12: Representative diagram of DNA-based assay to genotype a key mutation in An. gambiae CYP6P4. (a) CYP6P4 upstream region: alignment of sequences showing differences by resistance phenotype linked to the A/T (A-273-T) variant generating a restriction site for the PvuII restriction enzyme; a schematic representation of the CYP6P4 promoter PCR-RFLP illustration digestion of the PCR amplicon and genotyping results for F4 field-resistant from Nkolondom and Susceptible Kisumu crossing. CYP6P4 gene region: (b) alignment of sequences showing the differences in codon 144 (C-432-T) by resistance phenotype; a schematic representation of the CYP6P4-gene RFLP-PCR illustration digestion of the PCR amplicon and genotyping results for F4 field-resistant from Nkolondom and Susceptible Kisumu crossing; (c) alignment of sequences showing C168S-mutation differences by resistance phenotype; a schematic representation of the CYP6P4-gene AS-PCR illustration digestion of the PCR amplicon and genotyping results for F4 field-resistant from Nkolondom and Susceptible Kisumu crossing; Table  S1: Primers of the promoter region of CYP6M2 and CYP6P4; Table S2: Specific primers of CYP6M2 and CYP6P4 diagnostic assays; Table S3: Genotype and allele frequencies of 1014F kdr mutations in An. gambiae and An. coluzzii populations; Table S4: Frequency of key mutations found and genetic variability parameters of CYP6M2 per allele in An. gambiae; Table S5: Frequency of key mutations found and genetic variability parameters of the intergenic region between CYP6P4 and CYP6P5 per allele in An. gambiae; Table S6: Frequency of key mutations found and genetic variability parameters of the upstream region of mutation per allele in An. gambiae.  Data Availability Statement: All data generated or analyzed during this study are included within the article and its additional files. The sequences generated have been deposited in the GenBank database (study accession numbers: OP779727-OP779746; OP797835-OP797861; OP797862-OP797886; OP811275-OP811307; OP846110-OP846114; OP846115-OP846119; OP857148-OP857158; OP857159-OP857169).