Esdras Mahoutin Odjo^1,2*, Armel Djènontin^1,2, Aboubacar Sidick¹, Arsène Jacques Fassinou¹, Albert S. Salako¹, Brice A. Fanou⁴, Constantin J. Adoha^1,2, Boulais Yovogan^1,2, Christian S. T. Akpodji^1,2, Zinsou C. Koukpo¹, Antoine Abel Missihoun², Arthur Sovi^1,5,6, Casimir D. Kpanou¹, Linda Towakinou¹, Razaki Osse^1,7, Roseric Azondékon¹, Gil Germain Padonou^1,2, Martin Akogbéto^1†, Clement Agbangla^1,3†

¹Centre de Recherche Entomologique de Cotonou, Cotonou, Bénin.

²Faculté des Sciences et Techniques - Université d’Abomey-Calavi, Abomey-calavi, Bénin.

³Direction Générale de la Recherche Scientifique, Ministère de l’Enseignement Supérieur et de la Recherche Scientifique, Cotonou, Bénin

⁴Ecole Polytechnique d’Abomey Calavi, Université d’Abomey-Calavi, Abomey-Calavi, Bénin

⁵Faculty of Infectious and Tropical Diseases, Department of Disease Control, The London School of Hygiene and Tropical Medicine, London WC1E 7HT, UK.

⁶Faculté d’Agronomie- Université de Parakou, Parakou, Benin.

⁷Université Nationale d’Agriculture de Porto-Novo, Porto-Novo, Bénin

*Corresponding Author: Esdras Mahoutin Odjo, Centre de Recherche Entomologique de Cotonou, Cotonou, Bénin.

Abstract
From 2019 to 2021, pirimiphos-methyl, clothianidin + deltamethrin, and clothianidin alone were successively used to protect communities in the Donga and Alibori departments against mosquito bites. This study aims to analyze the genetic structure of Anopheles coluzzii and Anopheles gambiae over three years of intervention.

The vector population studied comes from adult mosquitoes collected monthly from 2019 to 2021 in the study communes. Mosquitoes were collected on the night of each visit in four houses per village and two villages per commune. Molecular analysis of populations was carried out on 3,134 sampled anopheles.

The frequencies of kdr-L995Fmutation were similar between the treated communes and the control commune, suggesting constant pressure from the pyrethroids and DDT used in agriculture. Anopheles coluzzii showed relatively lower frequencies than Anopheles gambiae for kdr- L995F mutation, linked to its preference for permanent breeding sites, which are less exposed to insecticides during the dry season. The ace-1- G280S mutation shows low frequencies in both species’ populations. Genetic differentiation is high for the ace-1 gene. Genetic differentiation varies for the voltage-gated sodium channel (vgsc) gene, suggesting variable gene flow between sub-populations.

The clothianidin + deltamethrin mixture could reduce the frequency of kdr-L995F. Genetic differentiation is very significant and similar between the ace-1 (G280S) gene sub-populations.

Keywords: Genetic differentiation, kdr-L995F, ace-1-G280S, Anopheles coluzzii, Anopheles gambiae, IRS, clothianidin, Benin

Introduction
Background
Efforts to combat malaria continue to face a convergence of challenges, mainly in Africa, which carries the heaviest burden of this devastating disease [1]. The COVID-19 pandemic, along with other humanitarian crises, pressures on health systems, budgetary constraints, rising biological threats, and the declining effectiveness of crucial malaria control strategies, are severely hampering progress towards global malaria control targets [1]. Malaria vectors are becoming increasingly resistant to insecticides used in public health for net impregnation and indoor residual spraying (IRS) [2].

In Benin, the species Anopheles gambiae and Anopheles coluzzii are the main vectors that dominate malaria transmission [3]. To control them in the northern region, where malaria incidence is high, the national malaria control program (NMCP) has turned to indoor residual spraying (IRS) alongside the use of long-acting insecticidal nets (LLINs) [4]. IRS involves spraying the inside walls of homes with residual insecticides to eliminate mosquito vectors and reduce malaria transmission [5]. From 2019 to 2021, three types of insecticide were successively used for large-scale community IRS in the departments of Donga and Alibori in northern Benin [4]. The used insecticides included pirimiphos-methyl 300 CS, a combination of clothianidin 500g/kg with deltamethrin 62.5 g/kg, and clothianidin 50 WG on its own. This periodic change of insecticides with different modes of action was aimed at optimizing the efficacy of IRS, which can be influenced by several factors, including the genetic variability of Anopheles populations.

This study aims to examine the genetic makeup of Anopheles coluzzii and Anopheles gambiae populations in the northern Benin region, specifically investigating the potential impact of IRS on the genetic diversity of these vector populations. The findings from these analyses will offer essential insights for enhancing the efficiency of malaria vector control methods.

Methods
Study area
Bembèrèkè commune in Borgou was selected as a control in comparison to the IRS-treated communes of Donga (Djougou and Copargo) and Alibori (Kandi and Gogounou) (Figure 1)

The Department of Borgou is bordered to the north by the Department of Alibori and to the southeast by the Department of Donga (Figure 1). Like Donga and Alibori, Borgou has a dry and rainy season, with annual rainfall between 900 and 1,300 mm. The average annual temperature in Borgou is around 26°C, and relative humidity varies between 30% and 70%. The soils are tropical ferruginous, ferritic, sandy clay, or sandy clay, and gneissic granite [6]. Borgou is a highly agricultural department. Around 66% of its population is involved in agriculture, covering 54% of the department's total surface area (13,962 km² / 25,856 km²).

Around 69% of households in the Alibori department are agricultural (74693/108351), corresponding to around 14% of Benin's farming population, making it the breadbasket of Benin [7]. Rainfall varies between 700 mm and 1,200 mm. The flat terrain, which is sometimes shaped by a Cretaceous sedimentary series or crowned by armored hillocks sloping down to the river Niger and hills of ferruginous sandstone, features sparse shrub savannah and gallery forests bordered by watercourses. The primary soils are of the ferruginous type on crystalline bedrock, very fertile alluvial soils in the Niger Valley and clayey, black silty soils in the lowlands, swamps, and very fertile gallery forests, where rice, market gardening, and yam are grown.

The average rainfall in Donga is between 1,200 mm and 1,300 mm. During the rainy season, the watercourses flood the low-lying areas, making them ideal for rice growing. The soils are crude mineral, tropical ferruginous, indurated, and hydromorphic. Gallery forests are formed by dense vegetation along watercourses. 59% of households in Donga were farmers (39,461/66,433), and 96.8% of farming households were involved in crop production and 3% in animal production [8].

Figure 1: Map of the study communes

Collection of mosquito samples
The Anopheles gambiae s.l. The population studied came from adult mosquitoes collected monthly in 2019, 2020, and 2021 by capture from human volunteers in the control commune (Bembèrèkè) and the IRS communes (Gogounou and Kandi to the north of Bembèrèkè, Djougou, and Copargo to the southeast of Bembèrèkè).

Two villages were chosen in each commune: one central and one peripheral. Local volunteers captured adult mosquitoes in four houses per visit, with one collector inside and one outside each household, conducting hourly collections. However, no insecticide-resistance bioassays were performed on the collected mosquitoes [9, 10].

Molecular characterization of vector
Collected mosquitoes were identified and separated at the species level based on morphological criteria according to established taxonomic keys [11, 12].

3,134 Anopheles gambiae mosquitoes collected from the field were subjected to molecular characterization. Genomic DNA extraction was performed on the mosquito's body, including the abdomen, wings, and legs. The molecular biology protocol outlined by Myriam and Cécile [13] was employed for this stage.

Species identification
The genomic DNA extracted was used for molecular identification of the species within the Anopheles gambiae complex.

All mosquitoes were subjected to PCR using the Scott et al. [14] protocol to identify the different species of the Anopheles gambiae complex. The primers used below should identify the following species: Anopheles gambiae s.s. (AG), Anopheles arabiensis (AA), Anopheles melas (ME), Anopheles merus (ME), Anopheles quadriannulatus A and B (QD).
UN: GTGTGCCGCTTCCTCGATGT
AG: CTGGTTTGGTCGGCACGTTT
AA: AAGTGTCCTTCTCCATCCTA
ME: TGACCAACCCACTCCCTTGA
QD: CAGACCAAGATGGTTAGTAT

PCR was performed in a 30-cycle program. Denaturation at 94 °C for 30 s, hybridization at 50 °C for 30 seconds, and elongation at 72 °C for 30 seconds. The amplified DNA copies were kept at a final temperature of 4 °C before migration on a 2.5% agarose gel with ethidium bromide used as an intercalating agent.

The technique of Santolamazza et al. [15] was used to distinguish the twin species Anopheles gambiae and Anopheles coluzzii. The PCR program includes an initial denaturation at 94 C for 5 minutes, followed by 35 cycles. Each cycle involves denaturation at 94 °C for 30s, hybridization at 54 °C for 30s, and elongation at 72 C for 30 s. A final elongation at 72 °C for 10 min is performed to allow complete amplification of the sequences. The following primers were used: 200X6.1FTCGCCTTAGACCTTGCGTTA and 200X6.1R-CGCTTCAAGAATTCGAGATAC.

The PCR products are maintained at 4°C before undergoing migration through 1.5% agarose gel electrophoresis, utilizing ethidium bromide as an intercalating agent.

Identification of the L995F mutation within the kdr gene
The presence of the resistance allele (L995F) of the kdr gene in samples collected from each study site was detected by PCR following the protocol described by Martinez-Torres et al. [16]. The following primers were used.

Agd1: 5′-ATAGATTCCCCGACCATG-3′;

Agd2: 5′-AGACAAGGATGATGAACC-3′;

Agd3: 5′ ATTTGCATTACTTACGACA-3′;

Agd4: 5′-CTGTAGTGATAGGAAATTTA-3′.

The amplification protocol consists of 40 cycles, each comprising an initial denaturation at 94°C for 1 minute, hybridization at 48°C for 2 minutes, and elongation at 72°C for 2 minutes. The PCR concludes with a final elongation step at 72°C for 10 minutes [16].

Identification of the G280S mutation within the ace-1 gene.
The Weill et al. protocol [17] was employed to detect the G280S mutation in mosquitoes. The PCR was conducted using the following primers:
Moustdir1 5′CCGGGNGCSACYATGTGGAA3′ and Moustrev1 5′ACGATMACGTTCTCYTCCGA3′.
The amplification program consisted of thirty cycles, comprising denaturation at 94°C for 30 seconds, hybridization at 52°C for 30 seconds, and elongation at 72°C for 1 minute. Subsequently, PCR products underwent digestion with the AluI restriction enzyme following the manufacturer's guidelines before being migrated onto a 2% agarose gel.

Statistical analysis
To calculate the allelic frequencies of the kdr-L995F and ace-1 G280S genes in each commune and by species, stratified 3×3×3 and 2×3×3 contingency tables were done in IBM-SPSS Statistics Subscription® (Build 1.0.01406). The Pearson chi-square (χ2) test was conducted to compare proportions. The allele frequencies of kdr-L995F and ace- 1R G280S were calculated using the formula:

Where F® is the frequency of resistance, n is the number of mosquitoes of a given genotype, RR is the homozygous resistant genotype, RS is the heterozygous resistant genotype, and SS is the susceptible genotype [18].

FIS Fixation Index
This parameter quantifies the variance between the population of individuals identified in the HO heterozygous state and the expected heterozygous (HE) rate. It is also called the pan mixing deviation and is calculated by the formula: FIS = 1 - (HO / HE) at the specified locus within a population of diploid individuals. This index represents the departure from the Hardy-Weinberg equilibrium. It varies from -1 to +1 and indicates the heterozygote deficit per population, per locus, and for all loci. FIS is positive when the population has a deficit of heterozygotes compared with the panmictic equilibrium and negative in the opposite case. If the assumptions of the Hardy-Weinberg model are met, it will be equal to HE. When the observed heterozygosity falls below the expected heterozygosity (HO < HE), it indicates an overrepresentation of homozygotes compared to anticipated in the studied populations. This excess of homozygotes can then be assimilated to a risk of inbreeding, drift, selection, and differentiation within populations.

Genetic differentiation (FST) of populations
The Hardy-Weinberg Equilibrium was assessed for each population using Genetics software version 4.7.5. The fixation index (FIS) was computed using Genepop software, following the method described by Weir and Cockerham [19].

An FST value nearing 0 indicates substantial genetic interchange among populations, implying minimal genetic differentiation and a panmictic population. Conversely, an FST value approaching 1 signifies significant genetic divergence between populations, indicating limited or no gene flow. According to Wright, an FST between 0 and 0.05 indicates little differentiation; an FST between 0.05 and 0.15 indicates moderate differentiation; an FST between 0.15 and 0.25 suggests significant differentiation, and above 0.25, the FST illustrates very significant differentiation [20].

The various communes concerned were considered sub-populations, and the genetic differentiation of the population (FST) was evaluated between the control commune and each IRS commune by year and by commune according to the years of the IRS campaign. The indices of genetic differentiation within populations (FST) were calculated using Genepop version 4.7.5 software.

Results
Allele frequencies of kdr-L995F mutation
Overall, analysis of these data reveals that the allelic frequency of the kdr-L995F mutation is relatively high in the whole population and varies according to the anopheline populations and the insecticide used. Although the insecticides used for spraying are not directly responsible for the kdr-L995F mutation, an increase in the frequency of this mutation has been observed in specific populations. (Figure 2).

Figure 2: Bar chart showing the kdr allele frequency (percentage with standard error bars in the first panel) in An. gambiae (A) and An. coluzzii (B) in the communes of the Alibori and Donga departments of Benin by species and IRS period PM: pirimiphos-methyl; CS: capsule suspensions; Clot + Del: mixture clothianidin 500 g/kg + deltamethrin 62.5 g/kg; Clo: clothianidin 50 WG; WG: Water dispersible granules

For Anopheles gambiae
The frequency of the kdr-L995F mutation varied across different years and insecticide treatments. In 2019, when pirimiphos-methyl 300 CS was used, the frequency ranged from 82% to 89%. For the mixture of clothianidin 500 g/kg and deltamethrin 62.5 g/kg in 2020, the frequency ranged from 70% to 94%. In 2021, when clothianidin 50 WG was used alone, the frequency ranged from 80% to 88% (Figure 2A). For each insecticide, no remarkable difference was observed between the permutation frequencies of the treated communes and those of the control commune. Anopheles gambiae showed a kdr-L995F mutation frequency of 89% with pirimiphos- methyl 300 CS in Djougou, compared with 77% with the Clothianidin+ Deltamethrin mixture. Apart from Gogounou, the clothianidin 500 g/kg + deltamethrin 62.5 g/kg mixture appears to have reduced the frequency of the permutation in the communes treated.

For Anopheles coluzzii
The magnitude of the kdr-L995F mutation varied across different insecticide treatments. For pirimiphos-methyl 300 CS, the frequency ranged from 65% to 78%, while for the mixture of clothianidin 500 g/kg and deltamethrin 62.5 g/kg, it ranged from 45% to 60%. Clothianidin 50 WG alone exhibited frequencies ranging from 58% to 95%. Notably, Anopheles coluzzii generally displayed a lower prevalence of the kdr-L995F mutation than Anopheles gambiae across most regions (Figure 2).

Allele frequencies of ace-1 G280S mutation
The data in Figure 3 shows the allelic frequencies of the G280S mutation of the ace-1 gene in Anopheles gambiae and Anopheles coluzzii in the control commune (Bembèrèkè) and in the IRS communes (Kandi, Gogounou, Djougou and Copargo) exposed to various insecticides. The mutation occurred infrequently across all the study communes, with a maximum value of 5%.

For Anopheles gambiae
Frequencies varied from one commune to another. The allelic frequency was 2% to 3% with pirimiphos-methyl 300 CS, 1% to 5% with the clothianidin + deltamethrin mixture, and 3% to 4% with Clothianidin alone. No notable difference (p > 0.05) was detected between the average allele frequency of the IRS communes and that of the control commune (Figure 3A).

Figure 3: Bar chart showing the ace-1 allele frequency (percentage with standard error bars in the first panel) in An. gambiae (A) and An. coluzzii (B) in the communes of the Alibori and Donga departments of Benin by species and IRS period PM: pirimiphos-methyl; CS: capsule suspensions; Clot + Del: mixture clothianidin 500 g/kg + deltamethrin 62.5 g/kg; Clo: clothianidin 50 WG; WG: Water dispersible granules

For Anopheles coluzzii
The mutation frequency was 0% to 3% with pirimiphos-methyl 300 CS, 1% to 3% for the clothianidin and deltamethrin mixture, and for clothianidin 50 WG alone.

FIS fixation index
The deviation from panmixis within the different populations was assessed (Tables 1 and 2). For the kdr L995F mutation of the vgsc gene, the Anopheles gambiae and Anopheles coluzzii populations showed a positive fixation index (FIS) (FIS > 0). Expected heterozygosity (HE) was higher than observed heterozygosity (HO). Fixation index values varied within populations between IRS campaigns from 2019 to 2021 with different insecticides. The minimum and maximum values were found in Kandi (0.35 < FIS <0.65) for Anopheles gambiae and in Kandi (FIS = 0.11) and Djougou (FIS= 0.57) for Anopheles coluzzii. These positive fixation indices reflected a heterozygosity deficit within the populations (Tables 1 and 2).

For the G280S mutation in the ace-1 gene, the fixation index (FIS) was negative overall in Anopheles gambiae and Anopheles coluzzii populations. In 2020, the Anopheles gambiae populations of Gogounou and Copargo tended to return to equilibrium under the Hardy-Weinberg hypothesis (FIS ≈ 0). The same is true for the different populations of Anopheles coluzzii from 2019 to 2020. The FIS values were low, showing an incipient excess of heterozygotes in the populations. (Tables 1 and 2).

Table 1. Observed heterozygosity (HO), expected heterozygosity (HE), fixation index (FIS) and number of individuals sampled of An. gambiae per population and per insecticide.

N: number of sampled individuals; HO: observed heterozygosity; HE: expected heterozygosity; FIS: fixation index; PM: pirimiphos-methyl; CS: capsule suspensions; Clot + Del: mixture clothianidin 500 g/kg + deltamethrin 62.5 g/kg; Clo: clothianidin 50 WG; WG: Water dispersible granules.

Table 2. Observed heterozygosity (HO), expected heterozygosity (HE), fixation index (FIS) and number of individuals sampled of An. coluzzii per population and per insecticide.

N: number of sampled individuals; HO: observed heterozygosity; HE: expected heterozygosity; FIS: fixation index; PM: pirimiphos-methyl; CS: capsule suspensions; Clot + Del: mixture clothianidin 500 g/kg + deltamethrin 62.5 g/kg; Clo: clothianidin 50 WG; WG: Water dispersible granules; NA: Not Applicable.

Genetic differentiation FST
Pairwise FST was calculated for both loci, and the results were compared between the control commune and each IRS commune for each insecticide (Table 3) and each commune according to insecticides (Table 4).

The FST values in Tables 3 and 4 indicate the degree of genetic differentiation between population pairs for two different loci (kdr- L995F and ace-1 G280S) over three different years (2019, 2020, 2021).

Anopheles gambiae showed significant genetic differentiation at the kdr-L995F locus (Tables 3 and 4) in 2019 and 2021 (FST = 0.24 and 0.70), followed by weak genetic differentiation in 2020 (FST = 0.02) between Bembèrèkè (control population) and Kandi (IRS population). Between Bembèrèkè and Gogounou, genetic differentiation was very high in 2019 (FST = 0.80) and moderate in 2020 and 2021 (FST = 0.09 and 0.20). Bembèrèkè, along with Djougou and Copargo, showed moderate genetic differentiation in 2019 (FST = 0.10 and 0.10), high differentiation in 2020 (FST = 0.32 and 0.35), then very high differentiation in 2021 for Djougou (FST = 0.93) and low differentiation for Copargo (FST = 0.03) (Table 3).

The FST values of Anopheles gambiae for the ace-1 G280S locus between Bembèrèkè and each IRS commune were globally between (FST = 0.6 and 1.0), indicating a significant genetic differentiation between the control commune and each IRS commune.

Overall, for Anopheles gambiae, genetic differentiation is more significant for the ace-1 G280S locus than for their L995F locus.

Anopheles coluzzii for the kdr-L995F locus (Tables 3 and 4) shows variations in genetic differentiation coefficients from one year to the next and between pairs of populations. Overall, the differentiation coefficients decreased from 2019 to 2021, resulting in very high genetic differentiation (FST > 0. 31) to significant genetic differentiation between Bembèrèkè and Gogounou with pirimiphos- methyl and Clothianidin alone (FST = 0.24, FST = 0.19) on the one hand, and between Bembèrèkè and Copargo (FST = 0.19) with the clothianidin + deltamethrin mixture on the other. Moderate genetic differentiation (FST = 0.09) was observed between Kandi and Bembèrèkè with the clothianidin + deltamethrin mixture. Weak genetic differentiation (FST < 0.0001) was observed with Clothianidin alone in 2021 between Bembèrèkè and Djougou and Bembèrèkè and Copargo.

For the ace-1 G280S locus, the genetic differentiation data between different pairs of populations show generally significant differentiation between populations and treatments.

Table 3. Genetic differentiation between An. gambiae and An. coluzzii populations in IRS communes and the control commune for the kdr-L995F and ace-1 G280S mutations

Bem : Bembèrèkè ; Kan : Kandi ; Gog : Gogounou ; Djo : Djougou ; Cop : Copargo ; PM : pirimiphos-methyl; CS: capsule suspensions; Clot + Del: mixture clothianidin 500 g/kg + deltamethrin 62.5 g/kg; Clo: clothianidin 50 WG; WG: Water dispersible granules; NA: Not Applicable.

Table 4. Genetic differentiation between An. coluzzii and An. gambiae populations in IRS and control communes for the kdr-L995F and ace-1 G280S mutations

P-methyl: pirimiphos-methyl; CS: capsule suspensions; Clot + Del: mixture clothianidin 500 g/kg + deltamethrin 62.5 g/kg; Clo: clothianidin 50 WG; WG: Water dispersible granules.

Discussion
The frequencies of the kdr-L995F mutation calculated in all the Anopheles gambiae and Anopheles coluzzii populations showed no remarkable difference between the treated communes and the control commune. This observation is thought to be due to the use of pyrethroid and organochlorine (DDT) insecticides in agriculture in these study areas [21–23]. These classes of insecticides target the voltage-gated sodium channel (vgsc) gene, which contains the kdr- L995F mutation [24, 25]. This heavy use of insecticides is thought to be at the origin of the selection of resistant individuals, thus increasing the frequency of the mutation in the various populations [26, 27]. Several studies have shown similar results in Benin and in several other countries in Africa and elsewhere [25, 27–31].

Furthermore, the frequencies of this mutation observed in Anopheles coluzzii populations are relatively lower than those of Anopheles gambiae over the years. This finding may be linked to the permanent roosting preferred by Anopheles coluzzii, which is dominant during the dry season [32, 33]. During this season, insecticide pressure is low without agricultural activities. However, the high allele frequency in Anopheles gambiae populations is linked to the use of insecticides, particularly pyrethroids [34, 35], which provide crop protection and, in turn, encourage the selection of resistance in these Anopheles gambiae populations [23, 30, 36]. Concerning the ace-1gene, the frequency of the G280S mutation was low in all Anopheles coluzzii and Anopheles gambiae populations during the three years of the study. These low allelic frequencies obtained in these populations are thought to be linked to the momentary use of these insecticides only during IRS, which does not favor the selection of resistant individuals. This low frequency of the resistant allele could also be linked to the high adaptation cost of the ace-1 mutation [37–40]. Previous studies have also shown the low frequency of this mutation in Benin [18, 41]. Several factors could explain the differences observed in panmixia. For the vast L995F gene, all Anopheles gambiae and Anopheles coluzzii populations showed positive FIS values, indicating a heterozygosity deficit (FIS > 0). On the one hand, selection favors homozygous resistant individuals, and on the other, inbreeding is due to preferential mating between individuals of the same genotype. Heterozygous individuals have a lower probability of survival or reproduction than homozygous individuals due to insecticide pressure. These observations were confirmed by the decrease in the frequency of heterozygotes observed in the population [42]. Concerning the ace-1 G280S gene, the fixation index was globally negative (FIS < 0) in Anopheles gambiae populations and mostly in Anopheles coluzzii populations, indicating an excess of heterozygosity. This excess of heterozygosity could be explained by the high genetic cost resulting in a drop in the frequency of the resistant allele. On the one hand, this could be due to natural selection, which is unfavourable to individuals carrying the mutation, resulting in a reduction in the transmission of the resistant allele to future generations. On the other hand, the low frequency of this mutation in populations may also be associated with this phenomenon despite the many years pirimiphos-methyl has been used in IRS. Other environmental factors and the dynamics of selection cannot be ruled out [43–45]. In some Anopheles coluzzii populations, the fixation index was substantially equal to zero, indicating a population at Hardy-Weinberg Equilibrium for the ace-1 locus despite the use of pirimiphos-methyl. This result could be due to the higher genetic cost in Anopheles coluzzii populations compared with Anopheles gambiae populations. These results could also be associated with a second- species risk of error.

Genetic differentiation is extreme in all populations of Anopheles gambiae and Anopheles coluzzii for the ace-1 gene. This result can be explained by the high genetic cost of this gene, which limits genetic flow between the different populations. This strong differentiation could also be explained by the distance separating the various sub- populations. Other ecological or behavioral factors that have yet to be elucidated could also explain these results. With regard to genetic differentiation at the level of the vast gene, there was weak, moderate, and strong differentiation depending on the sub-population pairs considered. The low genetic differentiation observed in 2020 between Kandi and Bembèrèkè (FST = 0.019) and in 2021 between Copargo and Bembèrèkè (FST = 0.026) for Anopheles gambiae suggests high gene flow between the sub-populations. It would also be the result of the absence of a genetic barrier, which would lead to random mixing of alleles.

The same findings were observed in 2021 for the Anopheles coluzzii sub-populations of Djougou-Bembèrèkè and Copargo-Bembèrèkè (FST < 0.001). After the IRS campaigns with the clothianidin + deltamethrin mixture, there was an overall decrease in the allelic frequency of the kdr-L995F mutation for the Anopheles coluzzii and Anopheles gambiae populations. This decrease in kdr-L995F frequency is accompanied by little genetic differentiation in the Anopheles gambiae and Anopheles coluzzii populations. This suggests that the use of the clothianidin + deltamethrin mixture in IRS would limit the spread of the resistant allele of the vgsc gene in Anopheles gambiae and Anopheles coluzzii. Additional research will enhance our comprehension of how the clothianidin + deltamethrin mixture influences the development of the kdr-L995F mutation in various subspecies of the Anopheles gambiae complex.

Conclusion
The lower frequency of the kdr-L995F mutation in Anopheles coluzzii, compared with Anopheles gambiae, is attributed to the preference of Anopheles coluzzii for permanent breeding sites, dominant during the dry season, characterized by low insecticide pressure in the absence of agricultural activities. Conversely, the elevated frequency observed in Anopheles gambiae is linked to the extensive application of pyrethroids and DDT, which are used for crop protection, thus promoting the emergence of resistance. The low frequency of the G280S mutation in populations of Anopheles coluzzii and Anopheles gambiae is explained by the momentary use of carbamates and organophosphates and by the high cost of adaptation of the G280S mutation. This high genetic cost could be behind the excess heterozygosity observed in specific populations, reducing the transmission of the resistant allele to future generations. The fixation index reveals a heterozygosity deficit for the vgsc-L995F gene, suggesting selection favoring homozygous resistant individuals and inbreeding. The high genetic differentiation for the ace-1 gene indicates limited gene flow between populations, probably due to the high genetic cost of this gene. The clothianidin + deltamethrin mixture could have a reducing impact (yet to be elucidated) on the frequency of the kdr-L995F mutation in Anopheles coluzzii and Anopheles gambiae populations.

Conflict of interest: The authors state that they do not have any conflicts of interest.

Data archiving
The data on the genotypes of the various individuals by species and by population and their analysis are available and can be consulted by the public. If you have any questions, please write to the corresponding author.

Funding:
This study was supported by the US President's Malaria Initiative (PMI) through funding for the PMI Indoor Residual Spraying (AIRS) project, the PMI Vector Link project and the PMI Evolve project.

Acknowledgments
We would like to thank the police and administrative authorities of the communes of Djougou, Copargo, Kandi, Gogounou and Bembèrèkè who facilitated the collection of adult mosquitoes in their localities.

Authors' contributions
EMO, MA, AD, and CA designed the study. EMO, AD, CA, AJF, RA, BAF, CJA, BY, AAM, AS, RO, ASS, GGP and MA critically revised the manuscript. EMO, GGP, ASS, ASi, AAM, AJF, BAF, CJA BY, RO, ZCK, CSTA, CDK, LT and CA carried out the field activities and the laboratory analysis. EMO, RA, ZCK and AJF analyzed the data. EMO drafted the manuscript. All authors read and approved the final manuscript.

Ethics approval
The study protocol was reviewed and approved by the full membership of the IECC (Grant No. IORG005698).

Abbreviations
ace-1: Acetylcholinesterase-1

CREC : Centre de Recherche Entomologique de Cotonou

IRS: indoor residual spraying

ITNs: Insecticide-treated bed nets

LLINs: Long-lasting insecticidal nets

kdr: knockdown resistance

NMCP: National Malaria Control Program

PCR: Polymerase chain reaction

PMI: U.S. President’s Malaria Initiative s.l.: Sensu lato

USAID: United States Agency for International Development

WHO: World Health Organization

References

1. World Health Organization (2022) Global messaging briefing kit, World malaria report 2022. 21.
2. Nalinya S, Musoke D, Deane K (2022) Malaria prevention interventions beyond long-lasting insecticidal nets and indoor residual spraying in low- and middle-income countries: a scoping review. Malar J. 21(1): 31.
3. Akogbéto MC, Salako AS, Dagnon F, Aïkpon R, Kouletio M, et al. (2018) Blood feeding behaviour comparison and contribution of Anopheles coluzzii and Anopheles gambiae, two sibling species living in sympatry, to malaria transmission in Alibori and Donga region, northern Benin, West Africa. Malar J. 17(1): 307.
4. PMI VB (2021) The PMI Vectorlink Benin 2021 End of Spray Report (EOSR).
5. Prevention C-C for DC and (2019) CDC - Malaria - Malaria Worldwide - How Can Malaria Cases and Deaths Be Reduced? - Indoor Residual Spraying. https://www.cdc.gov/malaria/malaria_worldwide/reduction/irs.ht ml. Accessed 1 Feb 2024
6. RGPH-4 (2013) cahier des villages et quartiers de ville du département du Borgou. 31.
7. RGPH-4 (2013) cahier des villages et quartiers de ville du département de l’Alibori. 26.
8. RGPH-4 (2013) cahier des villages et quartiers de ville du département de la Donga. 24.
9. Odjo EM, Impoinvil D, Fassinou AJYH, Padonou GG, Aïkpon R, et al. (2024) The frequency of kdr and ace-1 alleles in Anopheles gambiae s.l. before and during indoor residual spraying (IRS) implementation and four years after IRS withdrawal in three districts in Atacora, Benin. Parasit Vectors. 17(1): 115.
10. Odjo EM, Tognidro M, Govoetchan R, Missihoun AA, Padonou GG, et al. (2024) Malaria transmission potential of Anopheles gambiae s.l. in indoor residual spraying areas with clothianidin 50 WG in northern Benin. Trop Med Health. 52(1): 18.
11. Gillies MT, Coetzee M (1987) A supplement to the anophelinae of Africa south of the Sahara (Afrotropical Region). South African Institute for Medical Research. Johannesburg.
12. Gillies MT, De Meillon B (1968) The Anophelinae of Africa south of the Sahara (Ethiopian zoogeographical region), 2d ed. South African Institute for Medical Research. Johannesburg.
13. Myriam et Cécile (2003) Protocoles de biologie moléculaire en usage au lin. Institut de recherche pour le développement IRD.
14. Scott JA, Brogdon WG, Collins FH (1993) Identification of Single Specimens of the Anopheles gambiae Complex by the Polymerase Chain Reaction. Am J Trop Med Hyg. 49(4): 520– 529.
15. Santolamazza F, Mancini E, Simard F, Qi Y, Tu Z, et al. (2008) Insertion polymorphisms of SINE200 retrotransposons within speciation islands of Anopheles gambiae molecular forms. Malar J. 7: 163.
16. Martinez-Torres D, Chandre F, Williamson MS, Darriet F, Bergé JB, et al. (1998) Molecular characterization of pyrethroid knockdown resistance (kdr) in the major malaria vector Anopheles gambiae s.s. Insect Mol Biol. 7(2): 179–84.
17. Weill M, Malcolm C, Chandre F, Mogensen K, Berthomieu A, et al. (2004) The unique mutation in ace-1 giving high insecticide resistance is easily detectable in mosquito vectors. Insect Mol Biol. 13(1): 1–7.
18. Salako AS, Ahogni I, Aïkpon R, Sidick A, Dagnon F, et al. (2018) Insecticide resistance status, frequency of L1014F Kdr and G119S Ace-1 mutations, and expression of detoxification enzymes in Anopheles gambiae (s.l.) in two regions of northern Benin in preparation for indoor residual spraying. Parasit Vectors. 11(1): 618.
19. Weir BS, Cockerham CC (1984) Estimating F-Statistics for the Analysis of Population Structure. Evolution. 38(6): 1358–1370.
20. Wright S (1984) Evolution and the Genetics of Populations, Volume 4: Variability Within and Among Natural Populations. University of Chicago Press, Chicago, IL.
21. Yadouleton AW, Asidi A, Djouaka RF, Braïma J, Agossou CD, et al. (2009) Development of vegetable farming: a cause of the emergence of insecticide resistance in populations of Anopheles gambiae in urban areas of Benin. Malar J. 8: 103.
22. Yadouleton AW, Padonou G, Asidi A, Moiroux N, Bio-Banganna S, et al. (2010) Insecticide resistance status in Anopheles gambiae in southern Benin. Malar J. 9: 83.
23. Luc DS, Benoit A, Laurette D, Michel M (2016) Indirect evidence that agricultural pesticides select for insecticide resistance in the malaria vector Anopheles gambiae. J Vector Ecol J Soc Vector Ecol. 41(1): 34–40.
24. Kouadio FA, Wipf NC, Nygble AS, Fodjo BK, Sadia CG, et al. (2023) Relationship between insecticide resistance profiles in Anopheles gambiae sensu lato and agricultural practices in Côte d’Ivoire. Parasit Vectors. 16(1): 270.
25. Moss S, Pretorius E, Ceesay S, Hutchins H, da Silva ET, et al (2024) Genomic surveillance of Anopheles mosquitoes on the Bijagós Archipelago using custom targeted amplicon sequencing identifies mutations associated with insecticide resistance. Parasit Vectors. 17(1): 10.
26. Grigoraki L, Cowlishaw R, Nolan T, Donnelly M, Lycett G, et al. (2021) CRISPR/Cas9 modified An. gambiae carrying kdr mutation L1014F functionally validate its contribution in insecticide resistance and combined effect with metabolic enzymes. PLoS Genet. 17(7): e1009556.
27. Syahrani L, Asih PBS, Bowolaksono A, Dwiranti A, Zubaidah S, et al. (2024) Impact of a spatial repellent intervention on Anopheles kdr insecticide resistance allele in Sumba, Indonesia. Malar J. 23(1): 31.
28. Aïzoun N, Aïkpon R, Akogbéto M (2014) Evidence of increasing L1014F kdr mutation frequency in Anopheles gambiae s.l. pyrethroid resistant following a nationwide distribution of LLINs by the Beninese National Malaria Control Programme. Asian Pac J Trop Biomed 4(3): 239–43.
29. Sovi A, Govoétchan R, Ossé R, Koukpo CZ, Salako AS, et al. (2020) Resistance status of Anopheles gambiae s.l. to insecticides following the 2011 mass distribution campaign of long-lasting insecticidal nets (LLINs) in the Plateau Department, south-eastern Benin. Malar J. 19(1): 26.
30. Clarkson CS, Miles A, Harding NJ, O'Reilly AO, Weetman D, et al. (2021) The genetic architecture of target-site resistance to pyrethroid insecticides in the African malaria vectors Anopheles gambiae and Anopheles coluzzii. Mol Ecol 30(21): 5303–5317.
31. Boussougou-Sambe ST, Ngossanga B, Doumba-Ndalembouly AG, Boussougou LN, Woldearegai TG, et al (2023) Anopheles gambiae s.s. resistance to pyrethroids and DDT in semi-urban and rural areas of the Moyen-Ogooué Province, Gabon. Malar J. 22(1): 382.
32. Gimonneau G, Pombi M, Choisy M, Morand S, Dabiré RK, (2012) Larval habitat segregation between the molecular forms of the mosquito Anopheles gambiae in a rice field area of Burkina Faso, West Africa. Med Vet Entomol. 26(1): 9–17.
33. Etang J, Mbida Mbida A, Ntonga Akono P, Binyang J, Eboumbou Moukoko CE, et al (2016) Anopheles coluzzii larval habitat and insecticide resistance in the island area of Manoka, Cameroon. BMC Infect Dis. 16: 217.
34. Yadouleton A, Martin T, Padonou G, Chandre F, Asidi A, et al (2011) Cotton pest management practices and the selection of pyrethroid resistance in Anopheles gambiae population in northern Benin. Parasit Vectors. 4 :60.
35. Fassinou AJYH, Koukpo CZ, Ossè RA, Agossa FR, Azondékon R, et al (2019) Pesticides and the evolution of the genetic structure of Anopheles coluzzii populations in some localities in Benin (West Africa). Malar J. 18(1): 407.
36. Djègbè I, Boussari O, Sidick A, Martin T, Ranson H, et al. (2011) Dynamics of insecticide resistance in malaria vectors in Benin: first evidence of the presence of L1014S kdr mutation in Anopheles gambiae from West Africa. Malar J. 10: 261.
37. Djogbénou L, Chandre F, Berthomieu A, Dabiré R, Koffi A, et al. (2008) Evidence of Introgression of the ace-1R Mutation and of the ace-1 Duplication in West African Anopheles gambiae s.s. PLOS ONE. 3(5): e2172.
38. Assogba BS, Milesi P, Djogbénou LS, Berthomieu A, Makoundou P, et al. (2016) The ace-1 Locus Is Amplified in All Resistant Anopheles gambiae Mosquitoes: Fitness Consequences of Homogeneous and Heterogeneous Duplications. PLOS Biol. 14(12): e2000618.
39. Grau-Bové X, Lucas E, Pipini D, Rippon E, van 't Hof AE, et al (2021) Resistance to pirimiphos-methyl in West African Anopheles is spreading via duplication and introgression of the Ace1 locus. PLOS Genet. 17(1): e1009253.
40. Claret JL, Di-Liegro M, Namias A, Assogba B, Makoundou P, et al. (2024) Despite structural identity, ace-1 heterogenous duplication resistance alleles are quite diverse in Anopheles mosquitoes. Heredity. 132(4): 179-191.
41. Kpanou CD, Sagbohan HW, Dagnon F, Padonou GG, Ossè R, et al (2021) Characterization of resistance profile (intensity and mechanisms) of Anopheles gambiae in three communes of northern Benin, West Africa. Malar J. 20(1): 328.
42. South A, Hastings IM (2018) Insecticide resistance evolution with mixtures and sequences: a model-based explanation. Malar J. 17(1): 80.
43. Gillespie JH, Turelli M (1989) Genotype-environment interactions and the maintenance of polygenic variation. Genetics. 121(1): 129–138.
44. Szulkin M, Bierne N, David P (2010) Heterozygosity-fitness correlations: a time for reappraisal. Evolution. 64(5): 1202–17.
45. Llaurens V, Whibley A, Joron M (2017) Genetic architecture and balancing selection: the life and death of differentiated variants. Mol Ecol. 26(9): 2430–2448.