Genome-Wide Transcriptional Analysis and Functional Validation Linked a Cluster of Epsilon Glutathione S-Transferases with Insecticide Resistance in the Major Malaria Vector Anopheles funestus across Africa

Resistance is threatening the effectiveness of insecticide-based interventions in use for malaria control. Pinpointing genes associated with resistance is crucial for evidence-based resistance management targeting the major malaria vectors. Here, a combination of RNA-seq based genome-wide transcriptional analysis and RNA-silencing in vivo functional validation were used to identify key insecticide resistance genes associated with DDT and DDT/permethrin cross-resistance across Africa. A cluster of glutathione-S-transferase from epsilon group were found to be overexpressed in resistant populations of Anopheles funestus across Africa including GSTe1 [Cameroon (fold change, FC: 2.54), Ghana (4.20), Malawi (2.51)], GSTe2 [Cameroon (4.47), Ghana (7.52), Malawi (2.13)], GSTe3 [Cameroon (2.49), Uganda (2.60)], GSTe4 in Ghana (3.47), GSTe5 [Ghana (2.94), Malawi (2.26)], GSTe6 [Cameroun (3.0), Ghana (3.11), Malawi (3.07), Uganda (3.78)] and GSTe7 (2.39) in Ghana. Validation of GSTe genes expression profiles by qPCR confirmed that the genes are differentially expressed across Africa with a greater overexpression in DDT-resistant mosquitoes. RNAi-based knock-down analyses supported that five GSTe genes are playing a major role in resistance to pyrethroids (permethrin and deltamethrin) and DDT in An. funestus, with a significant recovery of susceptibility observed when GSTe2, 3, 4, 5 and GSTe6 were silenced. These findings established that GSTe3, 4, 5 and 6 contribute to DDT resistance and should be further characterized to identify their specific genetic variants, to help design DNA-based diagnostic assays, as previously done for the 119F-GSTe2 mutation. This study highlights the role of GSTes in the development of resistance to insecticides in malaria vectors and calls for actions to mitigate this resistance.


Introduction
Malaria is the deadliest vector-borne disease, killing more than 400,000 people every year [1]. Vector control interventions through the use of long-lasting insecticide nets and the implementation of indoor residual spray have led to a significant reduction in malaria incidence, between 2000 and 2015 [2]. This gain is under threat, as the most recent WHO World Malaria Report revealed there has been increase in annual case numbers since 2016. This malaria rebound is partly due to escalation of insecticide resistance in major malaria vectors such as An. gambiae and An. funestus. This was recently shown for a population of An. funestus for which high resistance level to pyrethroids was associated with a significant loss of efficacy of insecticide-treated nets including PBO-based nets [3]. The main mechanisms of resistance are target-site and metabolic resistance. The molecular basis of metabolic resistance is more complex and involves among others an overexpression and/or over-activity of major detoxification genes, such as cytochrome P450s (CYP450s), glutathione-S-transferases (GSTs) and the carboxylesterases [4][5][6], in addition to the recently described sensory appendages proteins [7]. Several studies, including genome-wide transcriptional analyses using microarray/qPCR and functional validation have linked GST with resistance in the major malaria vectors [8][9][10]. However, most of these studies have concentrated on GSTe2, neglecting the other GST epsilon genes, even though they have been shown to also consistently be overexpressed [11][12][13][14][15]. Indeed, among the 8 potential members of the GST epsilon class, several genes have previously been shown to be overexpressed in malaria vectors including An. gambiae [10] and An. funestus [14][15][16]. However, apart from few validations such as for GSTe4 in An. gambiae and An. arabiensis [17] little is known on the role of these genes in insecticide resistance. In the case of An. funestus, progress made has focused on the GSTe2 gene detecting a key resistance marker (L119F-GSTe2) now commonly used for resistance monitoring in An. funestus [16]. In addition to conferring pyrethroid/DDT resistance, it was also shown that the L119F-GSTe2-mediated metabolic resistance to pyrethroids/DDT is associated with negative effects on some life-history traits of field populations of An. funestus, supporting that insecticide resistance is associated with a fitness cost [18]. An experimental hut study, using the same marker in An. funestus population from Mibellon (Cameroon) had confirmed that presence of the L119F-GSTe2 was associated with resistance to DDT and pyrethroids [19] and was reducing the efficacy of bed nets [20].
RNA interference (RNAi) is one of the main approaches commonly used for in vivo validation of the role of detoxification enzymes in conferring resistance to insecticides in mosquitoes. This is done by injecting adult mosquitoes with double-stranded RNA (dsRNA) corresponding to the gene of interest. In turn, this induces mRNA degradation via the RNAi pathway and suppresses expression of the protein [21,22]. RNAi has been used to link overexpression of cytochrome P450s and GSTs to insecticide resistance in An. gambiae and Aedes mosquitoes [23][24][25]. Recently, RNAi has been used to establish the role of sensory appendage proteins in the leg of An. gambiae in conferring pyrethroid resistance [7] confirming that this method provides a robust approach to validate the contribution of specific genes such as those of the GST epsilon, to a specific phenotype.
In this study, using a genome-wide RNA-seq-based transcriptomic analysis, we detected candidate genes associated with DDT-resistance and DDT/permethrin crossresistance across different regions of Africa revealing the pre-eminence of GST epsilon genes. The role of epsilon class of GSTs in DDT/pyrethroid resistance was investigated in An. funestus population from Mibellon (Cameroon) using RNAi-mediated gene silencing. The results suggest that several members of this class in addition to GSTe2 contribute to the overall resistance observed in the field.

Mosquito Collection and Rearing
Indoor-resting female An. funestus were collected early in the morning (6:00 a.m.-8:00 a.m.), using battery-powered aspirators (John. W. Hock, Gainesville, FL, USA) in . The collection was done from randomly selected houses, following a verbal consent from the chief of the district and the household owners. Mosquitoes collected were kept in paper cups, being transported to the insectary in the Centre for Research in Infectious Diseases (CRID), Yaoundé, Cameroon. The F 1 generation was generated in the insectary from field blood-fed female mosquitoes using forced-egg laying method [26]. Briefly, F 0 gravid females were transferred into 1.5 mL Eppendorf tubes containing a wet filter paper, to enable them to lay eggs. After oviposition the parents that laid eggs were removed and used for species identification. Molecular identification to species level was carried out according to the protocol described priviously [27] using genomic DNA (gDNA) extracted from F 0 females identified morphologically as An. funestus. DNA extraction was done using the Livak protocol [28]. The eggs were pooled into bowls and supplemented with TetraMin™ baby fish food. The emerged F 1 female progenies were mixed in cages and 2 to 5-day old females used for insecticide bioassays and double-stranded RNA (dsRNA) injection.

Transcriptomic Profiling of DDT Resistance across Africa Using RNA-Seq
Transcriptional profiling of An. funestus populations was carried out to detect key candidate genes associated with DDT resistance across Africa. This was done using mosquitoes from four African regions: Central (Mibellon/Cameroon), West (Obuasi/Ghana), Southern (Chikwawa/Malawi), and East (Tororo/Uganda). Total RNA was extracted from pools of 10 female DDT-resistant mosquitoes (alive after 24 h exposure to DDT), unexposed mosquitoes (control) and lab susceptible colony (FANG) using the Arcturus PicoPure RNA Isolation Kit (Life Technologies, Carlsbad, CA, USA) according to the manufacturer's instructions. The FANG is a fully insecticide susceptible laboratory colony [29]. In addition, RNA was extracted from 10 permethrin-resistant mosquitoes (alive after 24 h exposure to permethrin) to study the gene differentially expressed across Africa when comparing DDT vs. permethrin-resistant mosquitoes.
RNA libraries were pooled in equimolar amounts using the Qubit and Bioanalyzer data. The quantity and quality of each pool were assessed by Bioanalyzer and subsequently by qPCR using the Kapa Illumina library quantification kit (Kapa Biosystems, Wilmington, MA, USA), on a Light Cycler LC480II (Roche, Basel, Switzerland), according to manufacturers' instructions. The pool of libraries was sequenced on one lane of the HiSeq 2500 (Illumina, San Diego, CA, USA) at 2 × 125 bp paired-end sequencing with v4 chemistry. Sequence library preparation, sequencing, initial processing and quality control were done by the Centre for Genomic Research, University of Liverpool, UK. Alignment to the reference sequence using the AfunF3.1 annotation. (https://vectorbase.org/vectorbase/app/record/dataset/DS_1a787d4361#pmids, accessed on 7 April 2020). Data were analysed as described previously [30][31][32]. Differential gene expression analysis was performed using edgeR and the Strand NGS program (Strand Life Sciences, version 3.0, Hebbal, Bangalore, India).

Investigation of Expression Profile of GSTe Genes in An. funestus across Africa
The expression profiles of GSTe1, GSTe2, GSTe3, GSTe4, GSTe5, GSTe6, GSTe7 and GSTe8 in DDT and permethrin-resistant mosquitoes was assessed across Africa (Cameroon, Benin, Malawi, and Uganda), using qRT-PCR. RNA was extracted from three biological replicates from each population of 10 each of DDT-resistant females (alive 24 h after exposure to DDT), permethrin-resistant (alive 24 h after exposure to permethrin), control (An. funestus mosquitoes not exposed to any insecticide), as well as susceptible laboratory colony (FANG). Briefly, 1 µg of the total RNA from each of the three biological replicates was used as the template for cDNA synthesis using Superscript III (Invitrogen, Carlsbad, CA, USA) with oligo-dT20 and RNase H, according to the manufacturer's instructions. The qRT-PCR amplification was performed as described [15] using the primers provided in Supplementary Table S2.1. The relative expression and fold change of each GSTe gene was calculated as previously described [33] by comparing expression in resistant, susceptible and control samples. The normalization was done with the ribosomal protein S7, RPS7 (AFUN007153) and actin5C (AFUN006819) housekeeping genes.  Table S2.2). Specific GST2, 3, 4, 5, 6, 7 and GSTe8 fragments were amplified by PCR from plasmid clones using KAPA Taq Kit (Kapa Biosystems, Wilmington, MA USA). Double-stranded RNA (dsRNA) was synthesized using in vitro transcription MEGAscript ®® T7 Kit (Ambion Inc., Austin, TX, USA) and purified using MEGAclear columns (Ambion). The purified products were concentrated by ethanol precipitation and the dsRNA was resuspended in nuclease-free water and stored at −20 • C. The successful construction of dsRNA was confirmed by running 3 µL of dsRNA-diluted products in 1.5% agarose gel in a Tris-acetate-EDTA (TAE) buffer.

Mosquitoes Injection and Susceptibility Bioassays
To explore the role of GSTe genes in conferring insecticide resistance, RNAi was performed on Mibellon An. funestus population, by injecting sequence-specific dsRNA to 2-3 days old F 1 female mosquitoes, followed by insecticide bioassay. A Nano injector (Nanoinject; Drummond, Burton, OH, USA) was used to inject dsGSTe2, 3, 4, 5, 6, 7 and dsGSTe8 into the thorax of 2 to 3 days old female An. funestus mosquitoes as described [34]. Briefly, mosquitoes, induced to sleep with CO 2 , were injected with 69 nL of either aliquot of above dsGSTes or dsGFP (control). Four days after injection, four replicates of 20 mosquitoes for each dsRNA were exposed to permethrin (0.75%), deltamethrin (0.05%) and DDT (4%) for 1 h following the WHO testing protocol [35]. Mosquitoes were transferred to holding tubes after exposure, supplemented with sugar and mortalities counted 24 h after the exposure. The susceptibility test was performed in triplicate with experimental mosquitoes comprising the mosquitoes injected with dsGSTes above, whereas mosquitoes injected with dsGFP and those not injected were used as controls.

Quantitative RT-PCR to Confirm the Knockdown Effect
For dsGSTe-injected and non-injected mosquitoes, RNA was extracted from 3 pools of 5 mosquitoes using TRIzol reagent (Gibco BRL, Gaithersburg, MD, USA). cDNA from each of the three biological replicates was synthesized using the Super-Script III (Invitrogen, Carlsbad, CA, USA) with oligo-dT20 and RNase H, according to the manufacturer's instructions. The cDNA from each replicate treatment was then used to assess the extent of RNAi by measuring levels of gene expression after injection by qRT-PCR. To assess the knockdown efficiency after injection and quantitative difference in the level of GSTes expression between injected and non-injected mosquitoes, a standard curve of each gene was established using a serial dilution of cDNA. The qPCR amplification was carried out in a MX3005 real-time PCR system using Brilliant III Ultra-Fast SYBR Green qPCR Master Mix (Agilent, Santa Clara, CA, USA). A total of 10 ng of cDNA from each sample was used as a template in a three-step program involving a denaturation at 95 • C for 3 min followed by 40 cycles of 10 s at 95 • C and 10 s at 60 • C and a last step of 1 min at 95 • C, 30 s at 55 • C and 30 s at 95 • C. The relative expression and fold-change of each target gene were calculated according to the 2 −∆∆CT Livak method [33], comparing expression in specific dsGSTe-injected samples to non-injected ones, after normalization with the housekeeping genes, RPS7 (AFUN007153) and actin5C (AFUN006819), as described above.

Data Analysis
All analyses were conducted using GraphPad Prism version 7.00, R 3.3.2. for Windows and Strand NGS program (Strand Life Sciences, version 3.0, Hebbal, Bangalore, India). Students' t-test was used to assess statistical differences between experimental and control groups.

RNAseq-Based Comparative Transcriptomic Profiling of DDT Resistance across Africa
To detect genes associated with DDT resistance in An. funestus mosquitoes Africa-wide, transcriptional profiling of mosquitoes from different regions of Africa was performed. This comprised populations from southern (Malawi), East (Uganda), West (Ghana) and Central (Cameroon) Africa. Priority was given to the comparison between genes upregulated in DDT-resistant mosquitoes (R) and the control (C, unexposed mosquitoes) because this comparison directly focuses on the difference between mosquitoes having the same genetic background (accounting for potential induction of expression), but differing in treatment received. Attention was also given to genes that were commonly upregulated in R vs Susceptible FANG colony mosquitoes (S) and C vs S as these genes are the ones expressed constitutively in natural mosquitoes' populations. The number of differentially expressed genes in each of the four populations and the FANG susceptible strain is shown in Venn diagrams ( Figure 1). Raw data from RNA-seq is deposited on sequence archive, with the following link: https://www.ebi.ac.uk/ena/browser/view/PRJEB24351, accessed on 10 January 2018.

Data Analysis
All analyses were conducted using GraphPad Prism version 7.00, R 3.3.2. for Windows and Strand NGS program (Strand Life Sciences, version 3.0, Hebbal, Bangalore, India). Students' t-test was used to assess statistical differences between experimental and control groups.

RNAseq-Based Comparative Transcriptomic Profiling of DDT Resistance across Africa
To detect genes associated with DDT resistance in An. funestus mosquitoes Africawide, transcriptional profiling of mosquitoes from different regions of Africa was performed. This comprised populations from southern (Malawi), East (Uganda), West (Ghana) and Central (Cameroon) Africa. Priority was given to the comparison between genes upregulated in DDT-resistant mosquitoes (R) and the control (C, unexposed mosquitoes) because this comparison directly focuses on the difference between mosquitoes having the same genetic background (accounting for potential induction of expression), but differing in treatment received. Attention was also given to genes that were commonly upregulated in R vs Susceptible FANG colony mosquitoes (S) and C vs S as these genes are the ones expressed constitutively in natural mosquitoes' populations. The number of differentially expressed genes in each of the four populations and the FANG susceptible strain is shown in Venn diagrams ( Figure 1). Raw data from RNA-seq is deposited on sequence archive, with the following link: https://www.ebi.ac.uk/ena/browser/view/PRJEB24351, accessed on 10 January 2018.

The Central Africa Population of Cameroon
The major detoxification gene families found to be overexpressed in Cameroon population are cytochrome P450s, glutathione S-transferases and carboxylesterases. When comparing resistant, susceptible and control mosquito cytochrome P450 CYP325A is over-expressed, followed by CYP6P9b, CYP6P5, CYP315A1. For GSTs, the epsilon and delta family are upregulated when comparing R-S vs. C-S expression profile. The most overexpressed GSTs are GSTe2, GSTe1, GSTe3, GSTe6, GSTd3 and GSTt2. Several genes from carboxylesterase classes, e.g., an unknown COE (AFUN002514), and COEBE3C (AFUN016311, a glutathione peroxidase (AFUN022201), an ATP-binding cassette transporter (AFUN019220), a UDP-glucuronosyltransferase (AFUN011266), and sulfotransferase family (AFUN016207) (SULT1B), were also found to be overexpressed (Table 1). Table 1. Detoxification genes differentially expressed in Cameroon An. funestus between different comparisons at false discovery rate (FDR) < 0.05 and fold change (FC) > 1.5 for genes induced upon exposure and constitutive expression (R-C) or FC > 2 for constitutive differential expression (C-S) and genes induced upon exposure and constitutive expression (R-S).

The Southern Africa Population of Malawi
The most overexpressed genes in the Malawi population when comparing R to S are the two P450s CYP6P9a and b (Table 3). Besides those two genes, many other P450s, including CYP325J1, CYP6M4, CYP6P2, CYP9K1, CYP314A1, CYP6N1 and cytochrome b5, are upregulated. For GSTs, the most overexpressed genes are from delta family, e.g., GDTD1 and GSTD11. The GST epsilon family include GSTe1, GSTe2, GSTe5 and GSTe6. In addition, theta family GSTt1 was also found to be overexpressed in exposed mosquitoes. As for Central and West Africa, we also have the overexpression of carboxylesterase (AFUN016265), sulfotransferase (AFUN016207) and ATP-binding cassette transporter (AFUN019220). Table 3. Detoxification genes differentially expressed in Malawi between different comparisons at FDR < 0.05 and FC > 1.5 for R-C or FC > 2 for C-S and R-S.

The East African Population of Uganda
Here, the most overexpressed genes are cytochrome P450s, with CYP4C26, CYP6P5, CYP6P4a, CYP306A1, CYP305A3 and CYP315A1. Some GST families are also significantly over-expressed when comparing R to C mosquitoes. These include GSTD1, GSTD3 and GSTe6. However, not many GST epsilon genes are overexpressed in Uganda compared to other Africa countries, e.g., Cameroon and Ghana. Many other genes families, such as carboxyesterases, sulfotransferase, NADHP, and ATP-binding cassette transporter, are expressed in Uganda (Table 4). Table 4. Detoxification genes differentially expressed in Uganda between different comparisons at FDR < 0.05 and FC > 1.5 for R-C or FC > 2 for C-S and R-S. Noticeably, the analysis of gene expression across Africa between DDT-resistant mosquito population revealed that GSTs are upregulated across the continent. The level of expression of GSTs is variable across the continent and three major families, the epsilon, the delta and the theta family, are the most overexpressed. Regarding the GST epsilon cluster specifically, we observed that all genes are up-regulated across Africa when comparing control mosquitoes versus DDT-resistant population, except for GSTe8.  Table S1.1).

qPCR Transcriptional Profiling of GSTe Genes in An. funestus across Africa
The validation of the expression profile of genes of the GSTe cluster across the continent was carried out using qPCR, comparing the expression level of GSTes between FANG, unexposed (control), permethrin-alive and DDT-alive mosquitoes after 24 h exposure. In Benin, the level of expression of GSTe2, GSTe3 and GSTe4 was significantly higher in mosquitoes resistant to DDT compared to the unexposed group and the permethrin-resistant mosquitoes ( Figure 2). In Uganda, all the GSTes were expressed but at comparatively the same level, with the most over-expressed one being GSTe6, followed by GSTe8. As observed with Benin samples, these two genes were more expressed in mosquitoes surviving DDT exposure. In Malawi, it was observed in ascending order that the most expressed GSTes were GSTe2, GSTe3, GSTe6, GSTe5 and GSTe8. However, except for the GSTe6, which was more overexpressed in mosquitoes resistant to permethrin, the level of expression of other GSTes was higher in DDT-resistant mosquitoes. For Cameroon, GSTe2, GSTe4 and GSTe3 were more overexpressed in DDT-alive females. Overall, these results showed that the GSTes are overexpressed mainly in DDT-resistant mosquitoes compared to permethrin and unexposed mosquitoes.
surviving DDT exposure. In Malawi, it was observed in ascending order that the most expressed GSTes were GSTe2, GSTe3, GSTe6, GSTe5 and GSTe8. However, except for the GSTe6, which was more overexpressed in mosquitoes resistant to permethrin, the level of expression of other GSTes was higher in DDT-resistant mosquitoes. For Cameroon, GSTe2, GSTe4 and GSTe3 were more overexpressed in DDT-alive females. Overall, these results showed that the GSTes are overexpressed mainly in DDT-resistant mosquitoes compared to permethrin and unexposed mosquitoes.

Confirmation of GSTe Knockdown Effect by qRT-PCR
To confirm whether the injection of dsGSTes did knock-down the expression of GSTe genes in vivo in mosquitoes, qPCR was performed using the cDNA from unexposed dsGSTes (injected) and non-injected mosquitoes with the primers of each GSTes cluster, using actin5C and RPS7 as housekeeping genes. As shown in Figure 3, we noticed a significant partial reduction in GSTes gene expression when comparing control non-exposed and double-strand GSTe-injected mosquitoes 4 days post-injection p = 0.0024 (GSTe2), 0.0014 (GSTe3), 0.0377 (GSTe4), 0.0422 (GSTe5), 0.0014 (GSTe6), 0.0387 (GSTe7) and 0.0241 (GSTe8). This low expression of all the GSTes in mosquitoes injected compared to the non-injected supports that in vivo injection of dsRNA significantly reduces the expression of GST epsilon genes in the Mibellon An. funestus mosquitoes. actin5C and RPS7 as housekeeping genes. As shown in Figure 3, we noticed a significant partial reduction in GSTes gene expression when comparing control non-exposed and double-strand GSTe-injected mosquitoes 4 days post-injection p = 0.0024 (GSTe2), 0.0014 (GSTe3), 0.0377 (GSTe4), 0.0422 (GSTe5), 0.0014 (GSTe6), 0.0387 (GSTe7) and 0.0241 (GSTe8). This low expression of all the GSTes in mosquitoes injected compared to the noninjected supports that in vivo injection of dsRNA significantly reduces the expression of GST epsilon genes in the Mibellon An. funestus mosquitoes. Figure 3. Confirmation of GSTe knockdown effect by quantitative RT-PCR between non-exposed double-strand injected and non-injected mosquitoes of the same age (a) dsGSTe2, (b) dsGSTe3, (c) dsGSTe4, (d) dsGSTe5, (e) dsGSTe6, (f) dsGSTe7, (g) dsGSTe8,. There is low expression of all GSTe injected mosquitoes compared to the non-injected mosquitoes. dsRNA injection significantly reduces the expression of the GST epsilon genes. * = p < 0.05, and ** = p < 0.01.

Discussion
Understanding of the dynamics of resistance development, the potential of some candidate genes to confer cross-resistance between insecticide classes and designing suitable diagnostic tools are crucial for malaria control. The major genes linked to metabolic insecticide resistance in the major malaria vector An. funestus, include CYP450s, GSTs and carboxylesterase. The characterisation of these major genes provides important information enabling the understanding and the dynamic of resistance development and how and where it spreads, facilitating its management.
In this study, comparative RNAseq-based transcriptomic profiling across Africa showed key differences in the level of expression of GSTs across Africa, including the epsilon, delta and the theta class. These genes are highly expressed in West and Central Africa, in contrast to southern Africa where GSTs are found to be less overexpressed, contrary to P450s, especially the duplicated CYP6P9a and -b genes, which are highly overexpressed in this region [32]. The greater over-expression of P450s genes observed when comparing transcriptomic profiling of DDT resistance across Africa do not necessarily indicate that this enzyme family plays the major role in DDT resistance but could be a result of the multiple resistance observed in these mosquito population with pyrethroids and carbamate resistance as reported. This is supported by the fact that previous studies, such as in An. gambiae, have showed little difference in mortality between bioassays done with DDT alone or after pre-exposure to the synergist PBO [36], showing that P450s may not be the main drivers of DDT resistance but more likely, GSTs in the absence of kdr as seen in An. funestus [37]. The difference in gene expression pattern observed between An. funestus populations is in line with previous studies done in An. funestus [15,32,38,39]. Transcriptional profiling using microarrays/qPCR has established higher overexpression of GSTs in populations from West Africa (Benin) compared with populations from Uganda and Malawi [15]. Among GSTs, the epsilon and delta families are the most expressed as seen also in other mosquito species [13,40], but the GST epsilon cluster is more consistently over-expressed across Africa as previously reported in An. funestus permethrin resistance [32]. All genes of the GST epsilon cluster are observed to be overexpressed in Africa in resistant population of An. funestus, but at different levels, except for GSTe8. This association between overexpression of GSTes and DDT resistance could be explained by the fact that GST plays a role in oxidative stress and its expression is elevated in the response of oxidative damage caused by xenobiotic [41]. Besides cytochrome P450 and GSTs, other gene families were involved in DDT resistance including carboxylesterases, sulfotransferase, ATP-binding factor, UDP-glucuronosyltransferase, and metalloproteinase, and the same pattern were observed in An. gambiae using microarray [42,43].
Silencing of An. funestus GSTe genes supported the role that these genes played in DDT, as well as cross-resistance they confer to pyrethroids. This is in agreement with the findings of Riveron and colleagues [16] who also revealed a cross-resistance between these insecticide classes for GSTe2 using GAL4/UAS transgenic expression in Drosophila and also the association studies of the L119F-GSTe2 genotypes and resistance phenotypes [14,16]. While the work of Riveron and other researchers [18] characterized GSTe2, in particular, in this study, other epsilon GSTs were investigated, confirming their role in DDT and pyrethroid resistance. This cross-resistance could be either by directly metabolising the insecticides or conferring protection from oxidative stress induced by pyrethroids using a mechanism of sequestration [44]. It could also act by detoxifying/scavenging the secondary product generated by reactive oxygen species or by directly metabolizing 4-hydroxynonenal, an end product of lipid peroxidation, through conjugation [45].
This study has shown that knockdown of GSTes in An. funestus significantly increase susceptibility to type I and II pyrethroids, suggesting that the overexpression of GSTes could confer permethrin and deltamethrin resistance. This observation agrees with the results of Lumjuan in 2011, where they proved that partial knockdown of GSTe2 and GSTe7 in Aedes aegypti significantly increased susceptibility to DDT and deltamethrin [46]. This study has also supported that gene silencing through RNAi technique is a good tool to validate the role of candidate genes in insecticide resistance notably for GSTes.

Conclusions
The results have identified genes associated with DDT resistance in the major malaria vector An. funestus Africa-wide. The gene families identified as overexpressed include the cytochrome P450s commonly known to confer resistance to wide range of public health insecticide in malaria vector, carboxylesterase and glutathione-S-transferase, which are the major focus of this study. This has established that GSTe genes are differentially overexpressed in resistant populations of An. funestus across Africa, with the genes consistently more overexpressed in Western and Central Africa compared to East and Southern Africa, consistent with the higher DDT resistance known in the An. funestus populations from West Africa. In addition, the GSTe2, GSTe3, GSTe4, GSTe5 and GSTe6 genes were shown to confer cross-resistance to permethrin, deltamethrin and DDT explaining the multiple resistance observed in the field, highlighting the complexity of resistance and challenges associated with malaria vector control. Further studies need to be performed to detect the genetic variants associated with the resistance conferred by the GST epsilon genes as done previously for GSTe2.

Conflicts of Interest:
The authors declare no conflict of interest.