{'fraction_non_alphanumeric': 0.027367225887281483, 'fraction_numerical': 0.02264941369927373, 'mean_word_length': 5.280304893041554, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Experimental evolution supports the potential of neonicotinoid-pyrethroid combination for managing insecticide resistance in malaria vectors

The introduction of neonicotinoids for managing insecticide resistance in mosquitoes is of high interest as they interact with a biochemical target not previously used in public health. In this concern, Bayer developed a combination of the neonicotinoid clothianidin and the pyrethroid deltamethrin (brand name Fludora Fusion) as a new vector control tool. Although this combination proved to be efficient against pyrethroid-resistant mosquitoes, its ability to prevent the selection of pyrethroid and neonicotinoid resistance alleles was not investigated. In this context, the objective of this work was to study the dynamics and the molecular mechanisms of resistance of An. gambiae to the separated or combined components of this combination. A field-derived An. gambiae line carrying resistance alleles to multiple insecticides at low frequencies was used as a starting for 33 successive generations of controlled selection. Resistance levels to each insecticide and target site mutation frequencies were monitored throughout the selection process. Cross resistance to other public health insecticides were also investigated. RNA-seq was used to compare gene transcription variations and polymorphisms across all lines. This study confirmed the potential of this insecticide combination to impair the selection of resistance as compared to its two separated components. Deltamethrin selection led to the rapid enrichment of the kdr L1014F target-site mutation. Clothianidin selection led to the over-transcription of multiple cytochrome P450s including some showing high homology with those conferring neonicotinoid resistance in other insects. A strong selection signature associated with clothianidin selection was also observed on a P450 gene cluster previously associated with resistance. Within this cluster, the gene CYP6M1 showed the highest selection signature together with a transcription profile supporting a role in clothianidin resistance. Modelling the impact of point mutations selected by clothianidin on CYP6M1 protein structure showed that selection retained a protein variant with a modified active site potentially enhancing clothianidin metabolism. In the context of the recent deployment of neonicotinoids for mosquito control and their frequent usage in agriculture, the present study highlights the benefit of combining them with other insecticides for preventing the selection of resistance and sustaining vector control activities.

The introduction of neonicotinoids for managing insecticide resistance in mosquitoes is of high interest as they interact with a biochemical target not previously used in public health. In this concern, Bayer developed a combination of the neonicotinoid clothianidin and the pyrethroid deltamethrin (brand name Fludora Fusion) as a new vector control tool. Although this combination proved to be efficient against pyrethroid-resistant mosquitoes, its ability to prevent the selection of pyrethroid and neonicotinoid resistance alleles was not investigated. In this context, the objective of this work was to study the dynamics and the molecular mechanisms of resistance of An. gambiae to the separated or combined components of this combination. A field-derived An. gambiae line carrying resistance alleles to multiple insecticides at low frequencies was used as a starting for 33 successive generations of controlled selection. Resistance levels to each insecticide and target site mutation frequencies were monitored throughout the selection process. Cross resistance to other public health insecticides were also investigated. RNA-seq was used to compare gene transcription variations and polymorphisms across all lines. This study confirmed the potential of this insecticide combination to impair the selection of resistance as compared to its two separated components. Deltamethrin selection led to the rapid enrichment of the kdr L1014F target-site mutation. Clothianidin selection led to the over-transcription of multiple cytochrome P450s including some showing high homology with those conferring neonicotinoid resistance in other insects. A strong selection signature associated with clothianidin selection was also observed on a P450 gene cluster previously associated with resistance. Within this cluster, the gene CYP6M1 showed the highest selection signature together with a transcription profile supporting a role in clothianidin resistance. Modelling the impact of point mutations selected by clothianidin on CYP6M1 protein structure showed that selection retained a protein variant with a modified active site potentially enhancing clothianidin metabolism. In the context of the recent deployment of neonicotinoids for mosquito control and their frequent usage in agriculture, the present study highlights the benefit of combining them with other insecticides for preventing the selection of resistance and sustaining vector control activities. www.nature.com/scientificreports/ Malaria, caused by the protozoan Plasmodium parasite and transmitted by Anopheles mosquitoes is a major public health problem in Africa with 213 million cases and 405 000 deaths in 2018 1 . Without effective vaccines and the development of antimalarial drug resistance by parasites, the fight against malaria mainly relies on vector control operations using insecticides 2 through insecticide-treated nets (ITNs) and indoor residual spraying (IRS) 3,4 . Among the insecticides used in public health, pyrethroids remain the most used on ITNs because of their strong effects against the target insects and their low toxicity to the environment and mammals. In addition to pyrethroids, other insecticide families such as carbamates, organophosphates and in some specific cases DDT are also used in public health to improve the effectiveness of the fight against malaria vectors 5,6 . However, malaria vector control is now threatened by the spread of pyrethroid resistance throughout Africa 7 . Resistance to pyrethroids can be the consequence of various mechanisms in mosquitoes, such as non-synonymous mutations affecting the voltage-gated sodium channel targeted by these insecticides (i.e. Knock Down Resistance 'kdr' mutations), a lower insecticide penetration through insect cuticle, its sequestration, or its biodegradation (metabolic resistance) 8,9 . Kdr mutations can confer resistance to pyrethroids and DDT. These mutations prevent the insecticide from binding to the voltage-gated sodium channel (VGSC) of the insect's nervous system. Kdr mutations affecting malaria vectors are widely distributed in Africa with two distinct mutations at position 1014 of the VGSC mainly associated with insecticide resistance 10,11 : the L1014F kdr west mutation and the L1014S kdr east mutation 12,13 . Metabolic resistance is also widespread in African malaria vectors and has been associated with resistance to multiple insecticide families. Increased insecticide metabolism can be caused by an increased activity of detoxification enzymes or structural modification affecting insecticide turnover. Detoxification enzymes include cytochrome P450 monooxygenases (P450s or CYP for genes), carboxy/cholinesterases (CCEs), glutathione S-transferases (GSTs) and UDP-glycosyl-transferases (UDPGTs) although other families can be involved 9,14,15 . P450s from the CYP6 gene family and epsilon GSTs have been mainly associated with metabolic resistance to pyrethroids in An. gambiae 14 . Overall, though target-site modifications and metabolic resistance appear to play a central role in pyrethroid resistance, additional resistance mechanisms such as cuticle modification, altered transport and sequestration have also been reported [16][17][18] .
Most African An. gambiae populations are resistant to pyrethroid insecticides and also show varying levels of resistance to other insecticides used for vector control (carbamates, organophosphates and organochlorines). Insecticide resistance management therefore will benefit from the introduction of new chemistries. The fastest way to achieve that is by repurposing chemicals already used in agriculture, especially those with different modes of action which are less likely to be affected by cross-resistance with current public health insecticides 19 . Among the insecticides which are being proposed for managing pyrethroid resistance, neonicotinoids have aroused a high interest as (i) they have not previously been used for public health, (ii) they show a good toxicity against mosquitoes and a low toxicity to mammals 20 and (iii) they target nicotinic acetylcholine receptors (nAChR) which constitute a novel target for malaria vectors 21,22 . In this situation, a new insecticide formulation combining the neonicotinoid clothianidin and the pyrethroid deltamethrin (8:1 w/w) under brand name Fludora Fusion was developed by Bayer to be an additional tool for insecticide resistance management in African malaria vectors. By combining two unrelated modes of action, this indoor residual spraying (IRS) combination is intended to slow down the selection of resistance alleles to the newly introduced neonicotinoid clothianidin. However, the combination of clothianidin with deltamethrin, to which resistance alleles are already circulating in natural populations, may limit its resistance management potential. Field trials suggested that Fludora Fusion shows a good residual activity against malaria vectors from different African countries 23,24 and a good efficiency against pyrethroid-resistant An. gambiae with a higher efficacy than deltamethrin alone 25,26 . Although the mortality induced by Fludora Fusion was similar to clothianidin alone, these field trials also reported a higher knock down and repellency effect with Fludora Fusion as compared to clothianidin alone, suggesting a potential benefit of deltamethrin in the combination 25,26 . However, the hypothesis that such a combination of two modes of action offers a different selection profile across multiple generations, or indeed a benefit, as compared to the individual insecticide components has not been tested in malaria vectors. As neonicotinoids are often used against crop pests in Africa, resistance alleles to this insecticide family may already be circulating at low frequency in natural mosquito populations neighboring agriculture areas. Such hypothesis is supported by a recent study suggesting the presence of neonicotinoid resistance alleles in Anopheles coluzzi populations neighboring agriculture areas 27,28 . In such context, the large-scale implementation of Fludora Fusion for vector control might lead to the rapid selection of these resistance alleles leading to a decreased efficacy of this novel formulation, particularly if pyrethroid resistance alleles are also present.
In this context, the primary objective of this work was to compare the dynamics of resistance between Fludora Fusion and its two individual insecticide components in An. gambiae and to assess the ability of the combination to hinder resistance selection. The secondary objective of this work was to investigate the associated resistance mechanisms using molecular approaches. An An. gambiae line carrying resistance alleles to multiple insecticide families at low frequency was created by crossing a multi-resistant strain originating from an intense agricultural area of Côte d'Ivoire and a susceptible strain. The resulting line was then used for selecting, across multiple generations, three lines with Fludora Fusion, deltamethrin alone or clothianidin alone. The evolution of resistance levels to each insecticide and target site mutation frequencies were monitored throughout the selection process. The cross-resistance profiles of the selected lines to various insecticides were also investigated. A comparative RNA-seq approach was used to characterize the impact of selection on each line. Both differential gene transcription levels and selection signatures were investigated. The role of particular P450s in clothianidin resistance were then further explored using RT-qPCR and in silico protein modeling. The outcomes of the present study are discussed in the context of the deployment of neonicotinoids for malaria control and the added value of insecticide combinations for insecticide resistance management.

Results
Resistance dynamics during the selection process. In order to maximize the range of phenotypic and genotypic variations during the selection process, a moderately resistant line (Tiassalé-S) was created by masscrossing the mutli-resistant line Tiassalé with a fully susceptible line and then used for insecticide selection (see "Methods"). As expected from the introgression of susceptible alleles by controlled cross, bioassays confirmed the low resistance level of the Tiassalé-S line at G0 to most insecticide families with mortality levels higher than 95% for deltamethrin, bendiocarb and fenitrothion (see Supplementary Fig. S1 online). However, DDT resistance was still significant with 46% mortality observed to 4% DDT. Monitoring the resistance of each selected line to its respective insecticide along the selection process showed a rapid rise of deltamethrin resistance in the Delta-R line as compared to the non-selected line Tiassalé-S (Fig. 1). Mortality following deltamethrin exposure dropped from 73 to 30% at generation G7 and further decreased to stabilize around 15% (P < 0.05 from G7 to G33). In the Clothia-R line, a significant increased resistance to clothianidin was observed from G13 as compared to the non-selected line (P < 0.05 from G13 to G33) with resistance stabilizing at a moderate level. In the Fludo-R line, a significant increased resistance to Fludora Fusion mixture was observed at G13 though this was not confirmed in the subsequent generations. Comparing the deltamethrin resistance level of each line at generation G17 confirmed the high resistance of the Delta-R line but also suggested a slight increased resistance in the Fludo-R line (P < 0.05) ( Fig. 2A). A slight increased clothianidin resistance was observed in both Clothia-R and Fludo-R lines (P < 0.05). A slight increased resistance to Fludora Fusion mixture was observed in Delta-R and Clothia-R lines (P < 0.05) but not in the Fludo-R line (P > 0.05). Assessing the resistance level of each line to other insecticides used in public health revealed no significant increased resistance to fenitrothion or bendiocarb in any line (Fig. 2B). However, DDT resistance alleles appeared further enriched in all selected lines with mortality levels dropping to less than 10% (P < 0.05).
Target sites mutations. The genotyping of the kdr west mutation (L1014F) revealed its association with deltamethrin resistance in the Delta-R line with an increased frequency from 47% at G0 to nearly 100% at G7 with most individuals being homozygotes for the resistant allele (P < 0.05, Fig. 3). Although a transient rise of kdr mutation frequency was observed in both Clothia-R and Fludo-R lines at G2, this was followed by a gradual decrease in the following generations supported by a decreasing frequency of resistant homozygotes (P < 0.05 from G11 and G17 for Clothia-R and Fludo-R lines respectively). The frequency of the G119S Ace1 mutation were estimated at 10%, 7%, 14% and 10% in the Tiassalé-S line (G0), Delta-R line (G17), Clothia-R line (G17) and Fludo-R line (G17) respectively (Supplementary Table S2). As expected because not targeted by any insecticide used for selection, no significant frequency variation was observed in the selected lines as compared to the Tiassale-S line (P > 0.05 for all comparisons) though a slight excess of heterozygotes was observed in the Clothia-R line(P = 0.03).
Gene transcription levels. RNA sequencing performed at generation G17 produced > 120 M reads across four replicates for each line. Though ~ 10% of these reads were filtered out based on sequence and mapping quality, this allowed detecting transcription signal in adult females for 10,829 genes. Among them, 2596 genes were found differentially transcribed in at least one selected line as compared to the non-selected Tiassalé-S line (FC ≥ 1.5-fold and adjusted P value ≤ 0.005, see Supplementary Table S3 online). A strikingly higher number of genes were differentially transcribed in the Delta-R (1439 genes) and Clothia-R (1974 genes) lines as compared to the Fludo-R line (462 genes, Supplementary Fig. S4), supporting a lesser adaptive response in the Fludo-R line as compared to the two other selected lines. Only 31 and 17 genes were found over-and under-transcribed in all selected lines respectively. Gene Ontology (GO) term enrichment analysis performed on genes significantly over-transcribed in selected lines revealed a significant enrichment of terms associated with P450 activity in the Clothia-R line and in a lesser extent in the Fludo-R line (see Supplementary Fig. S5). These included the terms 'monooxygenase activity' , ' oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen' , 'heme binding' and 'iron ion binding' . Terms associated with acetylcholine receptors were enriched in the Clothia-R and Delta-R lines including the terms 'acetylcholine-activated cation-selective channel activity' and 'acetylcholine binding' . Enrichment analysis performed on genes significantly under-transcribed showed no enrichment of functions classically associated with resistance except a slight enrichment of terms associated with P450 and esterase activities in the Fludo-R line.
Focusing on gene families associated with known resistance mechanisms revealed different gene transcription profiles in each selected line as compared to the non-selected Tiassalé-S line (Fig. 4). Five P450s belonging to CYP325, CYP6 and CYP4 families, six esterases, two UDPGTs (AGAP006775 and AGAP009562) and four cuticular proteins were over-transcribed in the Delta-R line. Among them only CYP4J10 was previously shown to be over-transcribed in association with pyrethroid resistance. The Clothia-R line showed an over-transcription of multiple candidate genes belonging to families commonly associated with insecticide resistance. This included 16 P450s, most belonging to the CYP6, CYP9, CYP12 and CYP4 families. Some of them (CYP6M1, CYP12F1, CYP4C27, CYP6P1, CYP4G16, CYP6M3, CYP6Z3, CYP6Z2, CYP4H24, CYP6Z1, CYP6P3, CYP6P4) were previously associated with insecticide resistance in An. gambiae or showed a high protein similarity with P450s conferring neonicotinoid resistance in other insect species (see "Methods"). Five carboxylesterases were also over-transcribed in the Clothia-R line, with COEAE2G previously associated with insecticide resistance. Nine transferases were also over-transcribed in the Clothia-R line, most of them belonging to the UDPGT family. In addition to detoxification enzymes, clothianidin selection was also associated with the over-transcription of six ABC transporters and six cuticular proteins, among which CPR21 and CPAP3-A1b were previously associated with resistance. In addition, clothianidin and deltamethrin selection led to the over-transcription of multiple nicotinic acetylcholine receptor subunits. Fewer candidate genes were over-transcribed in the Fludo-R line with  Resistance levels were compared using bottle assays and are expressed as mean % mortality ± 95% Wald confidence interval. Comparisons between each line and the Tiassalé-S line at each time point were performed using a Fisher test on mortality proportions (*p < 0.05). 1 Indicate that a new clothianidin solution was used for monitoring resistance levels at G13, G17 and G33, leading to higher clothianidin mortality rates as compared to previous generations. . All lines were compared at generation G17. Resistance levels to the insecticides used for selection were compared using bottle assays while resistance levels to other insecticides were compared using WHO test tubes equipped with papers impregnated with 0.5% bendiocarb, 1% fenitrothion and 4% DDT. Mortality rates are expressed as mean % mortality ± 95% Wald confidence interval and were compared to the unselected line using a Fisher test on mortality proportions (*p < 0.05).

Figure 3.
Dynamics of the L1014F kdr mutation along the selection process. For each line, the L1014F kdr mutation was genotyped in individual mosquitoes using the allele-specific qPCR TaqMan assay described in Bass et al. 66 . For each selected line, both genotype and allele frequencies are reported. N number of individuals genotyped. P significance of genotype frequencies changes as compared to generation G2 (Chi 2 test, ns nonsignificant, *P < 0.05). www.nature.com/scientificreports/ 0.96%). A similar trend was also observed using the Bayesian F ST -based approach also slightly more outliers were identified (from 0.95% to 1.95%). Summing up outliers/differential SNPs at the gene level and projecting them on the AgamP4 genome revealed different selection signature profiles between the three selected lines ( Fig. 5 and Supplementary Table S6 online). Only a few loci showed strong selection signatures in the Delta-R line, none matching with gene previously associated with insecticide resistance. No selection signature was observed near the Kdr mutation locus because the low transcription level of the VGSC gene impaired the detection of polymorphisms. The closest polymorphic SNPs located ~ 80 Kb upstream and ~ 140 Kb downstream of the VGSC gene did not show any significant selection signature. Conversely, the Clothia-R line showed multiple loci showing strong selection signatures essentially located on chr 2R, 3R and 3L. The two loci located on chr 2R at ~ 5.5 Mb and on chr 3L at ~ 19.8 Mb were not associated with any candidate gene but the locus located on chr 3R at ~ 6.9 Mb matched the CYP6M locus previously associated with insecticide resistance. Finally, the Fludo-R line showed some selection signatures but none of them were overlapping those observed in the two other lines. The two strongest selection signatures observed on chr 3R at ~ 10 Mb and on chr 3R at ~ 11 Mb did not match any candidate gene or known resistance locus.
Association of CYP6Ms with clothianidin resistance. Among the 105 P450s present in An. gambiae genome, 30 belong to the CYP6 gene family among which 14 were located in the CYP6M locus showing a strong and specific selection signature in the Clothia-R line (Fig. 6A). In addition to CYP6M2 known to metabolize pyrethroids, homology analysis showed that this cluster contains 8 CYP6 genes (CYP6Y2, CYP6Y1, CYP6M1, CYP6M3, CYP6M4, CYP6Z1, CYP6Z2, CYP6Z3) showing high protein homology with P450s conferring neonicotinoid resistance in other insects (see "Methods"). Among them, CYP6Y1, CYP6M1 and CYP6M3 were all over-transcribed in the Clothia-R line. First, the association between CYP6M1, CYP6M2 and CYP6M3 transcription levels and clothianidin resistance was examined at generation G31 using RT-qPCR. The transcription level of each gene was compared between the Tiassalé-S and the Clothia-R lines in combination with pre-exposure to the P450 inhibitor piperonyl  Fig. S7). Although not significant (P = 0.07), CYP6M1 showed a slight over-transcription in the Clothia-R line as compared to the Tiassalé-S line while CYP6M2 and CYP6M3 were clearly under-transcribed. The over-transcription of CYP6M1 in the Clothia-R line further increased in survivors (P = 0.0014), suggesting that its over-expression is associated with clothianidin survival. Exposing Clothia-R individuals to 4% PBO alone did not significantly induce CYP6M1 transcription (P = 0.19). Pre-exposing Clothia-R individuals to PBO before clothianidin exposure significantly decreased survival by 20% (P < 0.05). Such exposure also affected CYP6M1 transcription level in survivors suggesting a lower survival of individuals over-transcribing CYP6M1 following its inhibition by PBO. Then, polymorphisms affecting CYP6M1 showing the strongest selection signature were further investigated. This gene contained 28 differential/outlier SNPs, among which 11 were non-synonymous (Fig. 6B) leading to the following amino acid changes: V35I, I66L, G171D, E215D, F240L, D347E, D351E, M375T, A377S, F396L, S401A. Among them, E215D, F240L and M375T were located within CYP6M1 substrate recognition sites SRS2, SRS3 and SRS5 respectively (Fig. 6C). In addition, E215D and F240L matched positions previously shown to interact with the binding of the neonicotinoid imidacloprid in B. tabaci CYP6CM1vQ, D. melanogaster CYP6G1, N. lugens CYP6AR1 and Ae. aegypti CYP6BB2. These point-mutations were present at low frequency in the Tiassalé-S line (~ 1.6%) and strongly enriched in the Clothia-R line (~ 97%) while their frequency remained lower in other selected lines. Comparing in silico models obtained for the Tiassalé-S and Clothia-R CYP6M1 proteins docked with clothianidin allowed identifying residues most likely interacting with the ligand (Fig. 6D). This analysis suggested a 17% increase in the volume of the binding pocket in the Clothia-R variant (869 Å 3 ) as compared to the Tiassalé-S variant (738 Å 3 ). The presence of a hydrogen bond between clothianidin and the Val-369 (distance of 3.6 Å) was predicted in both variants though another hydrogen bond with Leu-485 was predicted in the Clothia-R variant due the rotation of clothianidin in the active site. This rotation also led to a lower distance between the Asp-215 and clothianidin suggesting an interaction of this residue with the ligand in the Clothia-R variant. This rotation also led to a different positioning of the nitro-and methyl-groups of clothianidin whereas the nitroguanidine moiety was still facing the heme, suggesting that preferred metabolites could be desmethylclothianidin, desmethyl-denitro-clothianidin or clothianidin-urea 29,30 . Overall, although the amino acid changes identified in the SRS regions between Tiassalé-S and Clothia-R variants did not significantly differed in their physio-chemical properties, in silico models suggest that they may affect the positioning of clothianidin within CYP6M1 active site and subsequent enzyme-ligand interactions.

Discussion
The pyrethroid-neonicotinoid combination Fludora Fusion was pre-qualified for IRS usage by the World Health Organization in December 2018 31 . Field trials performed in multiple African countries confirmed the long lasting effect of this combination on various sprayed surfaces and its good efficacy against pyrethroid-resistant malaria vectors 25,26 . Though resistance risk assessment is not mandatory in the public health insecticides evaluation scheme, understanding the resistance dynamics of novel insecticide-based product is essential for anticipating resistance management strategies and extend their lifespan. In this context, the aim of the present study was to investigate through a controlled-selection experiment the potential of the African malaria vector An. gambiae to develop resistance this novel insecticide combination across multiple generations and investigate the associated mechanisms. Because the operational advantage of this novel product essentially stands on the association of two active ingredients with distinct modes of action, our controlled selection experiment compared resistance dynamics to each insecticide separately and to their combination. As neonicotinoid resistance alleles may already be circulating in mosquito populations indirectly exposed to neonicotinoids used in agriculture 28,28 , a field-derived An. gambiae line originating from the agricultural area of Tiassalé, Côte d'Ivoire was used as a parental line selection.

Fludora Fusion mixture and its components select different adaptive responses. Previous stud-
ies showed that An. gambiae from the Tiassalé area (Côte d'Ivoire) are highly resistant to the pyrethroid deltamethrin, the carbamate bendiocarb, the organochlorine DDT and moderately resistant to the organophosphate fenitrothion 32 . This multi-resistance phenotype was associated with a high frequency of the Kdr west L1014F mutation (~ 80%), a moderate frequency of the ace1 G119S mutation (~ 50%) and the over-transcription of multiple P450s previously associated with insecticide resistance 33 . Bioassays performed on the Tiassalé-S line confirmed its lower resistance to insecticides following the introgression of susceptible alleles. This was confirmed by the lower frequencies of the Kdr west L1014F and the Ace1 G119S mutations with most individuals carrying them being heterozygotes. As predicted by the presence of pyrethroid-resistance alleles in the Tiassalé-S line, deltamethrin selection rapidly led to an increased resistance to this insecticide in the Delta-R line. This resistance phenotype was tightly associated with the increased frequency of the Kdr west mutation, confirming the key role of this mutation in deltamethrin resistance 9,12 . Only a few detoxification enzymes were over-transcribed in the Delta-R line not including key pyrethroid metabolizers such as CYP6M2 and CYP6P3 previously identified in Tiassalé 15,33 or any other detoxification enzyme commonly associated with insecticide resistance 14,34 . Such absence of key metabolic enzymes in the Delta-R line suggests that the selection pressure applied (> 50% mortality) rather selected for individuals carrying the Kdr west mutation than those carrying metabolic resistance alleles. Hence, the rapid selection of this major mutation could have prevented the selection of metabolic resistance alleles of lower importance in deltamethrin resistance 35 . This was supported by polymorphism data showing no selection signature at metabolic resistance loci in the Delta-R line. Overall, this supports that the selection of target site resistance mechanisms are favored under intense selection pressure 36,37 . No such increased resistance was observed with Fludora Fusion mixture after 33 generations, suggesting that in our experimental conditions, the clothianidin-deltamethrin combination prevented the selection of resistance alleles having a significant impact on the resistance phenotype. The gradual decreased frequency of the Kdr west mutation in the Fludo-R line also supports the role of clothianidin in precluding its selection by the deltamethrin-clothianidin mixture.
However, we showed that selection with the neonicotinoid clothianidin alone led to an increased resistance to this insecticide although the resistance level was moderate and remained stable. RNA-seq data generated   Although such experiments should be interpreted with caution because of the limited number of replicates, the over-transcription of CYP6M1 in the Clothia-R line observed by RT-qPCR and its PBO-dependent association with clothianidin survival also supports its contribution to the resistance phenotype. In terms of polymorphisms, CYP6M1 and was affected by multiple non-synonymous mutations with two of them (E215D and F240L) located at positions previously suggested to interact with the neonicotinoid imidacloprid in B. tabaci CYP6CM1vQ, D. melanogaster CYP6G1, N. lugens CYP6AR1 and Ae. aegypti CYP6BB2 [38][39][40][41][42] . In silico modelling of CYP6M1 variants suggested that structural differences associated with these amino acid changes enhance clothianidin metabolism in the Clothia-R line. Although this needs to be functionally validated, the present study suggests that both transcriptional regulation and variant selection of CYP6M1 contribute to clothianidin resistance in An. gambiae. From a larger perspective, the present study supports the key role of P450s in neonicotinoids resistance as previously observed in other insect species 41,[43][44][45][46] .
In addition to detoxification enzymes, RNA-seq data pointed out the striking over-transcription of multiple nicotinic receptor subunits in the Delta-R and Clothia-R lines though no selection signatures were observed at these loci. Furthermore, no differential non-synonymous SNP was detected within these genes despite a good coverage, suggesting the absence of neonicotinoid target-site mutations in these lines. It has been shown that the modification of synaptic Na + homeostasis by 1014F Kdr mutation lead to a decrease of intracellular Ca 2+ concentration which in turn affects synaptic transmission 47 . As neonicotinoids also affect intracellular Ca 2+ concentration 48 , it is probable that the over-expression of nicotinic receptors observed in the Delta-R and Clothia-R lines rather reflect physiological compensation mechanisms to altered synaptic transmission than direct resistance mechanisms. Finally, a significant increase of DDT resistance was observed in all selected lines as compared to the Tiassalé-S line. This increased DDT resistance is likely associated with the high frequency of the kdr mutation in the Delta-R line, but this is less likely in the Clothia-R and Fludo-R lines showing lower kdr frequencies. However, genes encoding detoxification enzymes commonly associated with DDT resistance such as GSTE2 49 , CYP6Z1 50 and CYP6M2 14,17,33 were not over-transcribed in these two lines, suggesting that other resistance alleles might explain the observed phenotype. For instance, the over-transcription of CYP4G16 in these lines in conjunction with multiple cuticle proteins is of interest as this particular P450 was associated with cuticle-based resistance in An. gambiae 51,52 .

Pyrethroid-neonicotinoid combination as a new tool for vector control in Africa.
The use of the same insecticide classes targeting only two biochemical targets (the AchE and the VGSC) for decades has led to the selection and spread of insecticide resistance throughout Africa 53 . Insecticide resistance is now considered as a significant burden for malaria control with Anopheles populations often being resistant to multiple insecticide families 7 . Pyrethroid resistance is of major concern as these insecticides are mainly used for impregnating bednets 54,55 . In this context, the development of novel vector control products that can be used for controlling pyrethroid resistance has been encouraged by WHO 56 . Among them, IRS formulations combining two insecticides with different modes of action such as Fludora Fusion are of high interest. Though this IRS combination proved to be efficient for several months against African malaria vectors 23-25 , its long-term efficacy and its impact on insecticide resistance dynamics remains to be confirmed.
Although our controlled selection experiment was performed on a single field-derived An gambiae line and our results may have been impacted by drift effects, the absence of resistance development to FludoraFusion mixture suggest that clothianidin/deltamethrin combination represents an added value to current insecticides used for controlling malaria vectors in Africa. In particular, the potential of clothianidin to prevent the selection of Kdr mutations appears of high interest for managing pyrethroid resistance 57,58 . In turn, the selection of P450s associated with clothianidin resistance was impaired when deltamethrin was added, suggesting the added value of combining these two insecticides in delaying the emergence of neonicotinoid resistance. Such results are also in favor of the concomitant use of pyrethroid-impregnated bednets and neonicotinoid-based IRS products for vector control within the same location.
The present study also suggests that the use of clothianidin alone can rapidly select for metabolic resistance alleles in An. gambiae. P450-mediated cross-resistance between pyrethroids and neonicotinoids has previously been observed in various insect species [59][60][61] and our results support that this may also occur relatively rapidly in An. gambiae. Field studies are now needed to understand if these resistance alleles are only present in agricultural areas where neonicotinoids are regularly used or have already spread throughout Africa, and if they will be selected along the deployment of neonicotinoids for malaria control. In this frame, the candidate genes identified here provides the first set of neonicotinoid resistance markers in An. gambiae.

Conclusions
The present study supports the potential of neonicotinoid-pyrethroid combination as a novel vector control tool having the ability to delay the selection of resistance. As such combination appears to show a good efficacy against pyrethroid resistant populations, their deployment within an integrated vector control framework and under careful monitoring may be of added value for managing insecticide resistance in African malaria vectors. www.nature.com/scientificreports/ Conversely, the deployment of neonicotinoids alone rises concerns as this may rapidly lead to the emergence and spread of neonicotinoid resistance alleles in malaria vectors as observed in other insect species.

Methods
Use of animals. Blood feeding of adult mosquitoes was performed on mice. Mice were maintained in the animal house of the federative structure Environmental and Systems Biology (BEeSy) of Grenoble-Alpes University agreed by the French Ministry of animal welfare (agreement n° B 38 421 10 001) and used in accordance to European Union laws (directive 2010/63/UE). The use of animals for this study was approved by the ethic committee ComEth Grenoble-C2EA-12 mandated by the French Ministry of higher Education and Research (MENESR). The study was carried out in compliance with the ARRIVE guidelines.

Mosquitoes.
A. gambiae larvae were collected in 2010 in the vicinity of Tiassalé in the south of Côte d'Ivoire, an intensive agricultural area which has high coverage of pyrethroid-impregnated mosquito nets together with an intensive use of pesticides for crop protection including pyrethroids, organophosphates, carbamates and neonicotinoids 62 . Previous studies showed that malaria vectors from this area carry multiple resistance alleles conferring resistance to various insecticide families including pyrethroids, DDT, carbamates and organophosphates 32,33,63 . After species identification, An. gambiae sensus stricto specimens were used to create the Tiassalé strain that was then kept without insecticide selection pressure at the Centre Suisse de Recherche Scientifique en Côte d'Ivoire (CSRS). In order to maximize the dynamic range of resistant allele frequency variations during the selection process, a moderately resistant line (Tiassalé-S line) created by mass-crossing the Tiassalé line once with the fully susceptible line Kisumu was used for the controlled selection experiment. The Tiassalé-S line was then maintained without selection pressure for two generations before selection. All the lines described thereafter were maintained in the LECA tropical insectaries under standard conditions (29 °C, 90% relative humidity, 14 h/10 h light/dark period). Larvae were bred in deionized water and fed with TetraMin fish flakes. Adults were fed on filter papers impregnated with a 5% honey solution and blood feeding of adult females was performed on mice.
Controlled selection. The Tiassalé-S line was divided into four distinct lines. The first line (Tiassalé-S) was kept without selection and used as control for the selection experiment. The three other lines were selected in parallel with the pyrethroid deltamethrin (Delta-R line), the neonicotinoid clothianidin (Clothia-R line) or a mixture of deltamethrin and clothianidin in equivalent proportion (8:1 w/w) to Fludora Fusion (Fludo-R line). Because filter paper impregnation is not appropriate for clothianidin (crystallization of the active ingredient) all insecticide selections were performed using 250 mL glass bottles impregnated with insecticide. All insecticides used for bottle impregnation were obtained as pure active ingredients from Sigma. Deltamethrin was diluted in 100% acetone while clothianidin and clothianindin/deltamethrin mixture were diluted in an acetonic solution containing 17% (v/v) of mero solvent (MERO EC733, an emulsifiable concentrate containing 81.7% of rapeseed fatty acid esters in ethoxy-7-tridecanol). Bottles were impregnated with 1 mL of insecticide solution. For each line, insecticide selection consisted in introducing at least 15 batches of 20 non-blood fed females into insecticide-impregnated bottles. Doses used for selection were initially calibrated in order to reach 60% mortality in the in the Tiassalé-S line as follows: deltamethrin 10 µg/mL; clothianidin 0.45 µg/mL; clothianindin/deltamethrin mixture 0.45/0.05625 µg/mL (w/w ratio of 8:1, thereafter designated as Fludora Fusion mixture). Exposure times were as follows: deltamethrin 15 min, clothianidin 20 min, Fludora Fusion mixture 15 min. These times were slightly adjusted through the selection process in order to maintain an equivalent selection pressure between lines (between 50-70% mortality). Mortality rates were recorded after a recovery time of 72 h in order to consider the slower effect of clothianidin and Fludora Fusion mixture. Survivors were then transferred into new cages and blood-fed to generate eggs for the next generation. Selection was performed for each line at each generation until G33 except for generations G3, G4, G8, G13, G15, G17, G20, G22, G27, G29 and G32 in order to limit drift effects as population size was low (< 200 emerging females) for one line.

Resistance monitoring through the selection process. The resistance profile of the parental line
Tiassalé-S to the four main insecticide classes was characterized before the start of the selection experiment (G0) using Standard WHO susceptibility assays. Insecticide exposures were performed according to WHO standard procedures 64 using test tubes equipped with filter papers impregnated with the following insecticides: deltamethrin 0.05%; DDT 4%; bendiocarb 0.5% and fenitrothion 1%. Exposure time was fixed to 1 h for each insecticide.
At least five lots of 20 five-day-old adult females were tested for each insecticide. Mortality rate was recorded after a 24 h recovery time during which mosquitoes were provided a 5% honey solution. Mortality rates were expressed as mean mortality ± Wald's confidence intervals. The resistance level of each selected line to its respective insecticide was monitored at generations G0 (parental Tiassalé-S line), G2, G5, G7, G9, G11, G13, G17 and G33. Insecticide testing was performed using CDC bottles impregnated with a constant dose of insecticide corresponding to the initial dose used for selection (see above). Exposures were performed on 5 batches of 20 non-blood fed females per insecticide. Mortality rates were recorded after a 72 h recovery time in order to consider the slower effect of clothianidin. The resistance level of the non-selected Tiassalé-S line to each insecticide was also monitored in order to account for variations across generations inherent to insecticide solution preparation. Mortality rates were expressed as mean mortality ± Wald's confidence intervals. Differences between each selected line and the control line Tiassalé-S at each generation were tested using Fisher exact tests on mortality proportions (N = 5). At G17, the resistance level of each line to the three insecticides were compared using the same methodology in order to assess cross www.nature.com/scientificreports/ resistance patterns between lines. The resistance levels of each line to DDT, bendiocarb and fenitrothion were also measured at G17 using standard WHO susceptibility assays as described above.
Target sites mutations. The frequencies of target-site mutations previously identified in the Tiassalé strain (Kdr west L1014F and Ace1 G119S) were monitored for each selected line through individual genotyping. Kdr L1014F mutation was genotyped at generations G0, G2, G7, G9, G11 and G17. Because the acetylcholinesterase was not targeted by any insecticide used for selection, the frequency of the ace1 G119S was only genotyped at G0 and G17. At each generation, genomic DNA was extracted from 30 adult females of each line using the cetyltrimetylammonium bromide (CTAB) method 65 . Genomic DNA was resuspended in 20 µL nuclease-free water, quantified using the Qubit DNA BR assay (Thermofisher Scientific) and diluted to 0.5 ng/µL for genotyping. The Kdr L1014F and Ace1 G119S mutations were genotyped using the TaqMan qPCR methods described in 66 . Quantitative PCR reactions were performed on a CFX96 Real Time system (Bio-Rad) with PCR cycles as follows: 95 °C for 10 min, followed by 40 cycles of 95 °C for 10 s and 60 °C for 45 s. For each mutation, individuals were scored as homozygous susceptible/resistant or heterozygous based on the intensity of the HEX/FAM channels at the end of the PCR reaction as compared to positive and negative samples of known genotypes. Changes in genotype frequencies were tested at each generation using a Chi 2 test. For Kdr mutation, changes were tested as compared to generation G2 as generation G0 showed an abnormal excess of heterozygotes likely resulting from the initial controlled cross between the resistant line Tiassale (carrying Kdr mutation at high frequency) and the fully susceptible line Kisumu.
RNA library preparation and sequencing. Differential gene expression between each selected line and the non-selected line Tiassalé-S was investigated at the whole transcriptome level at G17 using RNA-sequencing (RNA-seq). For each line, four pools of 30 three-day-old non-blood fed females not previously exposed to insecticide were used.

Selection signatures associated with insecticide selection. Selection signatures associated with
insecticide selection were investigated using the 145,008 bi-allelic SNPs (substitutions or indels) that were polymorphic (i.e. showing ≥ 5% allele frequency variation between the Tiassalé-S line and any selected line). A first approach consisted in comparing the mean allele frequency of each SNP between each selected line and the Tiassalé-S line across the four replicates using a Student's test followed with a Benjamini and Hochberg multiple testing correction 68 . SNPs showing a mean frequency variation between any selected line and the Tiassalé-S line ≥ 50% in either direction and a corrected P value ≤ 0.001 were considered as associated with insecticide selection (Differential SNPs). A second approach consisted in assessing F ST departure from neutrality using the Bayesian method implemented in BayeScan version 2.1 73 . A separated analysis was performed for each selected line consisting in contrasting the selected line versus the Tiassalé-S line across all replicates. Default settings were used except that prior odd was set to 1000 in order to increase stringency. Genes showing a Bayscan Q-value of zero were considered as 'Outliers' . For each selected line, the proportions of 'Differential SNPs' and 'Outlier SNPs' per gene were computed and plotted along chromosomes using gene centers as genomic coordinates.
Impact of clothianidin selection on the CYP6M locus. The association of CYP6M1, CYP6M2 and CYP6M3 with clothianidin resistance were further studied at generation G31 using RT-qPCR. These genes were chosen based on their RNA-seq expression profile, selection signature and their potential role in insecticide resistance. For each gene the following conditions were compared: Tiassalé-S unexposed (Tiassalé-S); Clothia-R unexposed (Clothia-R); Clothia-R surviving clothianidin exposure (Clothia-R surv); Clothia-R exposed to the P450 inhibitor piperonyl butoxide PBO (Clothia-R PBO +); Clothia-R pre-exposed to PBO and surviving clothianidin exposure (Clothia-R PBO + surv). PBO exposure consisted in exposing batches of 25 mosquitoes to 4% PBO for 1 h in glass bottles as described in 28 . Clothianidin exposure consisted in exposing batches of 25 mosquitoes for 1 h to the same dose of clothianidin used for selection (0.45 µg/mL). Mosquitoes from all conditions were sampled at 72 h post clothianidin exposure (6-days old). Such recovery time allowed considering the slow effect of clothianidin together with minimizing gene induction/repression effects that may follow insecticide exposure. For each condition, total RNA was extracted from four pools of 12 to 20 mosquitoes using Trizol (Life Technologies) according to manufacturer's instructions. Two µg total RNA were treated with DNase I (Invitrogen) and reverse transcribed using superscript III (Invitrogen) and Oligo DT 20 primer according to manufacturer's instructions. cDNA samples were then diluted to 1/100 for qPCR amplification. Each sample was amplified as three PCR replicates. Quantitative PCR reactions were performed on a CFX qPCR system (Bio-Rad). Each qPCR reaction contained 12.5 µL iQ SYBR Green Supermix (Bio-Rad), 0.75 µL of each primer (10 µM each), 6 µL nuclease-free water and 5 µL diluted cDNA template. PCR cycles were as follows: 95 °C 3 min followed by 40 cycles consisting of 95 °C for 15 s and 30 s at 60 °C. Specific primer pairs were used for each gene (see Supplementary Table S8) and amplification specificity was verified by melt curve analysis. Data analysis was performed according to the ΔΔ Ct method taking into account PCR efficiency 74 using the gene encoding the ribosomal protein S7 (AGAP010592) as control. Results were expressed as mean transcription fold change (± SD) across the four biological replicates as compared to unexposed Tiassalé-S individuals. Differences across conditions were tested using an ANOVA across all conditions followed by Fischer tests between pairs of conditions (N = 4). Differential SNPs impacting CYP6M1 (AGAP008209) were further examined. Amino acid changes resulting from differential non-synonymous SNPs were mapped to CYP6M1 protein sequence and aligned using ClutalW against the protein sequence of four CYP6s known to metabolize neonicotinoids: Drosophila melanogaster CYP6G1 39-41 , Bemiscia tabaci CYP6CM1vQ 42 , Nilaparvata lugens CYP6ER1 75 and Aedes aegypti CYP6BB2 38 . P450 substrate recognition sites (SRS) were identified from protein alignment as defined from CYP3A4 structure 76 . The impact of amino acid changes selected by clothianidin on CYP6M1 structure and clothiainidin docking was then investigated by comparing the protein models obtained for the two distinct CYP6M1 variants identified in the Tiassalé-S and Clothia-R lines. For each variant, the amino acid sequence used for protein modelling differed by the 11 amino acid changes showing a differential frequency higher than 50% (see Supplementary Fig. S9). CYP6M1 protein models were computed by the SWISS-MODEL server using the crystallized human P450 CYP3A4 as template (protein Accession 5VCD). The Rosetta algorithm was then used to refine protein structure and to compute the most probable positions for the heme and clothianidin in the active site of each variant according to Rosetta-computed force fields 77 . The volume of the binding pockets of each variant was estimated using POVME 2.0 78 .