Dynamics of pyrethroid resistance in malaria vectors in southern Benin following a large scale implementation of vector control interventions

Background Large-scale implementation of Indoor Residual Spraying and Insecticide Treated Nets has been implemented in Plateau Department, Benin between 2011 and 2014. The purpose of this study was to monitor the frequency and mechanisms of pyrethroid resistance in malaria vectors following the implementation of vector control tools for malaria prevention. Methods Anopheles larvae were collected in 13 villages twice a year from 2012 to 2014. WHO tube tests were used to assess the phenotypic resistance of each population to 0.05 % deltamethrin. Sibling species within Anopheles gambiae complex were identified by PCR techniques. Taqman and biochemical assays were performed to identify the presence of kdr mutations in individual mosquitoes and to detect any increase in the activity of enzymes putatively involved in insecticide metabolism (oxidases, esterase and glutathione-S-transferases). Quantitative real time PCR was used to measure the expression of three metabolic genes involved in pyrethroid resistance (CYP6P3, CYP6M2 and GSTD3). Results Anopheles populations showed < 90 % mortality to deltamethrin in all villages and at all time points. The 1014 F kdr allele frequency was close to fixation (> 0.9) over the sampling periods in both An. gambiae and An. coluzzii. Biochemical assays showed higher activities of alpha esterase and GST in field malaria vector populations compared to susceptible mosquitoes. qPCR assays showed a significant increase of CYP6P3, CYP6M2 GSTD3 expression in An. gambiae after a three-year implementation of LLINs. Conclusion The study confirmed that deltamethrin resistance is widespread in malaria vectors in Southern Benin. We suspect that the increase in deltamethrin resistance between 2012 and 2014 resulted from an increased expression of metabolic detoxification genes (CYP6M2 and CYP6P3) rather than from kdr mutations. It is urgent to evaluate further the impact of metabolic resistance on the efficacy of vector control interventions using pyrethroid insecticides.


Background
Global malaria vector control efforts rely on the use of Long Lasting Insecticide Nets (LLINs) and Indoor Residual Spraying (IRS). Twelve insecticides belonging to four chemical classes (organochlorines, organophosphates, carbamates and pyrethroids) are approved by the World Health Organization (WHO) for malaria vector control. All of these insecticides are neurotoxic and either target acetyl cholinesterase in the synapses or the voltage-gated sodium channel (VGSC). Pyrethroids are the only insecticides recommended by the WHO for LLINs because of their low mammalian toxicity, fast action, and high insecticidal activity [1]. Unfortunately, pyrethroid resistance has developed in most malaria vector species worldwide including Africa [2]. Indeed, two major mechanisms are known to confer pyrethroid resistance in malaria vectors: target site modification (kdr mutations) and increased metabolism of insecticides through detoxifying enzymes. The L1014F and L1014S substitutions in the para VGSC in the domains III-IV are known to decrease affinity of pyrethroids for this receptor [3]. Recently the mutation N1575Y has been described to potentiate the effect of the L1014F mutation [4]. The second resistance mechanism is called "metabolic" through higher catalytic properties and/or overexpression of carboxylesterases (COEs), cytochrome P450 mono-oxygenases (P450s) and Glutathione S-Transferases (GSTs) [5]. Some members of these families such as CYP6M2, CYP6Z2, CYP6P3 are known to contribute to pyrethroid detoxification in Anopheles mosquitoes [6]. Both metabolic and target site (kdr) resistance are present west Africa [7] and particularly Benin [8] and are suspected to reduce the efficacy of vector control intervention [9,10].
Since 2007 the National Malaria Control Program (NMCP) in Benin has implemented a nationwide distribution of LLINs. In 2011, 5,135,942 LLINs were distributed in Benin [11] followed by 5,663,220 in 2014 with a coverage of 97 % and 80 % respectively. In the study area, LLINs distributed in 2011 contained permethrin/ deltamethrin (pyrethroids) while those distributed in 2014 contained deltamethrin only (data from NMCP). After implementation of vector control interventions, it is essential to monitor any changes in susceptibility/resistance status of malaria vectors populations to public health insecticide (WHO 2015). Many cross-sectional studies were done to address the spatial distribution of pyrethroid resistance in Anopheles at a given time [12] but few longitudinal studies were conducted to address temporal changes in pyrethroid resistance phenotype following the implementation of vector control interventions.
The aim of this study was to investigate the dynamic of pyrethroid resistance in malaria vectors from 2012 to 2014 in 13 villages of Plateau Department and to characterize the mechanisms involved by monitoring changes in the frequency and expression of pyrethroidresistance markers. The outcomes shall help the NMCP to implement more effective vector control strategy against pyrethroid resistance populations.

Study area
The thirteen villages, selected on the basis of entomological and epidemiological criteria, were located in southeast Benin. They belonged to four districts: Ifangni, Sakete, Pobe, Ketou and were visited for mosquito collection ( Table 1). The study area is characterised by two rainy seasons from April to July and from September to October. The area is 3264 square km with a total population of 407,116 inhabitants (General Census of Population and Housing, 2002). Inhabitants of these villages are mostly farmers, traders, gardeners and fish breeders. Farmers grow cereals (maize, groundnuts, and beans), tubers (yams, manioc) and some vegetable crops such as tomato (Ifangni), chilli (Ketou). The fish breeding was only conducted in Itassoumba village (Ifangni) where tilapias and catfish were bred in large fishponds. These ponds provide a permanent breeding habitat for mosquitoes.

Mosquito collection
Mosquito collections were conducted during five rainy seasons. The larval sampling periods were June to August 2012 (termed June 2012), October to November 2012 (termed October 2012), May to July 2013 (termed June 2013), October to November 2013 (termed October 2013) and June to August 2014 (termed June 2014). In every survey, mosquito larvae were collected across the villages from several temporary breeding habitats, including household water storage [13]. Whenever possible, more than six larval habitats were examined per village. Mosquito larvae were transported to CREC, Cotonou and reared in an insectary under standard conditions (relative humidity 80 % ± 10 % and temperature 25°C ± 2°C). Adult mosquitoes were maintained with 10 % honey solution after emergence.

Insecticide susceptibility test
Bioassays were carried out on 2-5 old females using deltamethrin at the diagnostic dose of 0.05 % according to standard WHO procedures [14]. The laboratory susceptible reference strain of An. gambiae (Kisumu) was used to check the quality of the impregnated paper. After 1 h of insecticide exposure, mosquitoes were transferred to holding tubes and fed with 10 % honey solution. Mortality was recorded 24 h post-exposure. For each test, mosquitoes were also exposed to untreated paper to assess natural (control) mortality and to keep a batch of non-exposed mosquitoes for biochemical and molecular studies. Abbott's formula was used to correct the mortality when control mortality was between 5 and 20 %. After the test, legs were cut from control nonexposed mosquitoes for molecular determination and bodies were kept in RNA later at -20°C for mRNA expression. Another batch was frozen at -80°C for biochemical studies.

Molecular and biochemical assays
Genomic DNA was extracted using cetyl trimethyl ammonium bromide (CTAB) 2 % method modified from Doyle 1987 [15]. Briefly, mosquitoes were ground in CTAB 2 % then heated in a water bath at 65°for 5 min.
Chloroform was added to tubes, mixed by inversion, centrifuged and the upper phase transferred into another tube. DNA was precipitated with isopropanol and then washed once with 70 % cold ethanol. DNA was dried and suspended in distilled water. Species determination was performed by PCR [16]. The 1014 F and 1014S kdr mutations were detected by allelic discrimination Taqman assays as described by Bass [17] on field nonexposed females. Biochemical assay was used to quantify amounts of mixed function of oxidases (MFO), glutathione S-transferases (GST) and activities of non-specific esterase (NSE) using 30 female mosquitoes (non-exposed) for each village as described by Hemingway [18]. Each plate contained 10 unfed Kisumu adults used as the susceptible control. These biochemical tests were carried out on mosquitoes collected in June 2012 and June 2013 only due to insufficient sample size from others surveys.

RNA extraction and reverse transcription quantitative PCR
Pools of each Anopheles species were used to determine the relative gene expression of CYP6M2, CYP6P3 and GSTD3 by qPCR using SYBR Green. Total mRNA was extracted from batches of five mosquitoes (stored in RNA later) using Isolate RNA micro kit (Bioline) according to the manufacturer's instructions. Quantity and quality of mRNA were assessed using Nanodrop spectrophotometer (Nanodrop technologies). SuperScript III Reverse Transcriptase TM was used to synthesize first strand cDNA. The Kisumu was used as reference strain and the ribosomal gene RSP7 as housekeeping gene (shown to be consistent and with no differential expression between susceptible and resistant [19]). Three biological replicates were run for each sample and primers were designed on NCBI (http://www.ncbi.nlm.nih.gov/ tools/primer-blast/) ( Table 2). Real time PCR was run on Applied Bio systems ViiA7. Standard curves were generated using five times serially diluted cDNA sample to assess PCR efficiency. The PCR efficiency criterion was 100 ± 10 % for all of the genes and a single melting curve peak indicating specificity ( Table 2). The cDNA was diluted 10-fold in this concentration that fitted within the dynamic range of each qPCR and stored at -20°C.

Data analysis
Mann-Whitney test implemented in R 2.15.2 software is used to compare (i) mosquito mortalities between villages and sampling periods; and (ii) levels of enzymatic activity between the lab reference strain and field mosquitoes. A linear mixed model with villages as random effect implemented in R 2.15.2 software was used to test the effect of surveys and villages on mortality. WHO criteria for discriminating individuals for susceptibility/resistance status were applied: 98-100 % mortality indicating susceptibility; 90-97 % suspected resistance and < 90 % mortality-confirmed resistance [14]. Spearman's correlation test was used to investigate (i) the link between kdr allele frequency of field mosquitoes and mortality; (ii) between Anopheles species and mortality; and (iii) the fold expression of cytochrome P450 and the mortality by survey. All differences were considered significant for P-value < 0.05. The relative expression of target genes were determined according to ΔΔCt methods described by Schmittgen & Livak [20]. QPCR data were analysed using simple statistical randomization tests implemented in REST 2009 software.

Vector composition
Anopheles coluzzii and An. gambiae were found in sympatry in all villages. In June 2012, An. coluzzii was predominant in all sites (82 %) and An. gambiae and hybrids represented only 12 and 6 % of the species collected, respectively. Proportions of hybrids decreased over time to reach 0.3 % in June 2014. At the same time the proportion of An. coluzzii and An. gambiae fluctuated between 40 and 60 % (Fig. 1).  (Table 4). Kdr mutation was almost fixed in the study area and equally distributed between An. coluzzii and An. gambiae. Anopheles susceptibility to deltamethrin was not correlated with the kdr frequency over time (ρ = 0.146, CI 95% = 0.151-0.444, P = 0.336).

Metabolic resistance
Activities of esterases, glutathione S-transferases and mixed function of oxidases were measured using mosquitoes collected in June 2012 and June 2013. Results showed a higher GST activity in field Anopheles in June   (Table 3). Similarly, the level of α-esterase activity of field mosquitoes was significantly higher than that at Kisumu in June 2012. Conversely, we did not report higher activity of P450 in field populations at the two surveys compared to the susceptible reference strain (P > 0.05). The expression of detoxification genes was also measured by RT-qPCR in both species (Fig. 3). Only An. gambiae showed upregulation of metabolic genes compared to reference strain. GSTD3 was upregulated in June 2012 and June 2014, the mean fold changes (FC) were 1.8 and 49.3, respectively (P < 0.001) (Fig. 3). CYP6P3 and CYP6M2 were upregulated in June 2014 during the beginning of the rainy season; the FC were 4.8 and 44.42, respectively (P < 0.001). These genes were mostly down regulated in October 2012 and October 2013, suggesting a relationship between the level of expression of some P450 markers and the resistance phenotype to deltamethrin.

Discussion
This study addressed the dynamic of deltamethrin resistance in malaria vectors in the department of Plateau Benin, following a large scale implementation of malaria vector control tools. A combination of biological, biochemical and molecular assays were used to assess the frequency and mechanisms of pyrethroid-resistance in An. gambiae (s.l.) collected in 13 villages over three years. The results showed that An. coluzzii and An. gambiae were found in sympatry during all of the sampling periods but at various frequencies. A previous study demonstrated similar distribution of sibling species within An. gambiae (s.l.) complex in the study area [21].
Bioassay results showed significant variations of deltamethrin phenotypic resistance according to the surveys. Anopheles populations from the 13 villages were resistant to deltamethrin in June 2012 and June 2014 (< 90 % mortality) but mostly susceptible in October 2012 and 2013 (only 5 of 13 populations showed mortality < 90 %). This strong variation of phenotypic resistance in such a short period of time is difficult to explain knowing that environmental conditions did not change much between surveys (temperature 25 ± 2°C and humidity 80 ± 10 %). The frequency of the 1014 F kdr mutation was high in June 2012 in both species and did not increase much until June 2014. Consequently, it is unlikely that the kdr mutation alone Fig. 2 Phenotypic resistance to deltamethrin in Anopheles gambiae (s.l.) The scatter plot represents the mean mortalities with standard deviation based on 13 villages (a) and also according to surveys (b). The Mann-Whitney test showed significant time effect of sampling period on mosquito mortality while the variability between villages had no effect contributed greatly to the phenotype observed. A previous study questioned the causal association between kdr genotype and pyrethroid resistance (especially for type II pyrethroids) and suggested that the kdr genotype may not necessarily be the best predictor of resistance in malaria vectors. The 1014S kdr mutation originally from East Africa [22] was not detected in the present study. So far this mutation was found in Benin at a very low frequency in An. arabiensis [8].
In contrast, a significant increase in metabolic gene expression was reported in An. gambiae (s.s.) over the threeyear follow-up. Overexpression of three metabolic genes, CYP6M2 (> 40-fold), GSTD3 (> 40-fold) and CYP6P3 (> 4fold), were reported in June 2014, hence indicating a recent and strong selection pressure on this mosquito species.
Interestingly, the high gene expression was associated with low mortality rates as measured by WHO bioassays. This finding supports the involvement of these metabolic markers in deltamethrin-resistant phenotype. Downregulation of these genes in June 2012 despite the presence of resistance phenotype may be explained by the low number of An. gambiae collected in that survey (Table 4) and/or by the involvement of other (non-detectable) metabolic detoxification genes. The differential expression of metabolic markers between An. gambiae and An. coluzzi has been reported in Benin [23] and may reflect differential exposure to insecticides/xenobiotics at larval and/or adult stage.
CYP6M2 and CYP6P3 are regularly associated with pyrethroid resistance [24] and have been validated as     pyrethroid metabolizers [25,26]. No correlation was noted, however, between the oxidase activity and P450 gene expression. This can be explained by a lack of sensitivity and specificity of biochemical assays that employ generic heme peroxidase assays that are recognized by many members of the enzyme family [25]. The use of mixed populations (An. gambiae and An. coluzzii) during biochemical tests may have also affected the outcomes. Overexpression of CYP6M2 and CYP6P3 is widespread in West Africa including Benin, Nigeria [19], Ghana [26] and should be routinely monitored by national malaria control programmes. In addition, GSTD3 was found up regulated in Anopheles in June 2012 and 2014 hence correlating results of both, the biochemical test and bioassays. GSTD3 was found upregulated in DDT-resistant An. arabiensis and An. gambiae in Burkina Faso and Benin, respectively [8]. GSTs are regularly found overexpressed in many pyrethroid-resistant mosquitoes [27,28]. Some studies suggested their potential role against oxidative stress [29] and in pyrethroid sequestration [30]. Although the role of mosquito GSTs in pyrethroid resistance is likely, understanding the underlying mechanisms requires further investigations. In Benin, the use of pesticides for both vector control and agricultural practices (cotton crops and vegetable farms) are known to be a major source of selection pressure on malaria vectors [31][32][33]. For example, a randomized controlled trial conducted in southern Benin showed an increase in the kdr frequency from 20 to 80 % in An. gambiae eighteen months after the distribution of LLIN at community level [33]. We suspect that the increased coverage of deltamethrin-LLINs in 2014 in the study area (80 %, NMCP) has contributed to selection for pyrethroid-resistance metabolic markers. The role of agricultural practices in the selection of insecticide resistance in An. gambiae remains unknown. No significant changes in land use and agricultural practices were observed during the three-year follow-up but we could not record the type and amount of insecticides used by native farmers for crop protection. Similarly, the larval exposure to various xenobiotics (e.g. heavy metals, oils, fungicides, pollutants) that are known to modulate/enhance the metabolic detoxification profile of mosquitoes could not be investigated [34,35]. Clearly, much work has to be done to address the environmental factors contributing to resistance selection in malaria vectors in Benin.
The present study was part of a multidisciplinary project funded by the Bill & Melinda Gates Foundation that aims at addressing whether insecticide resistance can impact on the effectiveness of malaria vector control tools in Africa. Up to now the results of the first year follow-up did show substantial impact of insecticide resistance on the efficacy of LLINs in southern Benin [36] but analyses of years 2 and 3 are still ongoing.

Conclusions
This study monitored the levels and mechanisms of deltamethrin resistance in major malaria vectors in the Plateau department, Benin, after a three-year implementation of malaria vector control. The results showed that resistance to deltamethrin in malaria vectors was widespread and multifactorial. We also suspect that the increase in deltamethrin resistance between 2012 and 2014 resulted from an increased expression of metabolic detoxification genes (CYP6M2 and CYP6P3) rather than kdr mutations. It is now urgent to evaluate further the impact of metabolic resistance on the efficacy of vector control interventions using pyrethroid insecticides.