Introduction

The introduction of synthetic insecticides is often followed by the appearance of resistance phenotypes in field populations, leading to significant reductions in agricultural production1. There has been much debate about whether the evolution of resistance is caused by variation in a single gene (monogenic) or by the additive effects of many (polygenic)2, 3. Substantially more work has been dedicated to characterizing the monogenic variants, but such alleles arise in a genetic background where there is polygenic variation for tolerance to the insecticide2. Much still remains unclear about the relative contribution of different alleles to insecticide resistance, but D. melanogaster is uniquely placed to answer such questions, owing to the extensive genetic toolkit that has been developed in this model insect.

Imidacloprid is amongst the most widely used insecticides. It is derived from nicotine and is a member of the neonicotinoid chemical class. Neonicotinoids target nicotinic acetylcholine receptors (nAChRs) that have vital roles in neurotransmission and behaviour in insects4, 5. However, imidacloprid resistance via mutations in targets is not the most common resistance mechanisms, possibly due to associated fitness costs6. The overexpression of cytochrome P450 enzymes (P450s) more frequently underpins imidacloprid resistance7. Some members of the P450 superfamily can function as drug metabolizing enzymes (DMEs) with xenobiotic substrates, while others have vital roles in development using endogenous substrates8. P450s which are capable of metabolizing imidacloprid and conferring resistance have been identified in several species9,10,11; the Cyp6g1 gene of D. melanogaster has been particularly well studied12. Originally identified by mapping DDT resistance in the Hikone-R strain to a region containing a cluster of three P450 genes (Cyp6g1, Cyp6g2 and Cyp6t3), resistance was shown to be due to the overexpression of Cyp6g1 13. This overexpression was subsequently found to be caused by the insertion of the long terminal repeat of the retrotransposon, Accord, upstream of the gene14.

The expression of Cyp6g1 is highly variable in field populations due to the Accord insertion, copy number variation and further transposable element insertions15,16,17. The ancestral M haplotype contains a single copy of Cyp6g1 and expresses low levels of the gene compared to the more derived AA haplotype, which contains a duplication of the Accord-Cyp6g1 cassette in addition to several partial chimeric repeats of Cyp6g1-Cyp6g2 which are not characterized16. Further modifications of the AA haplotype resulted from the insertions of the transposable element HMS-Beagle and a P element upstream of Cyp6g1, creating the BA and BP haplotypes, respectively. These derived haplotypes have been associated with increased levels of Cyp6g1 expression and resistance to insecticides such as DDT and azinphos-methyl16, 17. Cyp6g2 expression correlates with Cyp6g1 expression in the Drosophila Genetic Reference Panel, but the contribution of this gene to resistance has not been shown17.

The structural modification of imidacloprid in biological systems includes both nitroreduction and oxidation reactions. Metabolites from both pathways have been detected in plants, animals and insects, but soil bacteria produce predominantly the nitroreduction metabolites18,19,20. In the case of D. melanogaster, the presence of nitroreduction metabolites is thought to be mostly due to endosymbiotic bacteria21. Insect P450s are thought to produce oxidative metabolites exclusively, and the metabolites formed by CYP6G1 have been the best characterized. Heterologous expression of Cyp6g1 in the tobacco plant Nicotiana tabacum produced the metabolites IMI-5-OH, IMI-Ole and IMI-diol22. These results were replicated when driving the expression of Cyp6g1 in D. melanogaster 23. The potential for the other 87 D. melanogaster P450s to be involved in imidacloprid resistance has yet to be tested; Cyp6g1 is so far the only P450 linked to imidacloprid resistance in this species. While it is possible that no other D. melanogaster P450 is capable of metabolizing imidacloprid, this appears unlikely given that many P450s are polyspecific, and multiple P450s are often upregulated in field resistant insects24, 25. Furthermore, the transcriptional response to xenobiotics is often regulated by transcription factors, such as Cap ‘n’ Collar and DHR96, that regulate the expression of many P450s26, 27.

A common method of describing the genetic basis of a trait is to associate genetic variation with phenotypic variation, attempting to identify the causative quantitative trait loci (QTLs) and transcripts. This approach assumes no a priori knowledge about the genes that influence a phenotype, and applying these techniques in model organisms with well characterized genetic resources further enhances detection power. The Drosophila Genetic Reference Panel (DGRP) exemplifies these capabilities. The DGRP is a collection of 201 fully sequenced inbred Drosophila stocks, which represents a snapshot of genetic diversity present in a single population from Raleigh, North Carolina, sampled in 201228, 29. Using the DGRP, a genome wide association study (GWAS) can be performed by testing the associations of the ~2.5 million genetic variants (most commonly single nucleotide polymorphisms; SNPs) across the DGRP genomes with phenotypic data for any quantitative trait. Further, sequencing of the DGRP male and female transcriptomes allowed for similar association studies to be performed with transcript expression level in a transcriptome wide association study (TWAS)30, 31. The DGRP has been used to understand the genetic basis of a wide variety of traits, including insecticide resistance17.

Here, the genetic basis of imidacloprid resistance in the DGRP is described. The Wiggle Index (WI) bioassay32 was used to estimate imidacloprid resistance by measuring acute imidacloprid response, and substantial variation was observed among DGRP genotypes. Quantification of imidacloprid and its metabolites via high-performance liquid chromatography coupled with mass spectrometry (HPLC-MS) showed differences between resistant and susceptible subsets of the DGRP, suggesting that differences in overall imidacloprid metabolism significantly contribute to the differences in resistance. Many QTLs and transcripts were associated with imidacloprid resistance, implicating several genes involved in Central Nervous System (CNS) development as well as the P450s Cyp6g1 and Cyp6g2. The subsequent deletion of Cyp6g1 from two laboratory strains showed no significant differences in imidacloprid resistance, while the same deletion from a resistant DGRP genotype significantly decreased resistance. These deletions allowed for the direct measurement of the contribution of different Cyp6g1 haplotypes to imidacloprid resistance. Cyp6g2 was also linked to imidacloprid resistance in the DGRP by observing increased expression of the gene in the metabolically relevant tissues (midgut and Malpighian tubules) in resistant genotypes. Transgenic overexpression confirmed its ability to metabolize and confer resistance to imidacloprid. These data suggest that genetic variation in CNS development and the expression levels of Cyp6g1 and Cyp6g2 contribute to imidacloprid resistance in field populations of D. melanogaster.

Results

Analysis of the Distribution of Imidacloprid Resistance

Measurement of imidacloprid resistance in the DGRP was accomplished using the WI32, which measures the acute (one hour) motility response of third instar larvae to insecticide exposure, at two doses (25 and 100 ppm). Imidacloprid resistance was quantified by relative movement ratios (RMRs), which reflect the motility of larvae after one hour of insecticide exposure relative to the motility of the same larvae before exposure (An RMR of 1 reflects no imidacloprid response, while an RMR of 0 implies the strongest possible response). Substantial variation in mean RMRs was observed among the 171 DGRP genotypes tested (Fig. 1; Supplementary Table S1; Supplementary Table S2). The mean RMRs of each genotype in the population showed significant correlation between the two doses (Adjusted R2 = 0.18; p-value <3.6 × 10−9, Supplementary Fig. S1), but not between the amount of initial motility and final RMR (Adjusted R2 = 0.001, p-value = 0.06; Supplementary Fig. S2). This suggests that the imidacloprid response was independent of larval motility measured in the absence of imidacloprid. The 25 ppm exposure produced a slightly left skewed distribution of RMRs and is discontinuous due to 3 extremely susceptible genotypes. The 100 ppm dose produced a more even distribution of RMRs. Broad sense heritability estimated for each dose was estimated to be H2 25ppm = 0.628, H2 100ppm = 0.699.

Figure 1
figure 1

Phenotypic Spread of the DGRP. Imidacloprid response in the DGRP was assessed using the Wiggle Index at 2 doses (A) 25 ppm and (B) 100 ppm. Relative movement ratios represent the amount of imidacloprid response with a value of 1 reflecting no response and 0 the most substantial response. At each dose there was a spectrum of phenotypic responses ranging from susceptible (blue) to resistant (red). Error bars represent the standard error of the mean.

Imidacloprid and Metabolite Quantification in the DGRP

To estimate the contribution of insecticide metabolism to the observed differences in imidacloprid response, the amount of imidacloprid and its metabolites (IMI-5-OH and IMI-Olefin) recovered from both larvae and the exposure media was quantified for resistant and susceptible subsets of the DGRP. This was performed under exposure conditions almost identical to those used to assess resistance in the DGRP21, and metabolic phenotypes were tested for correlation with the imidacloprid response measured at 25 ppm (Fig. 1A). The quantity of imidacloprid in larval bodies showed a significant positive correlation with RMR at 25 ppm among genotypes. The more imidacloprid found in the body of a genotype, the stronger the imidacloprid response (Fig. 2A). Additionally, the quantities of IMI-5-OH and IMI-Olefin recovered from the media showed significant negative correlation with RMR, suggesting that increased levels of metabolites in the media provided for a weaker response to imidacloprid (Fig. 2E,F). However, RMR did not significantly associate with the level of either metabolite in the body (Fig. 2B,C). These data suggest that imidacloprid metabolism is higher in resistant genotypes and that these metabolites are preferentially excreted from the body.

Figure 2
figure 2

HPLC-MS in Resistant and Susceptible DGRP Subsets. The amount of Imidacloprid (A,D), IMI-5-OH (B,E) and IMI-Olefin (C,F) recovered from larval bodies (AC) or the media (DF) is reported in parts per billion (ppb) from a subset of the most susceptible (blue) and resistant (red) DGRP lines. No data is presented for imidacloprid in the media, due to the relative abundance of this molecule in the media, which makes detecting changes impossible. Error bars represent standard error of the mean. Stars represent the significance of association of each phenotype with RMR among genotypes using Pearson’s correlation test (*p < 0.05; **p < 0.001; ***p < 0.0001).

A GWAS for Imidacloprid Response Yields Many Neuronal Candidate Genes

A GWAS was performed in order to identify the genetic basis of imidacloprid resistance. The genome wide association of the scores (−log p-values) of annotated genetic variation in the DGRP was uniformly distributed, suggesting test-statistics were not inflated (Supplementary Fig. S3). Manhattan plots showed only 30 variants that crossed the P ≤ 10−5 threshold and only one that fell below the Bonferroni threshold (Supplementary Fig. S4). Linkage disequilibrium between associated variants was low; only two minor linkage disequilibrium peaks were found among the associated variants.

The annotated function of the genes nearest to the significantly (P ≤ 10−5) associated genomic variants implicated a high proportion of candidates having roles in the development and function of the CNS (Supplementary Table S3). 52.6% (10/19) of these genes are reported to have enriched expression in the third instar CNS (expressed 2 fold or greater compared to all other third instar tissues) compared to 19.6% when all D. melanogaster genes are considered. A number of the candidate genes have not been annotated, precluding Gene Ontology term analysis, but several genes appear to have described roles in CNS development or function. A single SNP near the Sickie gene was the only variant which elicited a p-value below the Bonferroni threshold (p = 0.04) at 25 ppm and was also the only variant to be significantly associated at both 25 ppm and 100 ppm (Supplementary Table S3). Although this finding suggests a potential role for Sickie in imidacloprid response, the associated SNP is intronic and Sickie was not transcriptionally associated with imidacloprid response, precluding a firm understanding of how this variation influences imidacloprid response. However, these data collectively suggest that CNS development and function may contribute to imidacloprid response.

A Transcriptome Wide Association Study Suggests Cyp6g1 and Cyp6g2 are Influence Imidacloprid Response

RNA-seq data for 185 DGRP genotypes31 was used to associate the expression level of specific genes with imidacloprid response. Unlike the GWAS candidate list, the TWAS candidate list was not enriched for genes expressed in any particular third instar tissue (Supplementary Table S4). Furthermore, no pattern emerged with regard to any process or function. However, the well known DME Cyp6g1 was the most significant candidate at both doses (Fig. 3A, Supplementary Table S4). Cyp6g2 expression was also significantly associated at 100 ppm (Fig. 3B), and the expression of the two genes is highly correlated. There was no evidence that any variant from the GWAS was influencing the expression of any significantly associated transcripts, as no transcript expression QTLs (eQTLs) were present among the significantly associated GWAS variants.

Figure 3
figure 3

Transcriptional Association of Cyp6g1 and Cyp6g2. The association of the (A) Cyp6g1 and (B) Cyp6g2 transcripts with imidacloprid resistance is shown. Each plot compares the transcript’s expression reported by Huang et al.31 to the RMR of each genotype reported in the current study. Points are labelled according to their haplotype at the Cyp6g1 locus. A linear model is fit with a 95% confidence interval.

The Knockout of Cyp6g1 Displays Haplotype Dependent Effects

The imidacloprid response of three Cyp6g1 knockouts was compared to their matched controls in WI bioassay. RAL_517, a BA haplotype, showed significantly less imidacloprid response than RAL_517-Cyp6g1KO at both 25 and 100 ppm (Fig. 4A,B). These findings were not replicated when testing knockouts in the Canton-S and Wxac backgrounds, which carried M haplotypes. No significant differences were found between these knockouts and controls when exposed to 5 ppm imidacloprid, a dose used to detect potential response differences in the far more susceptible Canton-S and Wxac genotypes (Fig. 4C,D). These data suggest that Cyp6g1 makes a significant contribution to imidacloprid metabolism in BA haplotypes but not in backgrounds carrying an M haplotype.

Figure 4
figure 4

Imidacloprid Resistance in Cyp6g1KO Backgrounds. The effect of the removal of Cyp6g1 from 3 backgrounds on imidacloprid response is shown. Control lines are shown as dark colours while lighter colours represent knockouts. Removing Cyp6g1 from RAL_517 increased imidacloprid susceptibility at both (A) 100 ppm and (B) 25 ppm. No significant changes were observed when Cyp6g1 was removed from either Canton-S or Wxac using a discriminatory dose of 5 ppm. Error bars represent the standard error of the mean. Stars represent the significance of the difference between the two genotypes measured by the Students T-test corrected for multiple comparisons with the Bonferroni method (*p < 0.05; **p < 0.001; ***p < 0.0001).

Cyp6g2 Expression is Enriched in the Digestive Tissues of Resistant Genotypes and Metabolizes Imidacloprid

The potential for the other P450 genes adjacent to Cyp6g1 (Cyp6g2 and Cyp6t3) to contribute to imidacloprid resistance in the DGRP was tested, by quantifying the expression of all three genes in the digestive tissues (midgut and Malpighian tubules) in a subset of DGRP genotypes (two AA and two M haplotypes). All samples showed consistent expression of the housekeeper gene RP49 with the exception of the 3rd and 4th biological replicates of the RAL_360 genotype; these samples were excluded from the analysis. The remaining genotypes showed significantly higher levels of both Cyp6g1 and Cyp6g2 in the midguts and Malpighian tubules of AA haplotypes compared to M haplotypes (Fig. 5A–C; Supplementary Table S5). Although the upregulation of Cyp6g2 appears far weaker than for Cyp6g1, both genes showed larger differences between haplotypes when only the digestive tissues were considered compared to the whole body data reported previously (Fig. 5D–F, Supplementary Table S5)31. No significant patterns were observed for Cyp6t3. These data suggest that both Cyp6g1 and Cyp6g2 are upregulated in these digestive tissues in the AA haplotype, compared to the ancestral M haplotype.

Figure 5
figure 5

The Expression of P450s of the Cyp6g1 Locus in the Digestive Tissues. The expression of (A) Cyp6g1 (B) Cyp6g2 and (C) Cyp6t3 was quantified in the digestive tissues of third instar larvae in 4 DGRP genotypes. M haplotypes are represented by blue points and AA haplotypes by red. The whole body expression data presented by Huang et al.31 is presented for the same genes (DF). Error bars reflect the standard error of the mean. Significant differences were detected between all Cyp6g1 and Cyp6g2 measurements between M and AA haplotypes (ANOVA Tukey’s honestly significant difference) with p-values reported in Supplementary Table S5.

The relative capacity of Cyp6g1 and Cyp6g2 to confer resistance to imidacloprid was tested by overexpressing the genes in the digestive tissues using the HR_GAL4 driver14 and two newly created UAS genotypes, which contained each gene’s open reading frame in a common insertion site on the second chromosome. Compared to their controls, genotypes overexpressing either Cyp6g1 or Cyp6g2 showed significantly higher resistance to imidacloprid, with the magnitude of resistance conferred by Cyp6g1 being significantly higher (Fig. 6A,B; Supplementary Table S6). Although mRNA levels were not measured, the increased resistance conferred by Cyp6g1 relative to Cyp6g2 expressed from the same insertion site suggests that Cyp6g1 enzyme may have a higher capacity to confer resistance to imidacloprid. To further verify the capacity of Cyp6g2 to metabolize imidacloprid, HPLC-MS was employed to measure levels of imidacloprid and metabolites in a previously reported Cyp6g2 overexpression genotype33. While the HR_GAL4 x w1118 control produced relatively low levels of metabolites and had high levels of imidacloprid in the body, HR_GAL4 x UAS-Cyp6g2-3d larvae produced higher levels of imidacloprid metabolites and had less imidacloprid in the body (Supplementary Fig. S5A,E,F). These data indicate that Cyp6g2 can act as a DME against imidacloprid, although it is possible that its capacity to confer resistance may be less than that of Cyp6g1.

Figure 6
figure 6

The Overexpression of Cyp6g1 and Cyp6g2 Confers Imidacloprid Resistance. The HR_GAL4 driver was used to overexpress Cyp6g1 (magenta) and Cyp6g2 (orange) from a common insertion site and their imidacloprid resistance was compared to their background control (grey). The Wiggle Index bioassay measured imidacloprid resistance and was performed at (A) 20 ppm and (B) 40 ppm in order to assess the magnitude of each gene’s ability to confer resistance. Error bars reflect the standard error of the mean. Significant differences were detected between all genotypes at all doses (ANOVA Tukey’s honestly significant difference) with p-values reported in Supplementary Table S6.

Discussion

Of the candidates that were significantly (p < 10−5) associated with imidacloprid response in the GWAS, several were in or near genes that have annotated roles in CNS development or function (Supplementary Table S3). Representative of this group is Sickie, the only candidate gene to be associated below the Bonferroni threshold (p < 2.65 × 10−8). Sickie is orthologous to mammalian NAV2, which has been shown to regulate neuronal development34. Although originally identified as a regulator of Relish and posited to have a role in innate immune response35, Sickie has also been implicated in mushroom body development in D. melanogaster. Its expression is also highly enriched in the third instar larvae CNS, and genotypes carrying hypomorphic Sickie alleles showed axon growth defects36. The variant within Sickie reported in the current study is intronic, and therefore will not change the amino acid sequence. Intronic polymorphisms can regulate a gene’s expression, but Sickie expression was not significantly associated with imidacloprid response. It is possible that the the SNP found in our study is not itself causative, but rather linked to a causative variant nearby. Thus, the data presented here suggest variation within Sickie contributes to imidacloprid response, but the nature of this contribution will need to be addressed by future studies.

The biological mechanisms by which any variant of Sickie and other neuronal GWAS candidates might influence imidacloprid response are also unknown. It may be that variation in such genes changes the amount of nAChRs at the synaptic membrane that can be targeted by imidacloprid. Alternatively, changing the connectivity of neurons could influence how imidacloprid’s signal is propagated after the insecticide has bound its target. Transgenic techniques are often used to explore the biological function of genes implicated in GWAS, but this was not performed in the current study. Many of these genes are homozygous lethal if knocked out. Furthermore, it is not clear that gene knockout, knockdown or overexpression would reproduce a resistance phenotype because these alleles would differ from those observed in the DGRP and may carry severe fitness costs. If the specific SNP associated with resistance was introduced via CRISPR, the difference may be too small to detect. However, the enrichment for neuronal genes in the GWAS candidate list suggests a role for CNS function in imidacloprid response. This may be a mechanism common to neurotoxic insecticides as a GWAS for azinphos-methyl resistance in the DGRP also implicated a high proportion of candidates enriched in the CNS (37.1% compared to 19.6% genome wide)17. The higher neuronal proportion of CNS enriched candidates (52.6%) in the current study may be due to the WI measuring a behavioural phenotype, that may be more strongly influenced by CNS function. It should be noted that variation in genes encoding the nAChR subunits known to be targeted by imidacloprid, Dα1 and Dβ2 5, was not shown to impact the insecticide response. This may be due to an absence of suitable variation in the DGRP. Studies with laboratory mutants suggest that resistance via these genes is associated with a loss of function that may involve a fitness cost6.

The capacity to metabolize imidacloprid also contributes to resistance in the DGRP; consideration of a subset of 9 of the most resistant and susceptible DGRP genotypes revealed that more imidacloprid metabolism occurred in resistant genotypes compared to susceptible ones (Fig. 2). However, such differences were only apparent in the media, unlike the RAL_517-Cyp6g1KO, which showed significantly different metabolite levels in the media and body at both one and six hour time points21. This suggests that excretion is playing a critical role in the imidacloprid response within the DGRP. There appears to be a genetic component to excretion as the ratio between body and media metabolite levels varied between genotypes. In particular, RAL_509, the most susceptible genotype in this study, displays an interesting set of phenotypes. This genotype has a BA haplotype at the Cyp6g1 locus, and produces high levels of both IMI-Olefin and IMI-5-OH levels, but these metabolites are disproportionately retained in the body (Fig. 2). Although the genetic basis for excretion was not explored in the current work, transporter proteins, such as ATP binding cassette transporters, have been implicated in insecticide resistance recently and could underpin differences in excretion between DGRP genotypes37, 38. Further manipulation of transporters could reveal insights into how insecticides are excreted from the body.

The role played by metabolism in imidacloprid resistance was reinforced by the implication of Cyp6g1 and Cyp6g2 in the TWAS (Fig. 3; Supplementary Table S4). Although the involvement of Cyp6g1 in imidacloprid resistance was known previously13, the magnitude of the contribution of Cyp6g1 had not been tested. Removal of Cyp6g1 from M haplotypes did not affect WI response, while knockouts in a BA haplotype reduced WI response and imidacloprid metabolism21 (Fig. 4). It may be that in M haplotypes there is not sufficient Cyp6g1 to contribute to imidacloprid resistance. Previous studies have used RNAi to knock down the expression of Cyp6g1 in wild-type backgrounds and seen either a slight increase in susceptibility to DDT or no change39, 40. While lacking a clear understanding about the correlation between imidacloprid response and the amount of Cyp6g1 expression, our data show that Cyp6g1 does not contribute significantly to imidacloprid response in M haplotypes as measured with the WI. This does not rule out the possibility that such differences may be detected with other toxicological bioassays or with other insecticides.

The role of P450s apart from Cyp6g1 in imidacloprid resistance was unknown in this species prior this study. Detection of the same metabolites produced by CYP6G1 (IMI-5-OH and IMI-Olefin) in RAL_517-Cyp6g1KO suggested that other P450s metabolize imidacloprid in that background21. The ability of Cyp6g2 to both metabolize and confer resistance to imidacloprid when transgenically expressed suggests that this gene may be one source of the residual resistance in RAL_517-Cyp6g1KO (Fig. 6; Supplementary Fig. S5). Based on evidence from the literature it is widely believed that there are two groups of P450s, those involved in metabolism and those involved in development. The Cyp6g2 gene falls into both of these groups. In laboratory strains Cyp6g2 is specifically expressed in the corpus allatum within the ring gland41 and is implicated in the synthesis of juvenile hormone42. However, ectopic expression of this gene in digestive tissues showed that it was able to confer resistance to nitenpyram and diazinon33. The current work extended the substrate specificity of CYP6G2 to imidacloprid, finding that it was able to produce the same metabolites as CYP6G1 and confer resistance (Fig. 6; Supplementary Fig. S5).

Although Cyp6g2 expression was significantly associated with imidacloprid resistance in the TWAS, it was not known whether upregulation of this gene in resistant genotypes is restricted to its native ring gland specific expression pattern. Significantly higher levels of Cyp6g2 expression in the digestive tissues of AA genotypes (Fig. 5E) suggests that this gene may contribute meaningfully to imidacloprid metabolism within the DGRP. This increase is at least partially tissue specific, as increases in the digestive tissues were far higher than those in whole adult bodies (Fig. 5B,E). The most parsimonious source of the change in expression level and pattern is the presence of an Accord element, which could act as an enhancer, increasing the expression of both Cyp6g1 and Cyp6g2 in the digestive tissues. Cyp6g2 may represent the limit of Accord’s range as expression of the more distant Cyp6t3 did not appear to be influenced by the presence of Accord. The regulation of Cyp6g2 by Accord is not the only possible mechanism for of the observed expression in metabolic tissue. Other differences between AA and M haplotypes could influence Cyp6g2 expression and the expression pattern in BA haplotypes is unknown. Small sample sizes and the consideration of only AA and M haplotypes preclude the establishment of any definitive mechanism for regulating Cyp6g2 expression. Unfortunately, it is unlikely that the contribution of Cyp6g2 to imidacloprid resistance can be measured directly as was done for Cyp6g1. Previous work has shown RNAi of Cyp6g2 was lethal41 meaning that knockout experiments in RAL517 would likely be unsuccessful.

Much still remains unresolved about the relative contribution of different alleles to complex phenotypes such as insecticide resistance. While far more attention has been given to cases of monogenic resistance, all populations reflect a distribution of resistance levels among individuals that is governed by many loci2. This is true even in populations where resistance alleles have gone to fixation. Variation at the Cyp6g1 locus contributes significantly to imidacloprid resistance in the DGRP and is likely the largest single factor in determining the likelihood an insect survives an exposure. However, while removal of Cyp6g1 from the resistant RAL_517 genotype increased susceptibility to imidacloprid, the RAL_517-Cyp6g1KO genotype was still more resistant than approximately half of the DGRP genotypes (Fig. 4A). This indicates that Cyp6g1 significantly contributes to imidacloprid resistance, but also highlights the supporting role played by other genes such as Cyp6g2 and Sickie which likely contribute in smaller ways. Imidacloprid resistance in the DGRP can then be thought of as polygenic but with a single gene making a contribution far larger than the rest.

Methods

Fly Genotypes

All genotypes used in this study were ordered from the Bloomington Drosophila Stock Center (Bloomington, Indiana) or generated in this study (Supplementary Table S7). 178 of the 201 DGRP genotypes were tested in the the initial WI screen. A subset of 9 of the most imidacloprid resistant and susceptible DGRP genotypes were chosen for imidacloprid metabolism analysis via HPLC-MS. The Actin-Cas9 genotype (Bloomington #54590) and a genotype containing an attP landing site (attP40; 2L:5,108,448:5,108,448; Bloomington #25709) were used for used to knock out the Cyp6g1 gene in three different genetic backgrounds. Canton-S-Cyp6g1KO was generated in the Canton-S background, Wxac-Cyp61KO was generated from the Wxac background (Actin-Cas9 with the X chromosome replaced with the one from the w1118 genotype). RAL_517-Cyp6g1KO was created in the RAL_517 background, a genotype chosen due to its high level of imidacloprid resistance and Cyp6g1 expression among DGRP genotypes. RAL_517 carries a BA haplotype while Canton-S and Wxac both carry M haplotypes and were far more susceptible to imidacloprid.

The three Cyp6g1 knockouts generated here were created by using a transgenic CRISPR strategy described recently43. Briefly, Cyp6g1-sgRNA plasmids were made by first cutting the PCFD4 plasmid (Addgene #49411) with the restriction enzyme Bsb1. A separate fragment was generated by amplifying a portion of the PCFD4 plasmid with the Cyp6g1-PCFD4 primer set (Supplementary Table S8; 60 °C annealing, 1 minute extension), introducing Cyp6g1 sgRNAs into the PCR product. The cut plasmid and PCR product were then reassembled using the Gibson assembly kit (New England Biolabs) to make a circular PCFD4 plasmid with two Cyp6g1 sgRNAs under the control of two U6 promoters. Verification of this modification was accomplished by sequencing the plasmid using the PCFD4_seq primer (Supplementary Table S8).

This plasmid was then injected into a genotype expressing φ-31 integrase and which contained an attP landing site (attP40; 2L:5,108,448:5,108,448; Bloomington #25709), both of which facilitated the integration of the modified PCFD4 into the germline. Transgenic flies were identified by scoring the visible marker vermilion eyes which was restored to wild type upon successful PCFD4 integration. Chromosomes from Actin-Cas9 and the sgRNA expressing genotype were brought together in a crossing scheme which made near identical deletions of Cyp6g1 in each background (Supplementary Fig. S6A). The resulting deletion was confirmed by amplifying across the Cyp6g1 deletion using the Cyp6g1-KO primer set (Supplementary Table S8; 56 °C, 2 minutes extension). Sequencing this PCR product revealed the almost complete removal of the gene (Supplementary Fig. S6B), and the remaining transcript was predicted at 73 amino acids long (reduced from 524). Failure to amplify any full copies of the gene indicated that both copies of Cyp6g1 present in this RAL_517 were successfully deleted.

Overexpression of the genes Cyp6g1 and Cyp6g2 was achieved in the midgut Malpighian tubules and the fat body using the GAL4/UAS system44 and the HR_GAL4 driver14. In testing the capacity of CYP6G2 to metabolize imidacloprid, a previously reported UAS-Cyp6g2-3d genotype was used33. However, this genotype was created using a random insertion method, precluding a direct comparison of the impact of Cyp6g1 and Cyp6g2 on resistance. Hence, new UAS-Cyp6g1 and UAS-Cyp6g2 genotypes that share a common insertion site (attP40) were generated in this study. Open reading frames from each gene were amplified from the w 1118 background using Q5 high fidelity master mix (New England Biolabs). Each PCR fragment was A-tailed and cloned into the PGEM T-easy vector (Promega). The fragments were then cut out of this vector using the NotI enzyme and ligated into the pUASt-attB vector45. Plasmids were injected into a genotype carrying the same attP landing site genotype used for knockouts (attP40) and an X chromosome with a white- allele (Bloomington stock # 24749). These two genotypes were compared to their control, which was created by injecting an empty pUASt-attB vector into the same background.

The Wiggle Index

The response of D. melanogaster larvae to imidacloprid was measured using the WI assay, which estimates the insecticidal effect by quantifying the reduction in motility during insecticide exposure32. Third instar larvae of each genotype were picked, 25 per well, into a NUNC cell culture treated 24 well plate (Thermo-Scientific) preloaded with 200 µL 5% w/v sucrose (AR Grade; Chem Supply) in distilled H20. Larvae were filmed for 30 seconds at two time points, 0 min (before starting the exposure) and at 60 min after the addition of specific concentrations of imidacloprid (200 g L−1 Confidor®; Bayer Crop Science). Subsequently, the WI ImageJ script was used to quantify the total motion in each well at each time point. The ratio between initial and final motility was used to calculate Relative Movement Ratios (RMRs), which were averaged to estimate imidacloprid response for each genotype tested in this study. 178 DGRP genotypes were considered at doses of 25 and 100 ppm as was RAL_517 and RAL_517-Cyp6g1KO. Other Cyp6g1 knockouts were tested at 5 ppm. UAS-Cyp6g1 and UAS-Cyp6g2 were tested at 20 and 40 ppm.

Analysis of the WI and Candidate Gene Selection

Broad sense heritability (H2) of imidacloprid resistance was measured by comparing the variance of RMRs within genotypes to the variance between genotypes. The association of initial motility values (at 0 minutes) with final RMR was also tested in order to test the any confounding effects of starting movement on imidacloprid response. The correlation between 25 and 100 ppm RMRs was tested to observe the relationship between the two phenotypes. All these associations were tested by assessing the fit to a linear model using Pearson’s correlation test.

In order to implicate individual QTLs in imidacloprid resistance, mean RMRs for each genotype were used as input data for the Mackay lab DGRP GWAS pipeline29 found at (http://dgrp2.gnets.ncsu.edu/). Responses to the two doses (25 and 100 ppm) were considered as separate phenotypes. Similarly, a transcriptome wide association study, was performed by associating the phenotypic Wiggle data with expression data recently generated in 185 DGRP genotypes31. Analysis was performed using a modification of a recently reported pipeline17, which tested the fit of linear models (Pearson’s correlation test) correlating the expression level of an individual transcript and mean RMR among genotypes. Transcript levels for each gene were averaged between males and females as larvae were not sexed before testing in the WI assay. Furthermore, the expression quantitative trait loci (eQTLs) of each significant transcript31 were cross examined with significant SNPs from the GWAS in order to test if significantly associated genomic loci regulate the expression of any significantly associated transcripts.

Candidate genes from the GWAS and TWAS were chosen based on their significance of association with either the 25 or the 100 ppm phenotype. Any GWAS variant eliciting a p-value of less than 10−5 and any transcript with a p-value of less than 10−3 was considered as a candidate in the current study. In the case of the GWAS, variants were assigned to their nearest annotated gene according to the software SnpEff (built into the Mackay pipeline), and only p-values derived from simple regression were considered (mixed models were not considered). Candidate genes were further examined by testing enrichment in third instar larval tissues using the modEncode transcriptome datasets46.

Quantification of Imidacloprid and its Metabolites

Levels of imidacloprid and its metabolites were quantified via HPLC-MS following the method described by Fusetto et al.21. Quadruplicates of 200 third instar larvae from a given genotype were placed into 200 µL of 5% analytical reagent sucrose and were exposed to a 50:50 mix of 12C6:13C6 imidacloprid (99% analytical reagent) at a concentration of 25 ppm for 1 hour. Larvae were then recovered from the solution and washed 3 times with 3 mL of distilled H2O to remove chemical from the cuticle, and the exposure media was collected separately in order to estimate the amount of each chemical excreted. The compounds in each sample were quantified using high-performance liquid chromatography (HPLC) coupled with an Agilent 6520 Q-TOF Mass Spectrometer (Agilent Technologies, Inc., Santa Clara, CA, USA). These measurements were taken for a subset of the 9 most resistant and susceptible DGRP lines and the UAS-Cyp6g2-3d genotype33 as well as all relevant control lines. In the case of the DGRP subset, the correlation between the mean of each metabolic phenotype (each chemical in larval bodies or excreta) and the mean of each genotype’s WI RMR at 25 ppm (Fig. 1A) was tested among genotypes using Pearson’s correlation test. All other comparisons were made using a students t-test to compare genotypes, with a Bonferroni correction applied to correct for multiple testing. Values was reported in either parts per billion (ppb) or, in the case of UAS-Cyp6g2-3d, the area under the chromatogram peak was used as a relative quantity of the chemical.

Measurement of Cyp6g1 and Cyp6g2 Expression in Digestive Tissues

The expression of Cyp6g1, Cyp6g2 and Cyp6t3 was measured in digestive tissues (midgut and Malpighian tubules) of a subset of DGRP lines. 4 lines were chosen in total. Two carried the AA haplotype (RAL_138, RAL_360) and two carried the M (lines RAL_776, RAL_843). 20 midguts and the associated malpighian tubules from each genotype were used for a single biological replicate and 5 replicates were taken from each genotype. RNA was extracted from each sample using the TriSure (Bioline) protocol, and cDNA was generated using the Superscript III reverse transcriptase (NEB). All primers and qPCR parameters were reported previously33. Expression was quantified using the 2−ΔΔ CT method and normalized to the level of the genotype showing the lowest expression level tested (RAL_776). Differences in expression among genotypes were compared using analysis of variance (ANOVA) with Tukey’s honestly significant difference post-hoc test.

Data Analysis and Availability

All raw data generated in this study is available for download at https://github.com/shanedenecke/Denecke_et_al_2017_Imidacloprid_DGRP and was analysed in R47. An accompanying R markdown document, also available on Rpubs (http://rpubs.com/shanedenecke/Denecke_DGRP). This document provides R code to obtain the raw data and reproduce all the figures and tables presented in this study. Supplementary Figs S7 and S8 are not included in the markdown document because these do not contain analysis driven information. Unless otherwise stated all p values represent Student’s t-tests (for pairwise comparisons), ANOVA with Tukey’s post-hoc test (3 comparisons of three or more variables) and Pearson’s correlation test (for correlations and linear regressions).