Insecticide resistance characteristic of Anopheles vector species successfully controlled by deployment of pyrethroid and PBO long lasting insecticidal treated nets in Tanzania

Long lasting insecticidal nets (LLINs) are a proven tool to reduce malaria transmission, but in Africa efficacy is being reduced by pyrethroid resistance in the major vectors. A cluster randomized trial in Muleba district, Tanzania demonstrated that permethrin LLINs co-treated with piperonyl butoxide (PBO), a synergist that can block pyrethroid-metabolizing enzymes in the mosquito, had much greater efficacy than pyrethroid-only nets. Insecticide resistance profiles and underlying mechanisms were investigated in Anopheles gambiae and An. funestus from Muleba during the trial. Diagnostic dose bioassays using permethrin, together with intensity assays, suggest pyrethroid resistance that is both strong and very common, but not extreme. Transcriptomic analysis found multiple P450 genes over expressed including CYP6M2, CYP6Z3, CYP6P3, CYP6P4, CYP6AA1 and CYP9K1 in An. gambiae and CYP6N1, CYP6M7, CYP6M1 and CYP6Z3 in An. funestus. Indeed, very similar suites of P450 enzymes commonly associated with resistant populations elsewhere in Africa were detected as over expressed suggesting a convergence of mechanisms across Sub-Saharan African malaria vectors. The findings give insight into factors that may correlate with pyrethroid PBO LLIN success, broadly supporting model predictions, but revision to guidelines previously issued by the World Health Organization is warranted.


Introduction
The massive scale-up of insecticide treated bednets (ITNs) and in particular long-lasting insecticide-treated nets (LLINs) across sub-Saharan Africa has been the predominant factor in reducing malaria morbidity and deaths since the turn of the century [1]. Unfortunately, the number of malaria cases rose in several African countries in 2016 and 2017, and more widespread resurgence is possible [2]. Funding constraints in the most endemic countries is one factor holding back recent progress [3]. Although less easy to quantify, another key factor in this resurgence is resistance among the vectors to pyrethroids (used for all LLIN treatments [4]), which is now widespread.
Although less common than pyrethroid resistance, resistance to other insecticide classes used in vector control is emerging in many regions of sub-Saharan Africa [6][7][8]. To combat resistance to insecticides, and to pyrethroids in particular, the WHO has developed The Global Plan for Insecticide Resistance Management (GPIRM) [9]. Despite the prevalence of strong pyrethroidresistance, many malaria-endemic countries have yet to align their vector control strategies to those of the GPIRM, in part because of a continued dependence on pyrethroid-treated LLINs [5].
The advent of next generation LLINs -not solely treated with pyrethroids -has been urgently awaited. The first of these bi-treated nets combines a pyrethroid with a non-insecticidal synergist piperonyl butoxide (Py-PBO LLIN). The aim is to improve pyrethroid efficacy, primarily by inhibiting enzymes involved in insecticide detoxification processes [10].
The main mechanisms of pyrethroid resistance in An. gambiae involve mutations to the voltagegated sodium channel (Vgsc) target-site and metabolic resistance [11].
In An. gambiae, over expression of a handful of P450 genes from the CYP6 and CYP9 subfamilies have been repeatedly associated with pyrethroid resistance in field populations of An. gambiae in West and West-Central Africa [11]. When multiple mutations combine they may lead to high-levels of resistance and are likely to seriously threaten the efficacy of malaria control programs [12]. Fewer transcriptomic studies have been performed in East Africa and it is unclear whether a similar concentration of metabolic resistance on a few P450 genes occurs.
Although resistance-conferring Vgsc mutations are absent in An. funestus, metabolic resistance alone seems capable of producing high levels of resistance, which has been associated with control failure [13]. As with An. gambiae, a limited suite of CYP6 and CYP9 pyrethroidmetabolizing cytochrome P450s are involved in pyrethroid resistance in An. funestus, although the relative importance of specific genes varies geographically [14].  [17]. With higher costs and also limitations to supplies, the question of when and where to deploy Py-PBO-LLINs at an operational scale for malaria control is a crucial consideration for international agencies and national malaria control programmes.
Prior to the trial, An. gambiae was the predominant malaria vector in Muleba District with a high frequency of resistance to pyrethroids in diagnostic dose bioassays, attributed at least in part to nearfixation of the Vgsc1014S mutation [19]. Pyrethroid susceptibility in bioassays increased significantly if female mosquitoes were exposed to PBO indicating likely involvement of metabolic resistance mechanisms as well as the Vgsc mutation also referred to as kdr mutation. [20]. However, in Muleba, and in Tanzania generally, studies of specific genes involved in resistance have been limited to Vgsc mutations in An. gambiae [21] and have not been undertaken with An. funestus. Therefore, the present study, was aimed at characterizing molecular and metabolic resistance mechanisms present in intensely pyrethroid resistant populations of An. gambiae and An. funestus, and sought to describe the phenotypic and genetic insecticide resistance profiles of malaria vectors in Muleba.

Study area and mosquito collections
The study was conducted in Muleba district, on the western shore of Lake Victoria in Tanzania.
The area is characterized by high malaria prevalence and the presence of the two major malaria vectors An. gambiae s.s. and An. funestus [16].

Mosquito collection and identification
Indoor resting Anopheles mosquitoes were collected from houses in the different villages between 0600 and 0830 using mouth and/or prokopack aspirators. Blood fed female An. gambiae and An. funestus were kept for three days for blood digestion and used for all WHO diagnostic dose susceptibility and CDC synergist and intensity bottle assays. Three susceptible laboratory strains, namely An. gambiae s.s (Kisumu), An. coluzzii (Ngousso) and An. funestus (FANG) reared at the Liverpool School of Tropical Medicine were used for microarray and qPCR and used as reference strains for gene expression assays. The FANG strain that originated from Southern Angola, colonized at Liverpool since 2002 was used as a reference strain for An. funestus gene expression. For An. gambiae s.s gene expression, we used both Kisumu and Ngousso as reference strains to make analysis more stringent since the Kisumu strain has been colonized in laboratory for so long that it may not resemble a wild population very well.
Morphological identification of the field caught females Anopheles was done according to the key of Gillies and Coetzee [22]. The SINE-PCR method of Santolamazza et al. [23] was used to discriminate members of An. gambiae species complex tested during microarray, while TaqMan assays [24] were performed on a sub-sample exposed to CDC bottle and WHO bioassays.
Member species of the An. funestus group were identified using a cocktail PCR described by Koekemoer et al. [25].

Resistance assay with diagnostic concentration
Wild caught female An. gambiae s.l. and An. funestus were exposed for one hour in WHO cylinders to paper procured from the WHO recommended supplier Universiti Sains Malaysia treated with diagnostic concentrations of either permethrin (0.75%), lambda-cyhalothrin (0.05%), bendiocarb (0.1%) or pirimiphos-methyl (0.25%) [26]. Mosquitoes were then transferred to a holding tube, provided with 10% sugar solution and their mortality was recorded 24 hours later. Approximately 100 mosquitoes (25 per replicate) were used per test. Tests with control mortality exceeding 5% were excluded.

Synergist bioassays
Synergist assays with piperonyl butoxide (PBO) were undertaken to identify the potential role of elevated mixed-function oxidases in resistance in the study area. According to treatment mosquitoes were pre-exposed for one hour, either to 50 μg PBO treated or acetone coated bottles, and transferred for 30 minutes into bottles treated with 21.5µg permethrin or acetone. Mortality was recorded 24 hours post-exposure [27].

Resistance intensity dose response assay
To establish the intensity of resistance a dose response bioassay, using modified Centers for Disease Control and prevention (CDC) bottle bioassays, [26,27] were performed on female An. gambiae and An. funestus collected from Kakoma and Kabirizi villages in April-May 2016.
Wheaton bottles were coated with concentrations of permethrin that gave between 5% and 95% mortality (5 μg/ml to 860 μg/ml for An.gambiae s.l., 21.5 μg/ml to 215 μg/ml for An. funestus and 1.6 μg/ml to 21.5 μg/ml for susceptible An. gambiae Kisumu strain). Approximately 12 mosquitoes were aspirated into each bottle and knock-down recorded at the start and after 15and 30-minutes exposure. Mosquitoes were transferred to paper cups, provided with 10% sugar solution, and mortality recorded after 24 hours. Five to eight replicates were performed for each concentration alongside with a control bottle (coated with acetone).

Genotyping of kdr and GSTe2 mutations
The kdr and GSTe2 mutations were genotyped using F1 progeny from field-collected An. gambiae s.s. and An. funestus. To genotype the nucleotide variants leading to kdr mutations in the An. gambiae VGSC (L1014F or S), hydrolysis probe assays were undertaken as described by Bass et al. (2007) [28]. These used TaqMan primers and minor groove binding (MGB) probes (Applied Biosystems, UK) and SensiMix DNA kit (Quantace).
The qPCR was run on an MxPRO3005 thermal cycler and analysed from endpoint scatter plots using MxPRO software (Aligent technologies, Stratagene, USA). The qPCR cycling conditions for both L1014F and L1014S were 95°C for 10 minutes, followed by 45 cycles of 95°C for 15s and 63°C for 45s.
A further TaqMan assay was used to assess the presence and role of the L119F-GSTe2 mutation, previously associated with DDT and pyrethroid resistance in An. funestus [29]. Reactions were run and analyzed as above with the following conditions. The final volume contained 1× SensiMix (Bioline, London, UK), 800 nM of each primer and 200nM of each probe and the PCR cycling conditions included an initial denaturation at 95°C for 10 minutes, followed by 40 cycles of 95°C for 10s and 60°C for 45s.

Transcriptome analyses
F1 female and male An. gambiae s.s. and An. funestus were separated on the day of emergence and the males were discarded. Females were fed on 10% glucose solution until they were three days old. The mosquitoes were then killed instantly in ethanol and preserved in RNAlater, stored overnight at 4⁰C then transferred to -20⁰C for longer-term storage.

Candidate gene expression analysis
Quantitative PCR analysis was performed on genes identified as candidates from the microarray experiments. Primers were designed using the NCBI primer BLAST [31]. Complementary DNA of cDNA (at 2ng/ul) and 7.8 µl sterile-distilled water. The thermal profile was as follows: 1 cycle 95°C for 3 min, 40 cycles of 95⁰C for 10s then 60⁰C for 10s. Four biological replicates were run for each sample, with three technical replicates. Two endogenous normalizing genes, ribosomal S7 and elongation factor tau, were amplified for each sample to control for variation in cDNA quantity (S1 Table & S2 Table).

Data analysis
Diagnostic dose bioassay results were interpreted according to WHO criteria: A 24 hours mortality superior to 98% indicates susceptibility, mortality of 90 to 97% suspected resistance, and mortality of less than 90% confirmed resistance [26]. To assess resistance intensity, diagnostic concentrations which killed 50% (LC50) of wild An. gambiae and An. funestus were estimated by Probit analysis using Polo plus (version 1.0, Le Ora Software LLC). The LC50 values were used to calculate resistance ratios between each wild population and the An.
gambiae Kisumu susceptible reference strain, which was used as a comparator for both species, owing to unavailability of a susceptible An. funestus strain in the testing laboratory in Muleba.
Loess normalization of microarray data was performed by LIMMA 2.4.1 [32], with analysis of expression using the MAANOVA package [33], with data processing using custom R scripts [14].
Two of the An. funestus arrays were excluded owing to damage to a slide, but the fully interwoven loop design used still provided robust results under such circumstances [34].
Statistical significance of probes was determined using strict criteria based on a multiple testing corrected probability (q-value) threshold (q<0.0001), effect size (fold change, FC >2 or FC <-2) and, for An. gambiae, cross-experiment replication criteria vs. both susceptible colonies. Gene expression for each target gene, was normalized using that of the endogenous genes and was then analysed relative to the susceptible strains of An. gambiae or An. coluzzii (Kisumu and Ngousso) or An. funestus (FANG) using the 2 -ΔΔCt method, correcting for PCR efficiency variation [35].

Resistance assay with diagnostic concentration
In 2014 the mortality rate of An. gambiae exposed to permethrin and lambda-cyhalothrin in WHO cylinder tests was low, with only one of the ten bioassays recording mortality above 10% (Fig 2). Mortality to bendiocarb, the active ingredient used previously for IRS, was much more variable ranging from around 90% in three villages to <35% in Kikagate village. All An. gambiae were fully susceptible to pirimiphos methyl, the active ingredient used for IRS in the trial. The number of An. funestus collected were too few to test at this time.

Resistance intensity dose response assay
In the 2016 collections done in Kakoma, a dose that was 40 times (860 µg/ml) the diagnostic dose permethrin was required to produce a 24-hour mortality of 96.7% (88/91) in CDC bottle bioassays in An. gambiae s.l., and a dose that was 10 times the diagnostic concentration (215 µg/ml) was required to produce a mortality of 97.8% (87/89) in An. funestus (Fig 3). Although confidence intervals overlapped, this suggests a difference in mechanism or frequency of resistance between species, (Table 1). LC50=lethal concentration required to kill 50% of the population, CI=Confidence interval

Molecular mechanisms
For An. gambiae, the 1014S kdr mutation was fixed (N=227), whilst the 1014F mutation was rare, with only one homozygote detected with fewer than 10% F/S heterozygotes. No kdr mutations were detected in An. arabiensis (N=245) [16]. The mutation L119F-GSTe2 potentially linked to resistance was not detected in An. funestus (N=92).
Only An. funestus and An. gambiae were considered in the microarray experiment because An.
arabiensis is a minor contributor to malaria transmission in the study area [36,37].
Transcriptomic analysis of An. funestus from the village of Kabirizi, revealed 789 probes (out of approximately 60,000) as significant using stringent criteria, of which 375 were over-expressed relative to the FANG susceptible strain (Fig 4a). Of these, 48 could be identified as having a possible role in detoxification processes, with 47 of the probes coming from just five P450 genes (each represented by multiple probes). Two P450s, CYP6N1 and CYP6M7, were among the most statistically significant and highly expressed of all genes; expressed at levels 17 and 21-fold higher, respectively, than the FANG strain. The other significant genes represented by multiple significant probes were CYP6Z1 and CYP6M1 (each 4 to 5 times higher than FANG). Each of the genes was confirmed as being significantly over expressed by qPCR, and whilst fold-change values were lower than those in the microarray, CYP6N1 and CYP6M7 remained as the more highly expressed (Fig 4b).
In transcriptome analysis of An. gambiae from Kakoma village, 562 probes were consistently significantly differentially expressed in relation to both of the susceptible strains included (Fig   4a). Of these, 87 (65 of which were over expressed) were identifiable as possible detoxificationrelated genes, including members of the three major metabolic gene subfamilies P450s, GSTs and COEs, as well as transporter genes and alcohol dehydrogenases.
From these significant detoxification genes, P450s were the most common class, with CYP6M2, CYP6Z3, CYP6P3, CYP6P4, CYP6AA1 and CYP9K1 being the most notable (expression range between 4 and 11 higher than the susceptible strains). These P450 genes, along with VATPasethe most over expressed gene (≈88-fold) and also two highly under-expressed detoxification genes GSTE2 and CYP9J5 were chosen for qPCR analysis, along with an additional P450, CYP6M1, that was not significantly over expressed and which was included as a negative control. Excluding the VATPase gene, which showed more than 20-fold lower expression in qPCR than the microarray results, there was a broad agreement between datasets (Pearson's r=0.58, N=18). Most genes were expressed at similar or greater levels in qPCR than in the microarray (Fig 5a). were most strongly and consistently over expressed, followed by CYP6M2 and CYP9K1 (Fig   5b). Each of these An. gambiae P450 genes are known to metabolize pyrethroids.

Discussion
At both baseline and after two years of interventions, An. gambiae showed an extremely high frequency of resistance to both permethrin and lambda-cyhalothrin, with only one instance of permethrin mortality being greater than 10%. In 2017, An. funestus were also tested and showed moderate resistance level (31% mortality  [42,43] and demonstrated to metabolize pyrethroids [44,45].
Evidence that these genes appear key to resistance in An. gambiae has come from West and West Central Africa [11]. The P450s identified in these studies were also the key ones in Muleba, suggesting either convergent evolution or spread of mechanisms between both sides of the continent.
The most significantly over expressed genes in An. funestus were dominated by CYP6 subfamily P450s, most notably CYP6N1, CYP6M7, CYP6M1 and CYP6Z3, all of which have been previously associated with pyrethroid resistance [46,47]. The two best-known pyrethroidassociated P450s in An. funestus, CYP6P9a and CYP6P9b [48] were, however, not over expressed in Muleba, which thus appears to be a population in which CYP6M7 acts in their stead and which metabolizes pyrethroids with equally high efficiency [14]. The gene expression data for the two important Muleba malaria vectors thus provides strong evidence to validate the third WHO criterion, for deployment.
Interestingly, carbamate resistance increased significantly in An. gambiae between 2014 and 2017, despite cessation of carbamate use for control prior to the baseline collections. With acetylcholinesterase target site mutation absent in local An. gambiae [19], pyrethroid driven over expression of P450s might be linked to this change. As in a recent study in Zambia [42], pyrethroid resistance in An. funestus was also accompanied by bendiocarb resistance, suggesting that the mechanism underlying both pyrethroid and carbamate resistance is the same. Crossresistance between pyrethroids and carbamates is associated with CYP6Z1 [46] in An. funestus and with CYP6M2 and CYP6P3 in An. coluzzii, the closest sibling species to An. gambiae [49].
In Muleba the mechanistic profile of both An. gambiae and, especially, An. funestus, which lacks pyrethroid target site mutations, appears dominated by over expression of key candidate P450s.
A crucial point is that genes providing evidence for a P450-based resistance mechanism are not just those shown to be over expressed in the current or even previous studies but have been subject to functional validation to demonstrate their role in pyrethroid detoxification and/or resistance, which is the case for genes in both focal species here. Resistance frequency in Muleba was much greater than 90% in An. gambiae and in the cluster randomised trial Py-PBO LLINs were significantly more effective than standard LLINs [37]. Py-PBO LLINs could, therefore, still be recommended where resistance is greater than 80% and metabolic resistance is prevalent.

Conclusions
A posteriori resistance profiling shows that malaria vector populations in the study area had high resistance intensities yet malaria prevalence in the Py-PBO-LLIN was reduced by 44% compared to standard LLIN [16]. Revision of the condition specified by the WHO for the deployment of Py-PBO LLIN should be considered as evidence on the efficacy of these nets accumulates. to two susceptible strains Note: Two batches of candidate genes were run, so the normalising genes were run for each batch