Transcriptomic analysis of insecticide resistance in the lymphatic filariasis vector Culex quinquefasciatus

Culex quinquefasciatus plays an important role in transmission of vector-borne diseases of public health importance, including lymphatic filariasis (LF), as well as many arboviral diseases. Currently, efforts to tackle C. quinquefasciatus vectored diseases are based on either mass drug administration (MDA) for LF, or insecticide-based interventions. Widespread and intensive insecticide usage has resulted in increased resistance in mosquito vectors, including C. quinquefasciatus. Herein, the transcriptome profile of Ugandan bendiocarb-resistant C. quinquefasciatus was explored to identify candidate genes associated with insecticide resistance. High levels of insecticide resistance were observed for five out of six insecticides tested, with the lowest mortality (0.97%) reported to permethrin, while for DDT, lambdacyhalothrin, bendiocarb and deltamethrin the mortality rate ranged from 1.63–3.29%. Resistance to bendiocarb in exposed mosquitoes was marked, with 2.04% mortality following 1 h exposure and 58.02% after 4 h. Genotyping of the G119S Ace-1 target site mutation detected a highly significant association (p < 0.0001; OR = 25) between resistance and Ace1-119S. However, synergist assays using the P450 inhibitor PBO, or the esterase inhibitor TPP resulted in markedly increased mortality (to ≈80%), suggesting a role of metabolic resistance in the resistance phenotype. Using a novel, custom 60 K whole-transcriptome microarray 16 genes significantly overexpressed in resistant mosquitoes were detected, with the P450 Cyp6z18 showing the highest differential gene expression (>8-fold increase vs unexposed controls). These results provide evidence that bendiocarb resistance in Ugandan C. quinquefasciatus is mediated by both target-site mechanisms and over-expression of detoxification enzymes.

control into MDA programmes can reduce the required number of chemotherapy rounds and consequently the time frame to achieve the microfilaria (MF) prevalence threshold necessary for successful interruption of LF transmission 1,8 .
Although vector control has successfully reduced the burden of vector-borne diseases worldwide, the recurrent and extensive application of insecticides in endemic regions has also triggered an increase in the level of insensitivity to those insecticides approved for public health 9,10 . In addition, the limited number of insecticides available and the occurrence of cross-resistance between different classes is especially worrying for the sustainability of vector control 11,12 . Consequently, identification and monitoring of resistance patterns, and understanding the underlying mechanisms is crucial for extending the lifespan of currently available insecticides, as well as for planning more effective vector control programmes.
Insensitivity to insecticides in arthropods is thought to result mainly through mutations in target-site genes and/or overproduction of detoxifying enzymes 13,14 . Susceptibility studies in C. quinquefasciatus from diverse geographical regions have associated two main target-site mutations to resistant phenotypes. The L1014F mutation in the voltage-gated sodium channel gene, conferring kdr (knockdown resistance), has been associated with pyrethroid and DDT resistance, whilst the G119S mutation in the acetylcholinesterase (Ace-1) gene is linked to resistance to carbamates and organophosphates [14][15][16][17][18] . Metabolic resistance, which involves the over-expression, or increased catalytic capability of metabolic enzymes, is a less tractable mechanism since members of diverse gene families including carboxy/cholinesterases, glutathione S-transferases (GSTs) and cytochrome P450 monooxygenases (P450s) have previously been associated with resistance to different classes of insecticide in a range of vector species [19][20][21] . Over-expression of detoxification genes can be triggered by a range of mechanisms including gene duplication 22 , as observed for the resistance to organophosphates in C. quinquefasciatus mediated by esterases 23 , cis-regulatory elements 24,25 , trans-regulatory elements, or changes in post-transcriptional repression due to differential expression of miRNAs 26 .
Due to the diversity of genes or gene families involved in metabolic resistance, identification of candidate genes requires an agnostic survey of the patterns of gene expression associated with resistant phenotypes. Recently, studies have applied either microarray or RNA-Seq platforms to elucidate the relationship between gene expression and insecticide resistance [27][28][29] although to date, most of these whole-transcriptome studies in vector insects are restricted to mosquitoes of the genus Anopheles. Despite the role of Culex in transmission of several pathogens such as filarial worms and West Nile virus (WNV) 30 , and reports of high levels of insecticide resistance, few studies have addressed the relative impact of metabolic resistance in C. quinquefasciatus 31,32 particularly at a whole-transcriptome scale (although see 33,34 ).
In addition to their role as disease vectors, C. quinquefasciatus are nuisance biters and failure to effectively control this species can lead to the perception of malaria control failure which may ultimately lead to the rejection of controls (e.g. IRS and ITNs) 16,35 . In this study, we report the results of C. quinquefasciatus susceptibility bioassays for six insecticides (DDT, permethrin, deltamethrin, bendiocarb, fenitrothion and lambda-cyhalothrin) from Nagongera, Tororo District, Uganda. We then report and apply a novel 8 × 60K whole-transcriptome microarray to identify candidate genes associated with bendiocarb insecticide resistance, the active ingredient used in IRS control in Uganda at the time of collection 36 .

Materials and Methods
Sample collection. Mosquitoes were collected in Nagongera, Tororo, Uganda (0° 47′ 48.9978″, 33° 58′ 47.1″) between June and July 2012. Resting adult C. quinquefasciatus were collected exclusively inside houses using aspirators and transported to the insectary. From these collections, caught blood-fed females were maintained in individual Eppendorf tubes lined with moist filter paper to encourage egg laying 37 . From these, 64 females laid at least one egg batch. All egg batches were pooled and then floated simultaneously across five water-filled trays and emergent larvae fed on Tetramin fish food. Pupae from all trays were transferred to a single cage and adults allowed to emerge. Cages for emerged adults were changed every 3 days so that single cages contained only 3-5 day-old adults. All bioassays were performed over a single 10 day period on adults that emerged from the primary egg batches. Adult mosquitoes were fed ad libitum on 10% glucose solution and used for insecticide susceptibility testing. Genomic DNA from each female from which egg rafts were obtained to found the colony was individually isolated using a DNeasy kit (Qiagen) then used for identification of C. quinquefasciatus using a diagnostic PCR assay 38 . In addition to these field-collected mosquitoes, a laboratory colony of C. quinquefasciatus from the Tropical Pesticides Research Institute (TPRI) Tanzania was used as a susceptible reference strain 39 for the microarray study. The TPRI colony susceptibility to the insecticide bendiocarb (see below) was verified before the transcriptomic profiling by exposing 100 females (four replicates of N = 25) to bendiocarb in WHO assays (see below).
Insecticide susceptibility test. Bioassays were performed using test kits and insecticide-impregnated papers according to standard WHO methods (WHO) 40 , which briefly consist of exposing mosquitoes to insecticide impregnated papers fitted inside plastic tubes (125 mm × 44 mm) for 1 hour. Mosquitoes' susceptibility is assessed by the mortality rate after a 24 hour recovery period. For each insecticide tested, mosquitoes were assayed in four replicates, performed on consecutive days, using 25 non-blood fed females (3-5 day-old) emerged from the F1 generation. Mosquito rearing and bioassays were conducted under insectary conditions: 27 °C + 2 °C, 75% ± 10% relative humidity and 12 h:12 h light:dark cycles. Tests were performed with papers impregnated with diagnostic concentrations of six insecticides: DDT (4%), permethrin (0.75%), fenitrothion (1%), lambda-cyhalothrin (0.05%), bendiocarb (0.1%) and deltamethrin (0.05%). Mosquitoes were exposed for 1h, with the exception of bendiocarb and deltamethrin where, following the results of initial 1h exposures, additional four-hour exposures were used to increase the discrimination between resistant and unexposed sympatric mosquitoes. Control assays were performed with 25 mosquitoes exposed to non-insecticide treated papers.
Insecticide exposed mosquitoes were then transferred to clean holding tubes and provided with 10% glucose for a 24-hour period after which mortality was recorded, with dead (susceptible) mosquitoes collected and individually stored on silica gel (for downstream DNA analyses) whilst alive (resistant) mosquitoes had a hind leg removed (stored on silica) and the whole body stored in RNAlater (Sigma Aldrich). Specimens preserved on silica were used for DNA genotyping for the G119S mutation in Ace-1 while RNAlater was used to preserve RNA integrity before total RNA extraction, as described below. RNAlater stored mosquitoes were initially held overnight at 4 °C to allow the solution to penetrate the carcass before storage at −20 °C until. Synergist assays. Synergist tests were carried out using mosquitoes from the same generation utilised in the standard insecticide resistance assays following the procedure described above with an additional pre-exposure to three synergist compounds in order to investigate the potential mechanisms of metabolic resistance for bendiocarb and deltamethrin insecticides. For each synergist assay, four batches of 20-25 mosquitoes were either pre-exposed to impregnated papers (12 cm × 15 cm, Whatman grade no.1 filter paper) with 4% PBO (piperonyl butoxide -CYP450 synergist), 10% TPP (triphenyl phosphate -esterase synergist) or 8% DEM (diethyl maleate -GST synergist) for one hour, then exposed to bendiocarb (0.1%) or deltamethrin (0.05%) for four hours, followed by a 24h recovery period. Synergist only controls were run simultaneously.
Ace-1 genotyping of bendiocarb phenotyped mosquitoes. Prior to microarray analysis all mosquitoes were genotyped for the G119S mutation in Ace-1. DNA was isolated from the amputated hind leg of resistant mosquitoes with 50 µl of 10% Chelex 100 and 2 µl of proteinase K (10 mg/mL). The homogenate was incubated at 94 °C for 30 min followed by centrifugation at 6,000 rpm for 10 min to collect the supernatant.
A TaqMan assay, designed to detect the G119S mutation in the acetylcholinesterase gene 41 , was utilised for genotyping with reaction mixtures composed of 1 µl of genomic DNA, 1x SensiMix II probe, 400 nM of each primer and 100 nM of each probe in a final volume of 50 µl. Thermocycling was performed on the Stratagene MX3005P and consisted of 95 °C for 10 min and 40 cycles of 92 °C for 15sec and 60 °C for 1 min with endpoint discrimination.
Microarray. 8 × 60 K microarray construction and study design. The pattern of gene expression of bendiocarb resistant mosquitoes was investigated using a custom designed C. quinquefasciatus whole genome oligonucleotide microarray (Agilent ID 039759). An 8 × 60K microarray format with 60-mer oligonucleotide probes was designed to cover a variety of targets using eArray (http://earray.chem.agilent.com/earray/). The majority of this array encompassed three probe replicates (56,598 probes) for each of the 19,018 transcripts in the CpipJ1.3 gene build 42 . Additionally, we downloaded 205,396 ESTs from VectorBase. From these, we identified 1,987 contigs and 4,109 singleton sequences and designed a single probe for each (from these 1,987 contigs, only 1,935 had probes successfully designed and from 4,109 singletons 2,862 had successful probes designed). Additionally, three probe replicates were designed to each of four alternative GST transcripts not annotated in VectorBase 43 , and 25 additional replicated probe groups (10 replicates) to allow estimation of reproducibility (CV (coefficient of variation probes) (Fig. 1a). Nevertheless, it is important to bear in mind that the current array version encompassed genomic information exclusively from the CpipJ1.3 annotation and may therefore lack the resolution to pinpoint Interwoven hybridization loop design for comparison between bendiocarb exposed and non-exposed Ugandan field-collected mosquitoes and the TPRI susceptible strain. Circles represent pools of 10 females. C: Uganda non-exposed mosquitoes (sympatric control), R: Uganda Resistant mosquitoes, TPRI (Tropical Pesticides Research Institute): C. quinquefasciatus susceptible strain from Tanzania. specific gene sets linked to geographic background. Full details of the array design are given in ArrayExpress (http://www.ebi.ac.uk/arrayexpress/) with the accession number A-MTAB-649.
From the six insecticide susceptibility test groups, we chose bendiocarb selected mosquitoes for the microarray analysis since we observed only moderate mortality 58.02% (95% CI-confidence interval 51.53-64.25%) for this insecticide following 4h exposure ( Fig. 2) thereby allowing clear discrimination between resistant and unexposed sympatric mosquitoes and we also detected an increase in mortality in the presence of two synergists (See Fig. 2; TPP and PBO), indicating the likely involvement of metabolic resistance. Three experimental conditions were employed: Tororo_Resistant (following 4h bendiocarb exposure), Tororo_Control (4h exposure to control papers) and TPRI_Control (4h exposure to control papers). All mosquitoes were wildtype 119G homozygotes for Ace-1. Comparison between the bendiocarb selected samples and controls (Tororo_Control and TPRI_Control) was performed on three RNA pools per group using an interwoven loop design (Fig. 1b) as described by Vinciotti et al. 44 .
RNA extraction, labelling and hybridization. Total RNA was isolated from three pools of 10 female mosquitoes for each group using the RNAqueous-4PCR kit (Ambion) according to the manufacturer's instructions. Total RNA quantity was assessed using a Nanodrop spectrophotometer and RNA quality assessed on an Agilent Bioanalyzer. Each pool of RNA was individually labelled with Cyanine-3 and Cyanine-5 (Cy3 and Cy5) using the Low input Quick Amp Labelling Kit (Agilent Technologies) followed by purification through Qiagen RNeasy Columns (Qiagen) with quality and quantity checked using a Nanodrop and Bioanalyzer, respectively.
Before hybridization, 300 ng of Cy3 and Cy5 labeled cRNA was fragmented using the gene expression hybridization kit (Agilent) in a total volume of 25 µl including 5 µl of 10x blocking agent and 1 µl of 25x fragmentation buffer. The fragmentation reaction was incubated at 60 °C for 30 min, then chilled on ice for 2 min before addition of 25 µl of 2 x GE hybridization buffer Hi-RPM. Each array was hybridized using 45 µl of the fragmentation solution for 17 hours at 65 °C and 10 rpm. After hybridization, arrays were washed with wash buffers 1 and 2 for 1 min each, followed by acetonitrile for 10 sec and finally fixation solution for 30 sec. Arrays were scanned using the microarray scanner system (Agilent Technologies) and feature extraction performed using Feature Extraction software (Agilent Technologies) according to the manufacturer's recommendations. All arrays passed the Agilent quality control with QC score ≥10.
Data analysis. All microarray data analysis was performed in R version R.3.0.1 45 . Array normalization was carried out using the Limma 3.2.3 package 46 and data analysis performed using MAANOVA software version 3.0 47 to detect overall differential levels of gene expression across the three treatment groups. The top overexpressed genes were selected after an ANOVA F-test based upon a false discovery rate (FDR) decision criteria of log 10 (Q value) > 2.5. Within this subset of significantly differentially expressed probes, those that were significantly overexpressed in the Bendiocarb-exposed mosquitoes were pinpointed by comparing expression pattern in all pair-wise comparisons (see Table 1).
Additionally, the pattern of differential transcript level was also characterized through pair-wise comparisons between Ugandan resistant, sympatric control and the TPRI strain using a Student's t-test with genes considered differentially overexpressed where P < 0.05.
Functional characterization of differently expressed transcripts detected from the ANOVA and pair-wise comparison were then submitted to a Gene Ontology (GO) analysis to classify probes on their GO categories (cellular components, biological process and molecular functions). For this, a sub-set of overexpressed probes were selected based on the threshold of log 10 (Q value) > 2 of the ANOVA analysis to capture a broader selection of highly overexpressed probes. This was then submitted to VectorBase for functional annotation (GO term identification) using the Biomart tool. The REVIGO web server 48 was employed to summarize and visualize the distinct GO terms identified. Relatedness among GO terms was assessed using the uniqueness method, followed RT-qPCR validation. Reverse-transcription quantification PCR (RT-qPCR) using cDNA synthetized from the same pools of total RNA applied for microarray hybridization were used to confirm the expression profile of three out of 16 top candidate gene identified by the microarray (CPIJ020018 [=Cyp6z18], Cyp6n23, R2D2 Table 1) using two housekeeping genes -40S ribosomal protein S3 and β-tubulin -as endogenous controls. Specific primers to amplify PCR fragments with size ranging from 107 to 176 bp (see Supplementary Table S1) were designed using primer3 software 51 ; however, only for four genes (Cyp6z18, R2D2, 40S ribosomal protein S3 and β-tubulin) was it possible to design primers spanning exon junctions (primers for Cyp6n23 were exonic only). Specificity of primer sets was verified by identification of a single symmetrical amplicon peak following melting curve analysis. Additionally, PCR efficiency was verified using a 10-fold serial dilution of standard cDNA with only primer sets with an efficiency ranging from 90 to 110% taken forward for RT-qPCR reactions.
Two technical replicates of all RT-qPCR reactions were carried out for each gene in a total volume of 20 µl including 1 µl of cDNA (1:4 stock diluted), 10 µl of Brilliant II SYBR ® master mix (Agilent Technologies) and 100 nm of each forward and reverse primer. Amplification was conducted under standard qPCR reaction conditions on the Mx3500P qPCR system (Agilent Technologies). Gene expression quantification of the three selected genes was assessed according to the ΔΔ Ct method 52 .
Manual annotation of the CPIJ020018 gene region. During qPCR primer design for the top candidate gene CPIJ020018 (annotated as Cyp6z16 in VectorBase CpipJ1.3 assembly #1, annotation 1.3), we concluded that the available, automated gene annotation of this particular gene was unreliable. The genomic sequence of this region in VectorBase includes a region of about 810 bp in supercontig 3.2948 with no nucleotide sequence information, which spanned the automated annotation of Cyp6z16. We suspected additional coding sequence to lie within this region. To confirm this, we designed primers to span the complete region and amplified a 4.7 kb region covering the full length of the candidate gene genomic sequence. PCR reactions to amplify the Cyp6z16 genomic region were conducted in a final volume of 20 µl including 40 ng of genomic DNA, 1 x Phusion HF buffer, 200 µM each dNTP, 0.5 µM each primer Cx_6Z16F and Cx_6Z16R (see Supplementary Table S2) and 0.02 U/µl Phusion DNA polymerase. Reaction conditions were 98 °C for 30 sec, 30 cycles of 98 °C for 10 sec, 62 °C for 30 sec and 72 °C for 3 min, with a final extension of 72 °C for 5 min. PCR products were purified using the GeneJET PCR purification kit (Thermo Scientific) then cloned into pJET1.2 PCR vector (Thermo Scientific). Finally, the full length of cloned PCR product after plasmid purification using the GeneJET Plasmid Miniprep Kit (Thermo Scientific) was sent to Source BioScience, UK for Sanger sequencing using nine internal primers (see Supplementary Table S2). Sequence traces were analyzed using CodonCode Aligner version 4.2.2. Following removal of vector sequences a single contig was built from overlapping sequences. This contig sequence was then used for gene structure prediction and transcript annotation using the Augustus web interface 53 .

Gene
Transcript GO www.nature.com/scientificreports www.nature.com/scientificreports/ Transgenic expression in Drosophila flies. The full-length sequence of Cyp6z18 was codon-optimised by Genscript (Piscataway, NJ, USA) for Drosophila and synthesised with EcoRI-XbaI sites then cloned into the pUAST.attB vector (provided by Dr J Bischof, University of Zurich). Transformation was performed using the PhiC31 system with plasmids injected into the germline of D. melanogaster with a chromosome 2 attP landing site (y w M(eGFP, vas-int, dmRFP)ZH-2A; P{CaryP}attP40) by the University of Cambridge Fly Facility. A single transgenic line was generated and balanced. To induce expression of Cyp6z18, flies were crossed to the Act5C-GAL4 strain (y1 w*; P (Act5C-GAL4-w) E1/CyO,1;2) (Bloomington Stock Center, IN, USA). For each treatment, Cyp6z18 transgenic or untransformed controls, three crosses of 12-15 females to 6-7 Act5C-GAL4 males were performed. Shortly after pupae were seen, the parental generation was transferred into new rearing vials to continue egg laying. The date of appearance of newly-emerged offspring was observed for a total of 8 Cyp6z18 and 11 control rearing vials.

Frequency of Ace-1 resistant alleles in bendiocarb selected mosquitoes.
Both dead and alive mosquitoes following exposure to 0.1% bendiocarb (4h) were genotyped for the Ace1-119S mutation using a custom Taqman assay. The wild-type allele was observed at the highest frequency (Fig. 3a) (Fig. 3b). There was a highly significant association between the Ace1-119S allele and bendiocarb resistant phenotype ( Fig. 3c; P < 0.0001, Fisher's Exact test) with an OR of 25 (95% CI 3.37-186).
Gene expression profiling of bendiocarb selected mosquitoes. To identify candidate genes associated with resistance in bendiocarb selected mosquitoes, transcriptomic profiles of three groups of Ace-1 wild-type (119G) mosquitoes were compared: Ugandan bendiocarb exposure survivors, Ugandan Control exposed (sympatric control) and a control exposed fully susceptible TPRI strain. Among the three groups we identified 32 probes significantly differently expressed in the ANOVA analysis applying a threshold of -log 10 false-discovery rate (FDR) adjusted P value > 2.5, of which 50% of the probes had higher expression in the Ugandan exposed compared to both the sympatric control and TPRI ( Table 1). The list of top candidate genes included two P450s: www.nature.com/scientificreports www.nature.com/scientificreports/ CPIJ020018 (shown to be Cyp6z18 following reannotation -see below) and Cyp6n23, which displayed increases in gene expression compared to the TPRI strain of 8.75 and 4.37-fold, respectively.
For functional characterization of the most significantly differentially expressed probes, a list of probes with a − log 10 false-discovery rate (FDR) > 2 obtained from the ANOVA (Fig. 4a) comparison was submitted to Gene Ontology (GO) analysis. After removing duplicate probes, 358 unique genes were submitted to VectorBase for the GO term search using Biomart (for complete annotation see Supplementary Table S3.) In total 298 terms were obtained with the majority of extracted GO terms clustered on the molecular function and biological process categories (see Supplementary Table S4 for GO Term frequencies and description). Among the top 10 enriched terms for each category, GO terms linked to metabolic process, biosynthetic process, transport, membrane, integral component of membrane, nucleus, binding, cation binding and metal ion binding were observed with percentage of annotations ranging from 2.16 to 75.39% (1-81 genes, respectively) ( Fig. 4b-d). Metabolic process was the predominant term across all categories corresponding to 75.39% of enriched GO terms, while within the category of molecular function and cellular component we observed a large proportion (21.23%) of terms associated with binding functions (29 genes) and membrane 61.59% (29 genes), respectively. Further analysis utilising the 358 unique genes showed five clusters of GO terms with significant enrichment (enrichment score >1.3: Supplementary Tables S5 and S6). Cluster 1 had the highest enrichment score (5.15) and was associated with the GO term 'transferase activity' and Interpro IDs for Glutathione S-transferases ('Glutathione S-transferase, N-terminal' Benjamini-Hochberg P = 6.7 × 10 −6 ; 'Glutathione S-transferase, C-terminal' Benjamini-Hochberg P = 2 × 10 −5 ).
To explore further the transcriptomic profile of the bendiocarb resistant phenotype, pairwise comparison using Student's t-test was applied to compare exposed and unexposed Ugandan mosquitoes to the TPRI susceptible strain. Significantly up and down-regulated genes with a fold-change > 2.5 were also investigated by GO analysis. For down-regulated genes we observed similar figures between exposed and unexposed mosquitoes in contrast to up-regulated genes where we identified 33 genes exclusively in the exposed mosquitoes (Fig. 5a). Pairwise comparisons also identified eight significantly differentially expressed genes putatively associated with insecticide detoxification: three exclusively in the pools of exposed mosquitoes: GSTs (CPIJ018629-RA; fold-change 2.79 and CIPJ018632-RA, fold-change 2.0) and esterase (CPIJ013918-RA; fold-change 2.86) whereas for the other five: P450 (CPIJ020018-RA), GSTs (CPIJ010814-RA, CPIJ018624-RA, CPIJ018626-RA) and esterase (CPIJ013917-RA) were observed in both Ugandan exposed and unexposed mosquitoes.
In general, the GO term enrichment of the up and down-regulated genes for the Ugandan exposed and unexposed mosquitoes shows a similar list of terms (Fig. 5b) with the exception of three terms: metabolic process, phosphate-containing compound metabolic process and cation binding that were observed exclusively on the exposed mosquitoes with frequencies of 60.61, 16.39 and 17.07%, respectively (Fig. 5c). Nevertheless, this analysis excluded 39 and 18 up and down-transcribed genes (See Supplementary Tables S7 and S8), respectively, such as esterases (CPIJ013917-RA and CPIJ013918-RA) and heat shock proteins (CPIJ013880-RA and CPIJ005642) for instance, that were not associated with GO terms from VectorBase. www.nature.com/scientificreports www.nature.com/scientificreports/ Candidate gene validation by qPCR. The expression patterns of three candidate genes, randomly chosen from the list of top candidates, (both P450 genes Cyp6z18 and Cyp6n23, plus R2D2) were additionally assessed by qPCR. Satisfactory PCR efficiency ranging from 91.4 to 101.8%, (within the 10% acceptable variation) was obtained for primer pairs designed for candidate genes and endogenous controls ( Supplementary Fig. S1). Additionally, the primer sets were specific, with a single symmetrical amplicon peak in melting curve analyses ( Supplementary Fig. S2). For all three candidate genes (Cy6z17, Cyp6n23 and R2D2) we observed a good correlation between the microarray and RT-qPCR expression fold-change with ratios of up-regulation detected by both methods differing by less than 1.5 × (Supplementary Fig. S3).

Annotation of CPIJ020018.
After identification of CPIJ020018 as the top candidate gene in the microarray results, further analyses on the genomic and cDNA sequences available from VectorBase were conducted. This revealed an atypical gene architecture for a P450 gene with the presence of one intron >3Kb (Fig. 6a). Closer analysis of the genomic sequence indicated a region of 809 bp with no nucleotide information (Ns) internal to the CPIJ020018 genomic sequence. In silico analysis of a contig constructed after PCR amplification and cloning of the complete region encompassing the genomic region of interest (accession number MH822866), suggested the presence of two distinct P450s instead of the single Cyp6z16 predicted by the automated annotation in VectorBase. The two genes predicted are each composed of two exons separated by one intron (Fig. 6b and Supplementary Material 1). BLAST analysis of both predicted amino-acid sequences against Anopheles gambiae and Aedes aegypti sequences available on VectorBase and the Cytochrome P450 homepage 54 shows the top hits belong to CYP genes from the Z family. Following submission to the P450 nomenclature database the partial P450 is now labelled cyp6z16 and the full-length novel P450, which is interrogated by the microarray probes (hitting CPIJ020018 exon-2) is formally labelled cyp6z18 (note cyp6z17 is annotated is annotated on supercontig 3.3058). Together the gene prediction and BLAST analysis carried out indicate that the top candidate gene from the microarray analysis is cyp6z18 and not cyp6z16 (which is not interrogated by any probes due to incorrect annotation). cyp6z18 is closely related to other P450s which have previously been associated with metabolic resistance in both An. gambiae 55,56 and Ae. aegypti 57,58 (Fig. 6c). We note that whilst our microarray design was undertaken on the CpipJ1 assembly, the CpipJ2 assembly (available from April 2014) did not revise the assembly or annotation of this scaffold.
Transgenic expression of cyp6z18 in Drosophila flies. Drosophila melanogaster flies were transformed with a transgene encoding cyp6z18 for expression using the GAL4/UAS system. Individuals containing the www.nature.com/scientificreports www.nature.com/scientificreports/ cyp6z18 gene and untransformed controls were crossed to flies of the Act5C-GAL4 driver line to induce ubiquitous expression of cyp6z18. Expression of cyp6z18 negatively impacted fly fitness, evident in prolonged pre-adult development. The number of days between placing adult flies in a rearing vial and the first offspring emerging was measured. A significant difference in this value was detected, with an average of 23 days for cyp6z18-expressing flies (95% CI 21.21-25.04) and 16 days for untransformed control flies (95% CI 15.23-16.41) (unpaired t-test with Welch's correction P < 0.0001, N = 19). Because of this fitness cost, experiments to test the impact of cyp6z18 expression upon insecticide resistance were not pursued.

Discussion
Insecticide resistance is widespread in malaria vectors in Uganda 37,59-63 and known to result from both target-site and detoxification mechanisms [63][64][65] . In this study, we investigated the insecticide susceptibility status and possible mechanisms involved with insecticide resistance in the LF vector C. quinquefasciatus mosquitoes collected in a region of Uganda (Tororo) where both Anopheles gambiae and Anopheles funestus exhibit resistance to the insecticides currently used for malaria control 37,62,63 . We have shown previously 41 that the pattern of target-site mutations in Ugandan Culex varies regionally despite intense gene flow, suggesting a heterogeneous pattern of insecticidal selection pressures. Here we show that in Culex from Tororo, eastern Uganda, resistance to three of the four classes of insecticide occurs at high levels and, at least for bendiocarb, is mediated by both target-site and metabolic mechanisms. Although, the G119S mutation in Ace-1 is not at high frequency, it is strongly associated with resistance with an OR of 25 (95% CI 3.37-186). With such a strong resistance association it seems surprising that the 119S allele is at such low frequency 14% (95% CI 7.9-22.4) in the Nagongera population. This is a much lower frequency than observed in other vector mosquitoes where carbamates are routinely used in vector control 66 . Whilst insecticides linked to selection of Ace-1 resistant alleles have not been officially applied for vector control in the studied area for over one decade 67,68 indoor residual spraying (IRS) using bendiocarb was conducted in Tororo from December 2014 -February 2015 69,70 , after the period of collection of these samples, and hence this allele frequency is higher than expected. Nevertheless, an indirect source of insecticide exposure due to agricultural activities could also be shaping the insecticide resistance as described for C. quinquefasciatus populations from Iran 71 . However, the Ace-1 mutation is known to confer a fitness disadvantage which may explain the low frequency of the resistance associated allele 72 . This recent evolution of insecticide resistance in Ugandan Culex populations, possibly in response to malaria vector control, has also been suggested in other www.nature.com/scientificreports www.nature.com/scientificreports/ African countries such as Tanzania and Zambia, where moderate resistance has also been detected to bendiocarb, in contrast to the high levels of resistance detected to both pyrethroids and DDT 16,73 .
The potential involvement of metabolic resistance in bendiocarb resistant C. quinquefasciatus from Tororo is supported by the synergist assays, which show an increase in mortality to bendiocarb when mosquitoes were pre-exposed to either TPP (a synergist of esterases) or PBO (a synergist of cytochrome P450s). Together the Ace-1 genotyping and synergist assay data strongly suggested an alternative mechanism of resistance to the well-known Ace1-119S target-site mutation in the Ugandan bendiocarb resistant phenotype. Transcriptomic profiling of Ace-1 wild-type mosquitoes identified two P450s (CPIJ020018 and Cyp6n23) with the highest up-regulation in the resistant samples among the top candidate genes. Over-expression of both these cytochrome P450s (CPIJ020018 and cyp6n23) identified in our analysis are especially relevant as many genes belonging to this gene family have been associated with insecticide metabolism in a variety of vector species 28,74,75 .
The likely association of this candidate cytochrome P450s with the carbamate resistant phenotype is supported also by the synergism effect of the P450 inhibitor PBO 76 . Most CYPs previously associated with the insecticide resistance phenotype in mosquitoes, (e.g. cyp6p3, cyp9j32 and cyp6m10) are typically associated with metabolism of pyrethroids and DDT 56 (cyp6z1 and cyp6z2) have also been demonstrated to be capable of metabolizing the carbamate insecticide carbaryl 79 and upregulation of the An. funestus cyp6z1 is associated with carbamate resistance in this species 80 . Interestingly, both cytochrome P450s identified here belong to the CYP6 family, which includes most of the CYPs genes already described as insecticide metabolizers 23 . Our GO term enrichment analysis also indicated a likely impact of other metabolic genes such as GSTs, which from our data had the highest enrichment score in the cluster analysis, with two GST genes (CPIJ018629-RA and CIPJ018632-RA) overexpressed exclusively in exposed mosquitoes. These findings indicate that further functional analysis should be performed to pinpoint the role of GSTs in the carbamate resistant-phenotype of Culex. We note that other studies have already shown that genes belonging to this family have roles in resistance to DDT and pyrethroids in Culex mosquitoes 81,82 .
For the candidate gene cyp6z18, functional enzymatic characterization for insecticide metabolism, was applied as an in vivo approach using transgenic gene expression in Drosophila. Unfortunately, this approach did not succeed in validating this gene as the transgenic construction negatively affected Drosophila development. Further work using alternative approaches such as in vivo gene expression through CRISPR/Cas9 83 or in vitro metabolism assays e.g. 28 could be applied to confirm the role of the candidate gene in metabolic resistance. As shown herein and elsewhere 84,85 , investigation of resistant phenotypes in light of metabolic resistance is imperative for pinpointing the genetic mechanisms associated with evolving insecticide resistance. For instance, while our data show that both target-site and detoxification gene over-expression underly the resistance to bendiocarb in C. quinquefasciatus from Tororo, other studies have demonstrated that P450 overexpression alone is the primary mechanism in kdr-free populations such as in Anopheles arabiensis from Chad, Central Africa 84 and Malaysian Aedes albopictus 85 . Taken together, these findings highlight the relevance of identifying the mechanisms driving resistance in order to facilitate the development of field-applicable markers to detect increased resistance at the early stages of resistance development as well as to assist in decision-making for more effective control interventions. For instance, choosing either insecticide-only nets or insecticidal-synergistic nets (e.g. those incorporating piperonyl butoxide -PBO) 86 or the rational use of insecticides from distinct classes to tackle the threat of an increased burden of vector-borne diseases driven by reduced efficacy of ongoing vector control interventions 87 .
Further analysis of the genomic and transcriptomic sequences of the top candidate gene CPIJ020018 (annotated as cyp6z16 in VectorBase) suggests annotation inaccuracy. After re-annotation, our analysis indicates that CPIJ020018 consists of two distinct CYP genes (cyp6z16 and cyp6z18) instead of one as suggested by the automated annotation of VectorBase. Annotation problems as observed for cyp6z16 may potentially occur for other genes (see 88 for annotation deficiencies of recent Anopheles sequencing), thus a review and improvement of gene annotation for C. quinquefasciatus and other vector species such as Anopheles dirus, Aedes albopictus and Lutzomyia longipalpis recently included on the VectorBase genome browser, especially of gene families linked to insecticide resistance, would benefit future transcriptomic monitoring of insecticide resistance. Whilst, the increased usage of RNA-Seq for differential gene expression analyses may suggest this is unnecessary, the typical RNA-Seq pipeline maps reads to annotated gene sets and hence a robust annotation of these gene families would still be beneficial. C. quinquefasciatus remains the poor relation of disease vectors with respect to assembly and annotation of genome sequence, whilst the genome sequences of A. gambiae and A. aegypti have received important and successful genome updates 89,90 Culex remains assembled with a high number of scaffolds and hence gaps.
This microarray study also identified that four other detoxification genes (three glutathione S-transferases (GSTs) and one esterase) were also over-transcribed compared to TPRI, although they were not differently expressed between Tororo resistant and sympatric unexposed controls. High expression of these genes alone or in combination with other detoxification enzymes could also be a possible mechanism associated with the bendiocarb phenotype as reported previously 91,92 , although in our data no significant synergist effect by DEM, a GST inhibitor was observed, further characterization of these genes is warranted.
Although C. quinquefasciatus is a vector of important neglected tropical diseases, planning and management of vector control strategies for Anopheles species has received considerably more attention and funding. Whilst this could reflect that MDA is currently the primary intervention for LF eradication with vector control deemed secondary, a reduction of LF transmission by Culex species due to control programs aimed at anopheline species has been demonstrated 93,94 . Thus, in a possible integrated vector-control scenario the monitoring of insecticide resistance and determination of the mechanisms resistance for vector species of apparent secondary interest can be critical to effective integration 95 .