Cis- and trans-acting variants contribute to survivorship in a naïve Drosophila melanogaster population exposed to ryanoid insecticides.

Insecticide resistance is a paradigm of microevolution, and insecticides are responsible for the strongest cases of recent selection in the genome of Drosophila melanogaster Here we use a naïve population and a novel insecticide class to examine the ab initio genetic architecture of a potential selective response. Genome-wide association studies (GWAS) of chlorantraniliprole susceptibility reveal variation in a gene of major effect, Stretchin Myosin light chain kinase (Strn-Mlck), which we validate with linkage mapping and transgenic manipulation of gene expression. We propose that allelic variation in Strn-Mlck alters sensitivity to the calcium depletion attributable to chlorantraniliprole's mode of action. GWAS also reveal a network of genes involved in neuromuscular biology. In contrast, phenotype to transcriptome associations identify differences in constitutive levels of multiple transcripts regulated by cnc, the homolog of mammalian Nrf2. This suggests that genetic variation acts in trans to regulate multiple metabolic enzymes in this pathway. The most outstanding association is with the transcription level of Cyp12d1 which is also affected in cis by copy number variation. Transgenic overexpression of Cyp12d1 reduces susceptibility to both chlorantraniliprole and the closely related insecticide cyantraniliprole. This systems genetics study reveals multiple allelic variants segregating at intermediate frequency in a population that is completely naïve to this new insecticide chemistry and it foreshadows a selective response among natural populations to these chemicals.

An elaboration of the adage of Paracelsus (1493–1541) that “the dose makes the poison” is that there is a dose range of insecticides that kills some but not all insects in a population. By examining the genetic variation that contributes to survivorship on such discriminating doses, we can take a genetics approach to address a diverse set of questions relating to insecticide biology. Which genes have variants that affect survivorship, and how do they combine to provide the genetic architecture underpinning the trait? Do they provide insights into the mode of action of new insecticides? Do they suggest likely mechanisms by which insecticide resistance will arise? And what else do they tell us about past and future evolutionary responses to insecticides in pest and nontarget species?

Chlorantraniliprole (Rynaxapyr), is the first of the anthranilic diamides, a new class of insecticides. Unlike earlier insecticides that predominantly target neurotransmission, the anthranilic diamides are designed to target the ryanodine receptor, which is primarily involved in calcium homeostasis and muscle contraction (1, 2). Disruption of ryanodine receptor activity causes rapid incapacitation of the pest, leading to feeding cessation, lethargy, paralysis, and death (3, 4). Therefore, both the mode of action and the chemistry suggest that cross-resistance with older insecticides is unlikely.

Chlorantraniliprole was first sold in the Philippines in 2007, and worldwide soon after (5). Within years of introduction, resistance cases were reported in the diamondback moth Plutella xylostella (6, 7), and the tomato leafminer Tuta absoluta (8). While some of these cases can be attributed to mutations in the ryanodine receptor, the primary molecular target of these insecticides, there are others that suggest that resistance to this new insecticide class can arise through other means (9–12).

While Drosophila melanogaster is not a pest or a direct target of chlorantraniliprole applications, it is an organism of interest for two reasons. Firstly, D. melanogaster has long served as a model for insecticide resistance (13) and its status as a model organism more generally means that there is a wide variety of tools available to characterize genetic traits (14). Secondly, selective sweep analyses show that insecticides (particularly the organophosphates) have been major selective agents on D. melanogaster populations (15–18). These findings support the proposition that D. melanogaster can be used as a sentinel species for environmental pollutants, particularly insecticides (19, 20).

Like pest insects, D. melanogaster evolves insecticide resistance chiefly through target molecule insensitivity or detoxification enzyme adaptation, although other resistance mechanisms have been characterized (21). Resistance mutations in the target site genes typically diminish insecticide binding (e.g., ref. 22). Resistance mutations affecting detoxification enzymes can alter the protein sequence (e.g., ref. 23), but more generally increase transcriptional output through copy number variation (CNV) or cis-regulatory changes in the promoters of the resistance genes (e.g., refs. 24 and 25). Master regulatory genes that control, in trans, detoxification pathways have been reported in multiple arthropod species (26–33). Increased constitutive activation of these pathways has been shown to correlate with resistance in pest insects (32, 34) as well as in D. melanogaster (35), but so far natural variation underpinning such phenotypes has not been described at nucleotide resolution.

A powerful addition to the D. melanogaster toolkit is the Drosophila Genetic Reference Panel (DGRP) (36), which comprises 205 inbred lines derived from a single North American population. Each line in the DGRP has been bred for homozygosity and its genome sequenced, creating a “living library,” designed for genome-wide association studies (GWAS) that will associate genetic variants with phenotypes. The DGRP has been phenotyped for an extensive number of traits (37), including various insecticide phenotypes (16–18, 23, 38, 39) and “intermediate phenotypes” such as transcript abundance that enables eQTL to be mapped (40). Thus, the DGRP is becoming an important model for systems genetics of insects (41, 42), which, as we demonstrate here, enables deeper characterization of the genetic and regulatory mechanisms underpinning traits.

The DGRP lines were established from a collection of flies from the Farmers Market in Raleigh, NC in 2003 (36), before chlorantraniliprole became commercially available in 2007. Thus, they are naïve with respect to the completely novel class of chemistry of the group 28 insecticides. This provides a rare opportunity to examine the ab initio state of a potentially adaptive trait at unprecedented genetic resolution.

Here, we investigate the genetic architecture of chlorantraniliprole susceptibility across multiple doses. We interrogate associations between survivorship on food with varying concentrations of chlorantraniliprole and both genomic and transcriptomic variation. We find that an allele of large effect already segregates within this population, which is detectable only at higher levels of exposure. In addition, we show that phenotype to transcriptome associations reveal a completely different set of candidate genes, linked by a common trans-regulatory pathway.

Results

Phenotype to Genome Associations.

A total of 152 DGRP lines were scored for larval survivorship on six concentrations of chlorantraniliprole (Fig. 1 and ). The broad-sense heritability (H2) of the six concentrations ranged from 0.73 to 0.85, indicating a strong genetic component to chlorantraniliprole survivorship across the DGRP. GWAS on each of the concentrations identified 42–335 variants below the arbitrary genome-wide significance threshold (P < 1 × 10−5). The number of variants passing this threshold increased with concentration; however, this relationship failed to account for linkage disequilibrium (LD). Therefore, we examined the explanatory power of phenotype-associated DGRP variants for each chlorantraniliprole concentration using a multivariate genomic prediction model and found that more genes are required to explain the genetic architecture of the lowest concentration (0.5 μg/mL). For example, if the top 50 most-associated variants of the 0.5-μg/mL concentration are considered, they explain about the same amount of the phenotypic variation (R = 0.43) as the top five variants of the 5-μg/mL concentration ().

Fig. 1.

Chlorantraniliprole phenotypes of 152 DGRP lines. (A) The survivorship of each line on the different concentrations (colored). The lines are arranged in ascending order for each concentration. (B) LC50, calculated from a minimum of six single-concentration phenotypes.

These data were supplemented with screening at additional concentrations for some DGRP lines to estimate the concentration of chlorantraniliprole required to kill 50% of individuals in each line (LC50). The mean LC50 was equal to 1.4 μg/mL (SD = 2.07 μg/mL) while the maximum LC50 was 17 μg/mL (Fig. 1). A GWAS of the LC50 phenotype identified 931 associated variants below the arbitrary genome-wide significance threshold (P < 1 × 10−5; Fig. 2). Ninety-six of these remained significant after a Bonferroni correction for multiple testing (2.65 × 10−8; ). The strongest association was to a single nucleotide polymorphism (SNP) in an intron of Stretchin myosin light chain kinase (Strn-Mlck) (2R:11853686; P = 2.03 × 10−17; Fig. 2), and variants annotated to this gene accounted for 15 of the 96 Bonferroni-significant GWAS variants.

Fig. 2.

(A) Top associations with genomic variants. Manhattan plot of DGRP chlorantraniliprole LC50 associations. Strn-Mlck is a standout candidate, with a minimum P value of 2.03 × 10−17. (B) LD at LC50 GWAS top candidates sli and Strn-Mlck.

To account for the fact that the strong effect of Strn-Mlck variation may be influencing other GWAS associations, we fitted the effect of Strn-Mlck to the chlorantraniliprole LC50 data and ran a GWAS on the residuals. The gene with the most highly associated variants after Strn-Mlck in the original LC50 GWAS, sli, was maintained in the Strn-Mlck–corrected GWAS. This demonstrates that despite its proximity to Strn-Mlck, associations with variants in sli are not artifacts of LD with variants in Strn-Mlck (Fig. 2 and ).

sli is like other genes harboring highly associated variants in that it is involved in axon guidance, sarcomere organization, intracellular signaling, and regulation of cell growth. The R spider software (43) was used to identify networks of genes in the Kyoto Encyclopedia of Genes and Genomes (KEGG) or Reactome pathways that were enriched for variants associated with chlorantraniprole phenotypes. The top 100 variants from LC50 genomic prediction both before and after correcting for the Strn-Mlck association were considered. The uncorrected LC50 network contained 12 genes, five with GWAS associations, out of a total of 54 recognized by R spider (Monte Carlo simulation P = 0.01; ); the Strn-Mlck–corrected LC50 network contained 10 genes, five with GWAS associations, out of a total of 53 genes with annotated function, of which seven could be mapped to a reference network previously described in KEGG (P = 0.005; ). Three GWAS-associated genes appeared in both networks: sli, robo3, and norpA. Although not connected to either network based on Reactome or KEGG databases, Strn-Mlck can be linked to the corrected LC50 network, as it has been shown to be one of the serine/threonine kinases that phosphorylates the transcription factor foxo (44).

Phenotype to Transcriptome Associations.

Transcript-level variation from the DGRP transcriptome dataset (mean for males and females) (40) was associated with chlorantraniliprole LC50 for nine genes below Bonferroni significance [generalized linear model (GLM)] (P < 2.76 × 10−6; Fig. 3). These top candidates are enriched for members of genetically correlated transcription modules (40) for both males (module 79; seven transcripts) and females (module 21; six transcripts). Furthermore, eight of the nine Bonferroni-significant transcripts have also been shown to be induced through ectopic expression of the transcription factor Cap‘n’collar (cnc) and following exposure to phenobarbital (Fig. 3) (45). One top candidate, Cyp12d1, exhibits CNV within the DGRP such that two copies (Cyp12d1-p and Cyp12d1-d) are observed in 24% of lines. The DGRP has been explicitly genotyped for duplication of Cyp12d1 (46, 47) and the duplication was found to be correlated with Cyp12d1 transcription levels (Cyp12d1-p: male r2 = 0.12, female r2 = 0.14 and Cyp12d1-d: male r2 = 0.38, female r2 = 0.27). Thus, two distinct contributions to Cyp12d1 transcript level can be identified: that which is attributable to cnc regulation and that which is attributable to CNV at the locus.

Fig. 3.

P values of association between the chlorantraniliprole LC50 phenotype and each DGRP transcript level, grouped into sex-specific genetically correlated transcriptional modules. The transcripts within each module are ordered by rank to give “opera house” plots. Eight of the nine Bonferroni-significant associations (red line; P < 2.76 × 10−6) are affected by ectopic cnc expression.

Validation of Strn-Mlck and Cyp12d1-P in Chlorantraniliprole Survivorship Traits.

To support the involvement of Strn-Mlck in chlorantraniliprole survivorship, two DGRP lines (RAL-59 and RAL-399) that differ by both LC50 and Strn-Mlck variants, were crossed. Comparison of probit curves from parental lines and the F1 suggest that the more resistant allele is recessive (degree of dominance = −0.76 at LC50, −0.84 at LC90). A pair of restriction fragment length polymorphism (RFLP) assays showed that the Strn-Mlck haplotype identified in the GWAS was significantly enriched in chlorantraniliprole-screened F2 survivors over untreated controls (χ2 P = 9.97 × 10−7), indicating that this region is indeed associated with chlorantraniliprole survivorship. The importance of Strn-Mlck was also tested using three separate RNAi lines crossed to the Actin5c driver. In two of the three crosses, flies with Strn-Mlck knocked down showed significantly decreased LC50s relative to CyO siblings (Fig. 4).

Fig. 4.

Transgenic manipulation of candidate genes. LC50 values of combined reciprocal crosses for RNAi knockdown of Strn-Mlck using three TRiP lines HMS02663 (S1), JF02171 (S2), and HMS01665 (S3) crossed to the Actin5c driver. In two out of three of the crosses, Strn-Mlck knockdown results in a significantly decreased chlorantraniliprole LC50 relative to CyO siblings (C1–3). Overexpression of Cyp12d1-p (Y) significantly increases chlorantraniliprole LC50 relative to control (C4). Error bars represent 95% confidence of probit fit at LC50. *P < 0.05, **P < 0.01.

The involvement of Cyp12d1-p was tested using the GAL4-UAS system and the 6g1HR-GAL4 driver (48); flies overexpressing Cyp12d1-p in key metabolic tissues were 2.2-fold more tolerant to chlorantraniliprole than controls (Fig. 4). This system was also employed to test cross-tolerance to the closely related insecticide, cyantraniliprole. At three concentrations, survivorship of flies overexpressing Cyp12d1-p was significantly increased relative to controls (P < 0.05, two-tailed t test assuming unequal variances), suggesting that the enzyme acts on chemical moieties that the two insecticides have in common.

Induction of Cyp12d1.

It has been previously demonstrated that both Cyp12d1-p and Cyp12d1-d are up-regulated in response to phenobarbital and caffeine (45). To test if chlorantraniliprole or cyantraniliprole induces Cyp12d1 expression, Cyp12d1 transcript levels (Cyp12d1-p and Cyp12d1-d were not distinguished) were quantified after larvae from a laboratory strain (Canton-S) and two DGRP lines (chlorantraniliprole-resistant RAL-59 and chlorantraniliprole-susceptible RAL-399) were exposed to chlorantraniliprole, cyantraniliprole, phenobarbital, and caffeine. In accordance with DGRP transcriptome data from Huang et al. (40), we found Cyp12d1 transcript abundance in unexposed larvae significantly higher in the resistant DGRP line RAL-59 relative to susceptible line RAL-399 (Fig. 5). In all three lines Cyp12d1 transcript levels were increased significantly following exposure to phenobarbital and caffeine, but not following exposure to chlorantraniliprole or cyantraniliprole (Fig. 5).

Fig. 5.

Variation and inducibility of Cyp12d1 expression. Cyp12d1 expression (both Cyp12d1-p and Cyp12d1-d detected; normalized to housekeeper CG11322) in third-instar larvae from three DGRP lines unexposed (UN) and after exposure to four xenobiotic compounds; chlorantraniliprole (CH), cyantraniliprole (CY), phenobarbital (PH), and caffeine (CA). Error bars represent SD of the first derivative.

Discussion

Genetic Architecture Changes with Insecticide Concentration.

Theoretical considerations have led to the prediction that the genetic architecture of insecticide resistance would become less polygenic with increasing concentration of an insecticide (49). The DGRP allows identification of polygenes affecting survivorship on different concentrations of an insecticide in an unprecedented way. Concurring with expectations, smaller effect sizes were observed in the lowest concentration GWAS. However, relatively few variants passed the significance threshold in the GWAS at this concentration compared with those performed at higher concentrations. This is a consequence of the increasingly nonnormal distribution of the survivorship among lines as concentration increases and reflects the statistical mechanics and distribution assumptions of GWAS, and the degree of LD around associated variants at lower frequencies. We addressed this by modeling the additive contribution of small-effect alleles using a Bayesian linear regression. At the three highest concentrations, and the concentration inferred to kill 50% of flies (LC50), many of the associated variants are within large stretches of LD around Strn-Mlck. When the LD is accounted for, the number of genes associated with survivorship on high concentrations was indeed lower, consistent with theoretical predictions. The genetic architecture of lower-concentration survivorship is far more polygenic, with very few of the top associations showing strong LD relationships. There is also a loss of explanatory power at 0.5 μg/mL, with R2 plateauing at 69%, making it obvious that the set of associated variants is not capturing all of the variation in this phenotype. In addition, the heritability at the lowest concentration is smaller, suggesting other environmental effects are more pronounced at this concentration, or that epistasis, which is likely to be more prevalent with a greater level of polygenicity, is playing a role (50).

Strn-Mlck: A Novel Gene of Major Effect.

The GWAS presented here identified 15 variants in the Strn-Mlck gene as marking an allele of major effect. This genome-wide and unbiased approach was confirmed by two methods: RNAi knockdown supports the potential of this gene to contribute to the trait, and linkage mapping indicates that a naturally occurring variant of major effect occurs in this vicinity. A member of the Titan family, Strn-Mlck is highly complex, with a total of 33 exons, three start sites, and three poly-A sites, leading to 16 predicted isoforms. Myosin light chain kinases (MLCKs) are serine/threonine kinases whose substrate is the myosin regulatory chain, a component of the thick filaments in muscles. Most MLCKs are calmodulin dependent, meaning they are proteins that are sensitive to calcium homeostasis. Strn-Mlck’s association with the muscle and its sensitivity to calcium, the intracellular release of which is perturbed with ryanodine insecticides, strengthen the case that Strn-Mlck is involved. The alleles that increase survivorship may act downstream of the insecticide binding to the ryanodine receptor, and simply ameliorate the effect of the insecticide. However, there is an alternative hypothesis. Recently several genome-wide analyses of insecticide phenotypes have identified the involvement of genes in neurogenic pathways (38, 51, 52). In this study, we found that GWAS-associated variants are enriched for a network of neurodevelopmental genes which can be linked to Strn-Mlck. Thus, the highly associated variants in Strn-Mlck act by influencing the development of the nervous system such that some flies are less sensitive to the insecticide.

Transcriptomic Associations Identify a Known Detoxification Pathway.

Underscoring the power of the systems approach that is afforded by the DGRP, associations between chlorantraniliprole tolerance phenotypes and gene expression levels (40) implicated a completely different set of genes from GWAS candidates. A strong association was observed between chlorantraniliprole LC50 and the expression of a group of genes known to be coregulated in response to xenobiotics (45). These genes have previously been found to be regulated by cnc (45), a transcription factor well conserved across both vertebrates and invertebrates (53). cnc and its mammalian homolog Nrf2 are CNC-bZIP transcription factors that play a well-conserved role in oxidative stress response (54). While most studies characterize this pathway in terms of its induction in response to oxidative stress, the transcript-level variation identified in this study comes from unexposed flies, suggesting variation in the constitutive expression level of cnc-regulated transcripts among DGRP lines, a phenomenon that has previously been linked to insecticide resistance (34, 35, 54).

The constitutively up-regulated cnc-regulated transcripts identified in the DGRP include potential detoxification enzymes, including two cytochrome P450s (Cyp12d1 and Cyp6a8), a juvenile hormone epoxide hydrolase (Jheh1), two oxidoreductases (CG2065 and CG9360), and an ecdysone kinase-like enzyme (CG6908). Guio et al. (55) demonstrated that the insertion of a Bari1 transposable element upstream of Jheh1 introduces additional cnc binding motifs. This results in increased inducibility of Jheh1 and its tandem paralog, Jheh2, in response to oxidative stress, conferring increased resistance to the organophosphate insecticide malathion.

Cyp12d1 induction has been associated with exposure to multiple xenobiotics, including DDT, caffeine (56), pyrethrum (57), atrazine (58), and piperonyl butoxide (PBO) (60), and its overexpression increases survivorship to DDT, dicyclanil, and malathion (17, 60). Due to its polymorphic status, it has been assumed that the Cyp12d1 duplication is a recent event (61). Duplication frequency in the DGRP is correlated with higher expression levels in adults, a pattern that is also seen in other outbred populations (62). While the DGRP transcriptome data were measured in adults, our qPCR results suggest that this difference is also observed in larvae.

As the transcription output of a gene can be affected by numerous variants, including rare variants, genes that are missed in the GWAS may be identified by phenotype to transcriptome associations. Given that the transcriptome analysis implicated Cyp12d1, we examined its position in all of the datasets more carefully. We identified two cis-eQTL of Cyp12d1 (2R:6994376 and 2R:7007339) (40), associated with survivorship at 4-μg/mL and 1-μg/mL concentrations, respectively. Both these variants are in LD with the duplicated state of Cyp12d1 (R2 = 0.46 and 0.72, respectively), suggesting the effect of the Cyp12d1 duplication was indirectly detected in the GWAS. Huang et al. (40) also report a trans-eQTL for Cyp12d1 occurring in sli, the gene that ranks second in the LC50 GWAS after Strn-Mlck. Thus, trans-regulatory variation may be combining with cis-regulatory variation to affect the transcriptional output of Cyp12d1 and chlorantraniliprole survivorship.

The observation that an eQTL for a cnc-regulated gene maps to a region strongly implicated in the GWAS is intriguing. This study raises the possibility that there is a link between muscle function as implicated by the GWAS and the oxidative stress pathway as implicated by the transcriptome associations. There is some support for this in the literature: mutations in LamC, a LC50 GWAS candidate, trigger cellular redox imbalance, an enrichment of cnc in the cytosol of larval muscle genes, and an increase in the baseline expression of genes shown to be up-regulated by ectopic cnc expression (35, 63).

Conclusions.

Here we set out to explore the genetic architecture underlying a completely novel insecticide chemistry for which there was no expected adaptive precedence because the population sample was collected before chlorantraniliprole was deployed in the field, and is therefore completely naïve to this insecticide and indeed all of the insecticides in the new anthranilic diamide class. This contrasts to several recent studies where insecticide resistance loci have been found to feature genes of major effect that have been built through a series of adaptive substitutions (25, 64, 65). We explored the genetic architecture of chlorantraniliprole survivorship at multiple concentrations and surprisingly found that at higher concentrations there was clearly a gene of major effect. Standing variation in the neuromuscular gene Strn-Mlck increased the LC50 by ∼3 μg/mL. This is not the first time that molecular changes underpinning insecticide tolerance can be considered as exaptations (e.g., ref. 66); however this case unambiguously demonstrates that standing variation in a population is relevant to a future selective agent.

D. melanogaster is not a pest insect but it is common in orchards where these new insecticides are being applied. It may therefore be thought of a sentinel species (19), such that if, in the future, there is evidence for selection at Strn-Mlck or the loci we identified here, then this may be attributed to the use of these new insecticides. Given the precedent for the involvement of the cnc pathway in insecticide resistance in pest species (35), perhaps the factors which cause the constitutive activation of this regulatory hub present the greatest concern.

Materials and Methods

Fly Lines.

All DGRP lines and Transgenic RNAi Project (TRiP) lines (67) were obtained from the Bloomington Drosophila Stock Center. UAS-Cyp12d1-p and 6g1HR-GAL4 were obtained from the Batterham Lab at the University of Melbourne (48). Lines were maintained on cornmeal-yeast-agar media and were kept at 25 °C at constant light for at least one generation before use.

Chlorantraniliprole Phenotyping.

Altacor (350 g/kg chlorantraniliprole/Rynaxypyr) was obtained from DuPont Australia. Wettable granules were dispersed in water to make a stock solution of 100 μg/mL chlorantraniliprole. Chlorantraniliprole was administered through cornmeal-yeast-agar fly media. Insecticide was added to the media once it had cooled below 55 °C. Food dye was added in conjunction with the insecticide to ensure even dispersal through the food. Fifty first-instar larvae were then placed into fly vials, each containing 10 mL of chlorantraniliprole-laced media. This was done in triplicate for each line on each concentration.

A total of six concentrations for each DGRP line were used in the initial screens (0.5, 1, 2, 3, 4, and 5 μg/mL), with a number of lines having to be rescreened on higher concentrations (6, 8, and 12 μg/mL) to accurately determine the LC50. Survivorship was scored strictly 11 d after picking, with survivors being defined as any fly that successfully eclosed.

For calculation of LC50s for both DGRP and transgenic lines, linear models were fitted to concentration-mortality data on a log-probit scale using “glm” in the R statistical package and scripts from Johnson et al. (68). LC50 values and 95% confidence intervals were calculated using Fieller’s method from fitted linear models (69).

GWAS.

GWAS were performed on six single-concentration phenotypes and the LC50 for 152 DGRP lines using the DGRP webtool (70) (http://dgrp2.gnets.ncsu.edu). This involved uploading each of the seven phenotype files and for each running a pipeline that fits a linear model between the phenotype and each site variant among DGRP genomes. The linear model incorporates as covariates Wolbachia pipientis infection status, and five common chromosomal inversion genotypes that vary among the 205 DGRP lines, as described in ref. 70. A total of 1,887,900 variants were tested in each GWAS. To correct for the effect of the most associated variant on the LC50 phenotype, the genotype of this variant was fitted as a fixed effect in a linear model. The residuals were then extracted and submitted to the DGRP pipeline.

Genomic Prediction.

For genomic prediction analysis across the top variants associated with each concentration phenotype and the LC50, a Bayesian linear regression coupled with LASSO was implemented using the Bayesian linear regression package (71) in R, with the phenotype acting as the data vector and the genotypes, the incidence matrix for βL. Up to 500 of the top variants associated with each concentration phenotype (as ranked by P value) were used to explain the phenotype. Each model was implemented with 5,500 iterations, with a burn in of 5,000 iterations and a thinning interval of 50.

Network Based Analysis.

To examine annotated interactions between different GWAS candidates, R spider software (43) (http://www.bioprofiling.de/R_spider.html) was used to map GWAS candidates to KEGG and Reactome pathways.

Transcriptome data for 1- to 3-d-old adult flies from 185 DGRP lines were recovered from the DGRP website (http://dgrp.gnets.ncsu.edu/data.html) (40). Mean transcription level was calculated for each gene in each sex from two biological replicates, to give a mean level for each of the 18,140 transcripts measured by Huang et al. (40) in each DGRP line, for both males and females. The mean of male and female transcript levels was then calculated. A linear model was fit between mean transcription level of each gene measured by Huang et al. (40) and chlorantraniliprole ;LC50.

Strn-Mlck RFLP Mapping Crosses.

To examine the dominance of chlorantraniliprole tolerance and to confirm the association of GWAS-associated Strn-Mlck variants with survivorship, crosses were set up between DGRP lines RAL-399 and RAL-59. These lines have contrasting chlorantraniliprole phenotypes, are free of any of the major cosmopolitan inversions (70), and carry different states at the top GWAS variants in Strn-Mlck. Reciprocal crosses were performed to obtain two classes of F1 progeny to check for maternal and X-linked effects. Offspring from the reciprocal crosses were combined in equal numbers to establish an F2 population.

For the cross-typing, single-fly DNA extractions were performed. Two sets of primers for RFLPs were designed to distinguish “resistant” and “susceptible” Strn-Mlck haplotypes: EcoRV: 5′ TCAGTTCGTTGGTGTTCAGG 3′, 5′ ACTTCACGGTCAACCTGTCC 3′ and HpaII: 5′ GTTCAGATTCATCGCCATCC 3′, 5′ ACCTGCACTACCACGTACCC 3′. A total of 190 flies were scored on both assays. For each sample, 10 μL PCRs were set up using 5 μL of GoTaq Green (9PIM712, Promega), 1 μL each of forward and reverse primers (5 μL), and 3 μL of sterile water.

Knockdown of Strn-Mlck.

RNAi knockdown of different transcripts of Strn-Mlck was carried out using three different UAS lines from the TRiP Library (67), as no single line targets all of the predicted isoforms of Strn-Mlck. The lines, HMS01665, HMS02663, and JF02171 target exons 1, 4, and 30, respectively, and were reciprocally crossed to the ubiquitous actin5c driver. The presence of the CyO balancer chromosome in the Actin5c driver line meant that not all offspring of each cross would inherit the UAS-RNAi construct. For these crosses, the CyO-carrying offspring were used as controls for the UAS siblings.

Overexpression of Cyp12d1-P.

Cyp12d1-p was overexpressed using the GAL4/UAS system (72) and the 6g1HR-GAL4 driver described by Chung et al. (48). The 6g1HR-GAL4 virgin females, in which GAL4 is regulated by Cyp6g1 upstream sequence originating from Hikone-R line flies, were crossed to males carrying an additional copy of Cyp12d1-p under control of a UAS promoter (60). The w line was used as a control for UAS-Cyp12d1-p.

Cyp12d1 qRT-PCR Induction Assays.

For each replicate of each line, 50 third-instar larvae were added to plates containing cornmeal-yeast-agar media and either 2.5 μg/mL chlorantraniliprole, 0.5 μg/mL cyantraniliprole, 10 mM phenobarbital (dissolved in ETOH), or 1.5 mg/mL caffeine (dissolved in 80 °C water). Control plates were untreated. Larvae were allowed to feed for 4 h before being suspended in TRIsure (Bioline) and snap frozen in liquid nitrogen for storage at −70 °C. RNA was later extracted following the standard TRIsure protocol. cDNA was synthesized using M-MuLV reverse transcriptase (NEB) and random nonamer primers following the standard NEB reverse transcriptase protocol. qRT-PCR was performed on the Roche Light Cycler 480 using the following primers: Cyp12d1 (Cyp12d1-p or Cyp12d1-d): 5′ GGGGAAAACTACGATCAGCC 3′, 5′ CGGATTTCCTTAATGCGCTCT 3′ and CG11322 (housekeeper): 5′ TGGGCAGTGCCTTCTACATTT 3′, 5′ CGTACGCACCTCGCTTGTT 3′.