Insights into DDT Resistance from the Drosophila melanogaster Genetic Reference Panel.

Insecticide resistance is considered a classic model of microevolution, where a strong selective agent is applied to a large natural population, resulting in a change in frequency of alleles that confer resistance. While many insecticide resistance variants have been characterized at the gene level, they are typically single genes of large effect identified in highly resistant pest species. In contrast, multiple variants have been implicated in DDT resistance in Drosophila melanogaster; however, only the Cyp6g1 locus has previously been shown to be relevant to field populations. Here we use genome-wide association studies (GWAS) to identify DDT-associated polygenes and use selective sweep analyses to assess their adaptive significance. We identify and verify two candidate DDT resistance loci. A largely uncharacterized gene, CG10737, has a function in muscles that ameliorates the effects of DDT, while a putative detoxifying P450, Cyp6w1, shows compelling evidence of positive selection.

MANY insights into the genetic basis of adaptation have been gained by genetic analyses of resistance to man-made chemicals; whether that be bacteria to antibiotics (Zhang ), weeds to herbicides (Baucom 2016), fungi to fungicides (Schoustra ), or insects to insecticides (Crow 1957; McKenzie 1996). Typically, the genetics of such traits are thought to be strongly monogenic and thus distinct from most other traits whose variation is governed by the combined effect of multiple loci of small effect (Lande 1983; Roush and McKenzie 1987; Allen ). This contention is central to the debate over the genetics of DDT resistance in Drosophila melanogaster. As it is not considered a pest, this model insect may have never been exposed to the high selection intensities thought necessary for the evolution of major-gene-based resistance (Macnair 1991; McKenzie and Batterham 1994). Despite most genetic investigations supporting a polygenic architecture (Crow 1954; King and Sømme 1958; Dapkus and Merrell 1977; Shepanski ) or more precisely weak polygenicity (Dapkus 1992; i.e., at most one or two loci per chromosome arm), a gene of major effect, , was identified and described as “necessary and sufficient” for DDT resistance (Daborn ). Proposed solutions to this apparent contradiction between “polygenic” and “monogenic” schools include delineating a difference between field and laboratory evolved resistance (Ffrench-Constant 2013) or suggesting that is only inconsistently associated with resistance (Kuruganti ).

In the specific case of the long-term DDT-selected strain 91-R, there is a strong argument for the polygenicity of DDT resistance. 91-R is over 1500 times more resistant to DDT than another laboratory strain, Canton-S (Strycharz ); levels far beyond those attributable to alone. Much of the recent work aimed at identifying polygenes associated with DDT resistance has focused on contrasting these resistant and susceptible laboratory strains and has allowed the generation of transcriptomic (Pedra ), proteomic (Festucci-Buselli ), physiological (Strycharz ), and genomic (Steele , 2015) data sets of differentiated factors. These candidates, however, remain to be validated, and their contribution to either the selection response in the laboratory, or to field resistance, have yet to be quantified. Until now has remained the only DDT resistance locus with molecularly defined alleles that contribute to DDT-related genetic variation observed in field populations.

Previously Schmidt demonstrated that an allelic series at the locus is sweeping through D. melanogaster populations. The resistant variant described by Daborn , which harbored a partial Accord transposable element insertion upstream of the locus, was subsequently found to be invariably associated with a duplication of the locus (Schmidt ). It is still unclear to what extent the upregulation of observed at this allele is due to the Accord element (Chung ) or to the presence of multiple copies, or both. However, it is clear that the ancestral allele (Cyp6g1-M) is at very low frequencies throughout the world and has largely been replaced by haplotypes carrying the duplication and various transposable element insertions. The Cyp6g1-AA allele (which has two copies of the gene both bearing the Accord element) and the Cyp6g1-BA allele (where one of the duplicated copies contains an additional HMS-Beagle transposable element insertion) both confer resistance to DDT relative to the Cyp6g1-M allele, although they show only subtle phenotypic differences between one another. A fourth allele, Cyp6g1-BP, is at high frequencies in populations in northeastern Australia, and provides a further increase in resistance, in both a set of isochromosomal lines and in a field population (Schmidt ). In the latter, the Cyp6g1-BP allele accounts for ∼16% of the total variation in the DDT resistance phenotype.

The analysis of field variation at the locus provides two arguments suggesting that other loci contribute to DDT resistance in the field. First, as allelic variation does not account for all the 20–50% heritability expected for insecticide resistance (McKenzie 2000), there is much genetic variation that is not accounted for (20–68%). Second, as a succession of alleles provides a series of adaptive events at one locus, there has been time and selective pressure for other adaptive events arising elsewhere in the genome.

DDT affects the nervous system by targeting the voltage-gated sodium channel resulting in uncontrolled nerve firing. In D. melanogaster, the α-subunit of the voltage-gated sodium channel is encoded by , and target-site resistance to DDT has been described in mutants by Pittendrigh and Lindsay . However, variation affecting DDT resistance has not yet been identified in outbred populations of D. melanogaster, nor has it been implicated in the phenotypes of well-studied DDT-resistant laboratory strains (Dapkus and Merrell 1977; Steele ).

DDT-induced mortality is typically preceded by the loss of peripheral muscle control and manifests as a “knockdown” phenotype in which the fly lies paralyzed with appendages twitching, or exhibits uncontrolled flight. Here we use the Drosophila Genetic Reference Panel (DGRP; Mackay ; Mackay and Huang 2017), a resource for the dissection of quantitative traits, to characterize the genetic architecture of both knockdown and mortality caused by DDT exposure using predominantly field-derived genetic variation from a single population. This population has shown great phenotypic diversity for, and furnished genetic associations with, a diverse range of traits including insecticide resistance (Mackay ; Huang ; Battlay ). We performed genome-wide association studies (GWAS) on two DDT-related traits to gain insights into insecticide biology and insecticide-driven selection using a classic insecticide model.

Materials and Methods

Fly stocks

DGRP lines were obtained from the Bloomington Stock Center, and were maintained on corn meal media within a 22–25° temperature range. For the DDT assays, larval density was controlled by setting up vials with 10 mated females that were allowed to lay for 3 days, then transferred to fresh media for a further 3 days to establish an additional brood to provide enough progeny for replication. At this stage, two to three replicates per line were established. Adults were allowed to eclose for 2–3 days, and females were collected to food-containing holding vials after brief CO2 anesthetization. After 2–3 days’ recovery, adult females were assayed, thus the age range was 4–6 days for flies assayed.

DDT assays

Glass scintillation vials were inoculated with 200 μl of Acetone/DDT solution at a concentration of 0.5 μg/ml, giving a final contact assay amount of 100 μg of DDT. Cotton wool moistened with 10% sucrose solution was used to stopper the scintillation vials. Assays were set up between 7 and 8 am, roughly corresponding with the onset of the light phase of the diurnal cycle. Time of setup was recorded, and each vial was scored for the presence of flies exhibiting either knockdown or mortality at ∼1-hr intervals. Knockdown here is defined as either flies permanently seen to be in a prone position with jolting movement of legs or wings, or a prolonged inability to right themselves from a prone position after tapping of the scintillation vial on the laboratory bench. Mortality was determined as the point at which all external movement ceased. Data were recorded in data sheets with time points as columns and DGRP line number as rows. Assay vials contained 10 females, and there were at least three replicates per DGRP line. A total of 184 DGRP lines were assayed. In total, the assay time was 30 hr. Phenotypes were scored through 1–15, 24, and 30 hr.

Estimation of heritability

At least three replicates of ten flies per vial were analyzed per DGRP line, which allowed the calculation of broad-sense heritability (H2). Following the method of Clowers , heritability was estimated from the variance components of a linear model of the form: Phenotype = Population mean + Line effect + error. As these lines are inbred, the components of the Genetic Variance (e.g., additive or dominant), cannot be separated, thus these estimates are of broad-sense (H2) heritability. An ANOVA was performed in Statistics Package in Social Sciences (SPSS), and the variance components were estimated as the Mean Sum of Squares. Total phenotypic variance was estimated as Genetic Variance + Environmental Variance. H2 was thus estimated as Gv/Gv+Ev. This was calculated for each discrete time point.

GWAS

Phenotypes were submitted to the Mackay Laboratory DGRP2 website for GWAS (http://dgrp.gnets.ncsu.edu/; Huang ). To perform our own GWAS for DDT resistance traits, we used the PLINK GWAS program (PLINK v1.9, https://www.cog-genomics.org/plink2; Chang ). For genotypes, we initially downloaded the DGRP freeze 1 SNP tables (http://www.hgsc.bcm.tmc.edu/projects/dgrp/freeze1_July_2010/snp_calls/). These tables encoded “haploid” variant calls, with heterozygous positions encoded with the relevant ambiguity codes. To make PLINK-compatible files, we used perl scripts to modify these SNP tables into diploid genotypes. We were surprised by a large number of sites that contained more than two alleles. The initial DGRP freeze 1 GWAS tool removed all such sites, but assuming that many of these represented low-frequency sites or errors, we ranked each variant by frequency and kept the two most frequent ones. The one major difference between our initial results and that of the DGRP freeze 1 GWAS tool was the triallelic variant at .

However, our analyses here are based on the newer DGRP freeze 2 release. For this we downloaded the dgrp2.bed file (PLINK format) from the DGRP2 website (http://dgrp2.gnets.ncsu.edu/data/website/dgrp2.bed). Thus, these genotypes are identical to those used in the DGRP webtool, except for choices in variant filtering as indicated in the PLINK command line below. It should also be noted that now no heterozygous genotypes are present, and that the two most common variants at triallelic sites are also present in this bed file.

We used the following base PLINK command to perform GWAS (linear regression of genotype phenotype association assuming an additive model of allelic effects) after initially filtering for SNPs with MAF ≥ 0.05 and at least 70% genotyping rate (given that a triallelic site with equal proportions of all three genotypes would have a rate of ∼66% in the DGRP freeze 2 data):

plink–allow-extra-chr–allow-no-sex–bfile dgrp2–covar dgrp2_ESTRAT_PCA20.txt–linear–map3–no-fid–no-parents–no-pheno–no-sex–aperm 5 1000000

The “aperm” option was used to generate empirical estimates of the P-value for each SNP. To calculate the family-wise error rate (FWER) for a given P-value or level of significance, we generated 5000 random permutations of the 4-hr knockdown (4 hrkd) and 24-hr mortality (24 hrm) phenotype data and performed the same linear association model as for the observed phenotypes. For each permutation we record the lowest observed P-value. The FWER is then the number of random permutations that generated a P-value lower than or equal to that observed for each SNP, i.e., that chance of at least one false positive at this level of significance. After filtering, 1,776,058 SNPs were tested for each PLINK GWAS.

We performed principal components analysis on the DGRP data using PLINK after prefiltering genotypes for minor allele frequency (>0.05), missing rate (<0.70), and linkage disequilibrium (LD) pruning (r2 < 0.2) using the indep pairwise command, with both window and step size of 500 variants. To control for confounding cryptic relatedness in the DGRP, we used the first 20 principal components as covariables in all GWAS, as indicated above with the “–covar dgrp2_ESTRAT_PCA20.txt” option in the PLINK command. As inversions and Wolbachia infection status can also influence the phenotypes of the DGRP lines, we used the phenotypes adjusted for these factors outputted from the DGRP2 website. Thus, our results for the DGRP freeze 1 and DGRP freeze 2 differ due to now correcting for population structure, adjusting phenotypes for inversion and Wolbachia infection status, and the removal of heterozygous genotypes in the DGRP freeze 2. Finally, for annotating variants we used the annotations made available on the DGRP2 website (http://dgrp2.gnets.ncsu.edu/data/website/dgrp.fb557.annot.txt).

To test specific variants at , we used PCR primers described in Schmidt and/or local alignments of Illumina and 454 reads to determine the genotypes. For , our DGRP freeze 1 GWAS indicated the presence of a triallelic site at 2R:6,174,944[V6]. We used the reprocessed DGRP genotype data from the Drosophila Genome Nexus (Lack ) files to extract genotype data for this position (described further under the iHS and nSL tests methods), as the genotypes for the third variant state are removed from the DGRP freeze 2 data. For both these genes, the sites of interest are triallelic in the DGRP. We generated files with the ancestral allele and derived alleles, both singularly and combined after collapsing to the same allelic code, included to test for associations with DDT phenotypes. In the particular case of , we introduced a variable site to code Cyp6g1-M, Cyp6g1-AA, or Cyp6g1-BA alleles.

iHS and nSL tests of recent positive selection

We downloaded the haploid FASTA alignments for 205 DGRP2 genomes, described by Lack , that are contained in the Drosophila Genome Nexus (DGN; http://johnpool.net/genomes.html). We applied the supplied masking filter for identity by decent (IBD) regions, but not the masks for admixture. For ease of downstream processing, we converted these FASTA alignments into Variant Call Format (VCF) files using custom bash and perl scripts. Note, these scripts kept the genotype information at triallelic sites.

After IBD masking and the prefiltering of excessively heterozygous regions in the DGN data (using the scripts available from http://johnpool.net/genomes.html), the genomes varied drastically in the proportion of sites with missing data. We excluded chromosomes with >20% missing data, leaving us with 134 lines for chromosome 3L. For other chromosomes, we randomly sampled lines without replacement to n = 134. In the compilation of VCF files, we excluded sites with >15% missing genotypes in the DGRP2. There were of course missing data remaining after these two steps, still rendering them unsuitable for haplotype tests. We filled in missing genotypes by randomly sampling from called genotypes, in proportion to their population frequency. In the main, this procedure will reduce local LD, and thus is conservative vis-a-vis LD-based tests of selection, and is similar to the method employed by Garud .

We used two tests of selection that both utilize signatures of LD around a beneficial allele. This LD can be measured by the extended haplotype homozygosity (EHH) statistic (Sabeti ). EHH summarizes the probability that two chromosomes bearing the beneficial allele are also identical by state at a nearby neutral variant. Increased LD will thus also increase the EHH statistic. The integrated haplotype score (iHS; Voight ) extends EHH by calculating the area under the EHH curve (iHH, and thus integrated) for both the derived and ancestral variant, and then finding the ratio of iHHDerived/iHHAncestral. If LD for both the derived and ancestral allele is approximately the same, this ratio should be close to 1, whereas large departures could indicate the action of recent positive selection. As LD is also influenced by the age of mutations, iHS scores are standardized in bins of derived allele frequency, as under a neutral model frequency is indicative of allele age (on average). nSL (number of segregating sites by length; Ferrer-Admetlla ) is similar in spirit to iHS, but instead of EHH, the base statistic is the average count of SNPs for which two haplotypes are identical. We use it here because this statistic may be less biased toward regions of low recombination and has an increased ability to detect soft sweeps (i.e., the beneficial allele is segregating on more than one distinct haplotype at the time of selection).

To calculate both iHS and nSL, we used the versions implemented in selscan with default parameters for each (Szpiech and Hernandez 2014). We kept biallelic SNPs with 5% < DAF < 95%. The ancestral state was estimated as the homologous D. simulans reference genome allele (aligned to D. melanogaster, also downloaded from http://johnpool.net/genomes.html). Both statistics were normalized in 1% frequency bins per chromosome, and a P-value per site was calculated from the empirical distribution of normalized scores (|iHS| |nSL|) per chromosome. Genetic map positions were estimated using the recombination maps generated by Comeron .

For the triallelic sites at and , we manually constructed haplotype files containing the set of one of the derived haplotypes and the ancestral haplotypes. These sets are a subset of those for biallelic sites, and P-values for each of these were calculated as above.

For the extended haplotype structure for the _GLY allele, the EHH curve does not decay below 0.05 within the default window size, a reflection of both increased haplotype structure but also a region of increased heterozygosity downstream of , which is filtered out in most DGRP2 individuals and introduces a gap >200 kb. To compute an iHS value for this variant we relaxed the default parameters, so that the max gap is 250 kb. Manual inspection of the EHH curves (Figure 4) suggests that the high iHS statistic for this variant is due to haplotype structure and not this gap in the estimated haplotypes.

Figure 4

Extended Haplotype Homozygosity (EHH) plots: (A) around the Cyp6w1 locus using the triallelic second site of codon 370 as the focal variant and (B) around the Cyp6g1 locus using transposable element insertion sites as the focal variant.

CG10737 RNAi

Line CG10737-KK [Vienna Drosophila Resource Centre (VDRC) ID 106383] that has an upstream activator sequence (UAS) followed by a hairpin sequence matching CG10737 with no recorded off-targets was obtained from the VRDC. It was crossed to Mef2-GAL4 (Bloomington ID 27390). Line y, w[1118]; P(attP,y+,w[3′]) (VDRC ID 60100) was used as control and also crossed to the Mef2-GAL4 driver. The offspring of these crosses were evaluated for the DDT resistance phenotypes, 3hrkd, and 24hrm. The screens for each cross were performed over a range of DDT concentrations using replicates of 10 females aged 4–6 days old. For every DDT concentration, a minimum of 80 females were screened for 3hrkd and a minimum of 50 females were screened for 24hrm. The relationship between the 3hrkd resistance phenotypes and DDT concentration was evaluated using XLSTAT to perform a Probit analysis, a log normal regression of the knockdown data (Sakuma 1998). Dosage mortality curves were constructed for the relevant crosses and their controls and the ED50, the dose that will theoretically knockdown 50% of the population, was estimated. To ascertain that the gene had indeed been knocked down, total RNA was extracted from 20 to 30 adult females using the TRIsure-reagent. The extracted RNA was treated with RNA-free DNase (New England Biolabs, Ipswich, Massachusetts) and residual genomic DNA contamination was assessed for each sample by attempting to amplify a PCR product from the RNA. The samples were converted to complementary DNA (cDNA) using 1 µg of RNA in 20 µl reaction with Anchored Oligo (dT)20 and MuLV Reverse Transcriptase (NEB) according to the manufacturer’s instructions.

Cyp6w1 transgenic overexpression

UAS-Cyp6w1 constructs were designed using the y; cn bw sp; reference strain coding sequence as a template. Constructs containing each state of triallele identified in the 24 hrm GWAS were created by mutating the VAL370 codon in the reference sequence to ALA370 and GLY370. All three constructs were manufactured by Biomatik USA (Wilmington, DE) and were supplied cloned into plasmid pUASTattB.

The UAS-Cyp6w1 constructs were transformed into the y{vas-int.Dm}ZH- 2A w*; (EPS) M{36P3-RFP.attP}ZH-86Fb recipient strain, which has a defined integration site on Chromosome 3, using the attP–attB system and QC31 integrase (Bischof ). To generate heterozygous Cyp6w1/tubulin-GAL4 or serrate flies, virgin female homozygous UAS-Cyp6w1 flies were collected on light CO2, stored at 22–25° for 4 days to ensure virginity, and then 10 females per vial were crossed to tubulin-GAL4/serrate males. At the same time, vials with either 10 UAS-Cyp6w1 or tubulin- GAL4/serrate females were established. Adults for all these fly lines were allowed to eclose for 2–3 days, and females were collected to food-containing holding vials after brief CO2 anesthesia.

Each line was tested on a minimum of five doses of DDT. A minimum of five replicates of 20 individuals per dose per line were treated in the bio-assay. Glass scintillation vials were inoculated with 200 μl of Acetone/DDT solution with concentrations ranging from 5 × 10−5 to 1.0 μg/μl, giving a final contact assay range of 0.01–200 μg of DDT. Cotton wool moistened with 10% sucrose solution was used to stopper the scintillation vials. Flies were assayed at 25° for 24 hr. Data for each tested line were individually analyzed using PriProbit (Sakuma 1998) to estimate the LC50 and generate data for plotting dosage mortality curves and 95% confidence intervals.

91-R and 91-C genome examination

bam files containing alignments of 91-R and 91-C sequencing reads to the y; cn bw sp; reference genome were obtained from the sequence read archive (Leinonen ; SRR1237973; SRR1237974). Alignments at relevant loci were visualized using IGV software (Thorvaldsdóttir ) to determine presence of SNPs and discordant reads consistent with structural variation.

Data availability

Raw phenotypic data are available in the Supplemental Material, File S1.

Results

One hundred eighty-four lines of the DGRP were exposed to a dose of 100 μg of DDT and monitored over 30 hr. Two phenotypes were recorded; percentage knocked down and percentage dead. The median of the 50% knockdown time for the DGRP lines was 15 hr, while the median of 50% mortality time was 18 hr.

We chose to focus on 24hr mortality (24hrm), as it is a standard assay in D. melanogaster insecticide resistance literature (Daborn ; Schmidt ). We also observed minimal zero-dose control mortality and maximal broad-sense heritability (H2 = 0.8) at this time point. Of our knockdown timepoints, we chose 4 hrkd for further analysis, as the H2 of our assay reached an early peak at this time point and it thus provides a robust yet contrasting data set to the 24hrm (Figure 1A and File S1). Based on the number of lines phenotyped and genotyped (n = 179), the number of replicates (r = 3), and the broad-sense heritability of the trait (H2 = 0.8), the method of Mackay and Huang (2017) was used to estimate the power to detect variants for a range of effect sizes. Variants with an effect of 12% of H2 (∼0.1 units of the trait’s SD) would be detected with 50% chance and variants with an effect of 22% of H2 (∼0.18 units of the trait’s SD) with 95% chance, respectively (File S4).

Figure 1

DDT-Knockdown and mortality. (A) The mean percentage 4-hr knockdown (4hrkd) and 24-hr mortality (24hrm) for each DGRP line (based on at least 3 replicates of 10 individuals each) with lines ordered from lowest to highest for each trait. (B) The scatterplot shows the correlation between 24-hr mortality and 4-hr knockdown. Flies that are knocked down early can recover and exhibit late mortality.

The 4 hour knockdown (4hrkd) phenotype may reveal the intrinsic variation in insecticide uptake routes and nervous system sensitivities, while the later 24hrm phenotype may encompass variation in the physiological response to the insecticide such as inducible detoxification mechanisms. Four-hour knockdown and 24hrm show distinctly different distributions within the DGRP. The population mean for 4 hrkd is 21% with SD of 19%, while 24 hrm has a mean mortality of 58%, with SD 29%. The genetic correlation between the two traits was estimated to be 0.42 (Figure 1B).

We performed GWAS on our 4hrkd and 24hrm phenotypes using both the DGRP2 webtool (Huang ) and a custom implementation using PLINK (version 1.9). We did not observe a significant (ANOVA P > 0.05) effect on DDT-induced knockdown or mortality of Wolbacchia infection (knockdown: r2 = 0.013 and mortality: r2 = 0.010), inversions with a frequency >5% (knockdown: r2 range = −0.145–0.176 and mortality: r2 range = −0.264–0.136), or genome size (knockdown: r2 = 0.005 and mortality: r2 = 0.004). All four GWAS detected variants with P-values less than the arbitrary genome-wide significance threshold (1 × 10−5; File S2; Mackay ). None of the variants associated with 24 hrm were also associated with the 4 hrkd phenotype. Only 24 hrm with the PLINK pipeline showed a variant with FWER <0.05, in a cytochrome P450 gene . Using a less conservative cutoff of P < 1 × 10−6, 4hrkd showed eight variants. Of note, one of these is in CG10737, a gene that has previously been implicated in DDT resistance.

CG10737

Two SNPs separated by 30 nucleotides occur in the intron of and are highly associated with 4hrkd in both PLINK and DGRP2 GWAS; they are in complete LD with each other and have slightly different P-values of association because of missing data (2R:19,169,399[V6], 2R:19,169,429[V6]; File S2). is predicted to encode a protein with C1 and C2 domains that are capable of diacylglycerol binding and calcium-dependent targeting, respectively, and is therefore postulated to be involved in intracellular signal transduction (Mitchell ; Attrill ). A previous study highlighted that CG10737 is one of five “lipid metabolism” genes shown to be differentially expressed among two DDT-resistant lines and a susceptible laboratory line. Specifically, was expressed at significantly lower levels in the DDT-resistant Wisc-1 line than the susceptible Canton-S line (Pedra ). This led us to hypothesize that a reduction of transcript would lead to increased resistance to DDT-knockdown. To test this hypothesis, we crossed a VDRC UAS-CG10737-KK line to knockdown this gene using a Mef2-GAL4 driver, which targets the muscles. The justification for this driver is fourfold: (i) according to the modENCODE data sets is highly expressed in the muscle-associated tissues (crop, hindgut, heart, carcass, and male accessory gland) and is considered muscle-enriched in two microarray analyses (Zhan ; Weake ); (ii) genes that are coexpressed with are known muscle genes (Myosin heavy chain, myosin light chain 1 and 2, troponin, and myosuppressins); (iii) chromatin co-immunoprecipitation experiments show that MEF2 binds to the CG10737 cis regulatory region (Sandmann ; Cunha ); and (iv) Bonn (2010) demonstrated that antibodies to CG10737 stained less intensely in the somatic and visceral mesoderm when in a mef2 null background.

The progeny of the Mef2-GAL4 × UAS-CG10737-KK cross were assayed for DDT-knockdown resistance over a range of doses. The dose of DDT required to knockdown 50% of flies (ED50) was more than twice as high for UAS-CG10737-KK lines as it was for the control lines (109 μg/vial vs. <47 μg/vial; Figure 2). These findings support the hypothesis that transcript reduction increases resistance to DDT-knockdown.

Figure 2

RNAi knockdown confirms that lowering CG10737 transcript abundance increases DDT resistance and that this effect can be mediated by manipulating muscle expression. Flies in which CG10737 is knocked down in muscles using the mef2-GAL4 driver crossed to the VDRC CG10737-KK line are significantly more resistant to DDT than control flies, which are the result of substituting the relevant genetic background line for either mef2-GAL4 (w) or the CG10737-KK line (60100) in the cross.

Cyp6w1

This gene is the standout candidate as judged by significance of association with variation in 24hrm. The associated SNP (2R:6,174,944[V6]) results in a coding sequence change in an enzyme belonging to a family capable of detoxification and which is most abundantly expressed in the fly fat body, a known detoxification tissue. It is also expressed in the spermatheca, heart, head, and carcass (Chintapalli ), and Pedra found it to be overexpressed in the DDT-resistant 91-R line relative to the Canton-S susceptible line. The variant site associated with DDT 24hrm in is particularly interesting as three different allelic states at this site are observed within the DGRP, and each state encodes a different amino acid. Such triallelic sites, which represent 2.5% of variable sites in the DGRP, get excluded by many analyses (either explicitly, by excluding sites with more than two states, or implicitly, because the most abundant two states combined are not enough to pass abundance thresholds).We surmise that this latter factor is the reason we do not recover this variant using the Mackay DGRP2 webtool (Huang ; http://dgrp.gnets.ncsu.edu/), as <75% of lines have a called genotype for this position. In contrast, we used a lower missing genotype threshold of 60% in our PLINK GWAS pipeline. We initially identified 2R:6,174,944[V6] as a biallelic TC transition resulting in a Val370Ala substitution. By investigating this site in the Drosophila genome DGRP2 data, we found a third state, G, that is at a frequency of 14% in the DGRP, and results in a Val370Gly substitution. Combining the effect of the C and G alleles, the P-value for association with DDT mortality is 7.88 × 10−7. Individually the G allele is insignificant genome-wide (P = 0.015), but taken together these results suggest that it too could influence DDT resistance.

To test what effect the three amino acid states at site 370 of Cyp6w1 have on DDT resistance, three transgenic lines were generated such they that only differed at that single site. The three lines were crossed to a GAL4 driver line under the control of tubulin promoter sequence, so that the UAS- isoforms were expressed in a ubiquitous manner. All three of the transgenes were at least 12 times more resistant to DDT than their intracross sibs when under the control of the tubulin-GAL4 driver [Figure 3; LD50 (the dose predicted to result in 50% mortality) for overexpression lines: Cyp6w1_VAL = 14 μg/vial, Cyp6w1_GLY = 12 μg/vial and Cyp6w1_ALA = 28 μg/vial]. This indicates that overexpression of in itself can result in a DDT resistance. However, an additional effect is seen for the Ala370 amino acid substitution, as flies overexpressing this isoform are at least two times more resistant than either Val370 or Gly370 counterparts (Figure 3). Thus, the Ala370 substitution appears to increase DDT resistance, supporting the results of our 24hrm GWAS.

Figure 3

Toxicology of transgenic lines overexpressing the three allelic forms of Cyp6w1 (far right) and various controls. All three allelic forms of Cyp6w1 confer DDT resistance relative to controls, and Cyp6w1_ALA confers a further significant increase in DDT resistance relative to the other two alleles.

As our transgenic analyses showed that overexpression can yield resistance and transcripts are more abundant in Pedra , we examined transcription in females from a DGRP transcriptome data set (Huang ). However, there was no strong correlation between transcription level and DDT mortality or knockdown data, nor between transcription level and the triallelic states of . Similarly, we did not find any of the other eight detoxification or four lipid metabolism genes that were focused on by Pedra to have strong correlations between the DGRP transcription data sets and either DDT traits we examined (File S3).

Selective sweep tests at DDT-associated loci

To assess whether 4 hrkd- and 24 hrm-associated (P < 1 × 10−5) variants from the PLINK GWAS exhibited evidence of being the targets of selection, we used the integrated Haplotype Score (iHS) test (Sabeti ; Voight ). We also calculated the related statistic nSL (number of segregating sites by length; Ferrer-Admetlla ), which defines haplotype lengths via pairwise differences and may have higher power to detect soft selective sweeps. The only derived variant yielding a significant iHS or nSL score was the Glycine allele of (|iHS| = 3.82, P-value = 0.00085; Figure 4; |nSL| = 3.04, P-value = 0.0028). The other derived variants did not show significant departures from expectation, including the second derived allele at – the alanine allele that confers increased DDT resistance in our transgenic experiments.

To ascertain the power of these iHS and nSL approaches with this data set we also examined , the genomic region around which Catania , Schlenke and Begun (2004), and Garud noted a reduction in genetic diversity indicative of a selective sweep. While there is an allelic series at this locus (Schmidt ), the two predominate alleles in the DGRP are Cyp6g1-AA and Cyp6g1-BA, and the EHH for both of these alleles and for the ancestral Cyp6g1-M allele was determined. Clearly both derived alleles exhibit large extended haplotypes. Both of these are significant compared to ancestral Cyp6g1-M allele (Figure 4; Cyp6g1-AA: |iHS| = 2.45, P-value = 0.018; |nSL| = 2.61, P-value = 0.011; Cyp6g1-BA: |iHS| = 1.96, P-value = 0.027; |nSL| = 2.47, P-value = 0.012).

Further evidence that insecticide resistance alleles may be adaptive comes from elevated levels of population differentiation (Taylor ). We examined allele frequencies in African, European, Australian, and North American populations (Pool ; Martins ; Reinhardt ; Bergman and Haddrill 2015; Table 1); most populations indicate segregation of each of the triallelic states, except for Portugal, which does not exhibit Ala370. The frequency of the variants, however, shows marked population differentiation, perhaps indicating geographical variation in selection intensity by DDT or other insecticides.

Table 1

Population frequencies of Cyp6w1 codon 370 trialleles

Continent	Population	Data Source	Sample Size	Allele Frequency
				VAL	ALA	GLY
Africa	Combined	Pool et al. (2012)	139	0.9	0.07	0.3
Australia	Queensland	Reinhardt et al. (2014)	17	0.53	0.03	0.46
	Tasmania	Reinhardt et al. (2014)	15	1	0	0
Europe	France – A	Bergman and Haddrill 2015	50	0.78	0.02	0.2
	France – B	Pool et al. (2012)	8	0.5	0	0.5
	Portugal	Martins et al. (2014)	12	0.95	0	0.05
North America	Maine	Reinhardt et al. (2014)	16	0.69	0.25	0.06
	North Carolina	Mackay et al. (2012)	162	0.54	0.32	0.14
	Florida	Reinhardt et al. (2014)	16	0.63	0.31	0.06

Discussion

The DGRP is a powerful resource for identifying natural variation contributing to a range of phenotypes (Mackay and Huang 2017). However, it is of particular interest to the study of insecticide resistance due to the important role insecticides have apparently played in the evolutionary history of this population; two of the strongest peaks of selection in the DGRP, genome wide, are and (Garud ).

is the molecular target of organophosphate and carbamate insecticides, and within the DGRP segregate three of the four substitutions in known to confer resistance to these insecticide classes in D. melanogaster (Menozzi ; Battlay ). has likewise been demonstrated to provide resistance to insecticides, including the neonicotinoids imidacloprid and nitenpyram (Daborn , 2007; Joußen ) and the organophosphate azinphos-methyl (Battlay ). There is also evidence that may contribute to resistance in other organophosphates including malathion, parathion, and diazinon (Kikkawa 1961; Pyke ). But the best studied of ’s insecticide substrates remains DDT (Daborn , 2007; Joußen ; Schmidt ); derived alleles of have previously been shown to contribute significantly to DDT resistance, with the Cyp6g1-BP allele explaining ∼16% of the total variance in DDT resistance in a field population from Queensland, Australia (Schmidt ).

While DDT is not used in much of the world today, DDT resistance is a trait of interest for two important reasons: first, it may help to explain the strong and recent sweeps at the locus (Schmidt ; Kolaczkowski ; Garud ), and second it was the subject of the classic work by Crow (1956), Sokal and Hunter (1954), and others on microevolution. We performed GWAS of two DDT resistance phenotypes in the hope of dissecting the quantitative nature of a selective response to DDT exposure.

Cyp6g1 was not identified in the GWAS

The association of allelic variation with either the 4 hrkd or 24 hrm DDT phenotypes is not significant, even at the genome-wide significance threshold. This is partly attributable to the sweep itself; resistant Cyp6g1-AA and Cyp6g1-BA alleles are present in all but nine DGRP lines (all nine were phenotyped in this study), while the Cyp6g1-BP allele, which confers the most drastic increase in DDT resistance (Schmidt ), is absent from the population. Although much of the power to detect the effect of resistance is gone from the DGRP, Battlay have previously shown that a haplotype in LD with the ancestral Cyp6g1-M allele was strongly associated with susceptibility to a low dose of the organophosphate insecticide azinphos-methyl in DGRP larvae. One reason was associated with the azinphos-methyl phenotype and not our DDT phenotypes could be the relative contributions of overexpression to these traits: Daborn demonstrated a threefold increase in adult DDT LD50 when was overexpressed using the GAL4-UAS system, while Battlay observed a 6.5-fold increase in larval azinphos-methyl LD50 in the same crosses. It is also possible that the dose of the DDT used in this study did not maximize the chances of recovering the variation present in the DGRP. Steele also reported an absence of signal in their comparison of 91-R and 91-C genome sequences. This is due, however, to Steele limiting their study to analysis of open reading frames. We revisited genome sequence alignments generated by Steele and found that evidence of the Accord LTR insertion into the 5′ UTR of (Daborn ) as well as copy number variation in both and (Schmidt ), both identifiers of resistance alleles, are present in 91-R but not 91-C.

The two detoxification enzymes implicated in DDT mortality, and , could interact in parallel or in series. Epistatic interactions between alleles at the two loci were examined but, unfortunately, we did not have the power to test for such interaction in the DGRP because there is only one DGRP line (Ral-486) that is homozygous for the susceptible Cyp6g1-M allele and homozygous for the resistant Cyp6w1-Ala allele. We have, however, assessed pairwise epistatic interactions of the GWAS nominally associated variants listed in the combined mortality and knockdown lists, and none were significantly extreme to pass the Bonferroni correction threshold (∼0.000015) for either the knockdown (lowest value: 0.0002) or the mortality (lowest value: 0.0015) traits.

Variation at CG10737 contributes to 4-hrkd but not 24-hrm

The GWAS presented here demonstrate that the two DDT-associated traits under examination have distinct genetic architecture, and for each trait we have generated transgenic evidence to support the involvement of a top candidate in DDT phenotype variation. has previously been implicated as a DDT resistance candidate, being one of 158 genes shown to be significantly differentially regulated between DDT-resistant 91-R and a standard susceptible laboratory stock (Canton-S; Pedra ). Moreover, inspection of the 91-R and 91-C genome sequences generated by Steele shows that both GWAS SNPs, associated with increased knockdown susceptibility, are present in 91-C but absent in 91-R. Using a completely independent approach we have found that these noncoding variants in , present in ∼10% of the DGRP, are strongly associated with 4 hrkd. Our RNAi knockdown of increased DDT resistance, consistent with Pedra observation that the transcript was less abundant in the resistant 91-R line. We note that the significantly associated GWAS variants also lie in a region that modENCODE found to bind the Trithorax-like () GAGA factor, suggesting an untested mechanism linking the naturally occurring polymorphisms to transcriptional abundance.

An important insight gained through these experiments is that they demonstrate a role of in the muscles of flies, and a relationship between DDT and muscle biology. knockdown in the muscles resulted in a two- to fourfold increase in resistance, demonstrating that DDT perturbation can be ameliorated by genetic manipulation limited to the muscles.

The molecular function of has not been well characterized. It is a gene with 14 coding exons and produces at least 12 alternatively spliced mRNA isoforms. It putatively encodes a transmembrane protein with C1 and C2 calcium/lipid-binding domains and it has been annotated as “protein kinase c-like” although it lacks a catalytic domain. C1 domains have cysteine motifs that form zinc fingers and typically bind diacylglycerol (Colón-González and Kazanietz 2006), while C2 domains typically bind calcium ions and are associated with protein–protein and protein–membrane interactions (Ochoa ). Given the important role that calcium plays in synapse and muscle signaling, perhaps a decrease in the expression level of makes flies less sensitive to Ca2+ ion signaling. It thereby may reduce the expression of the DDT-knockdown phenotype by mediating the effect of overstimulation of muscle cells caused by the uncontrolled firing of nerves due to DDT’s effect on the voltage-gated sodium channel, Para, rather than influencing nerve functionality directly.

While is associated with DDT-knockdown in the GWAS, it was not significantly associated with 24hrm GWAS. Furthermore, a reexamination of the DDT-knockdown phenotype for DGRP lines over the initial time course of the experiment found that variants in this gene were only strongly associated with knockdown over a narrow time window (3–5 hr), suggesting that the gene’s role in knockdown resistance is transient. Additionally, in the RNAi experiments, the dose that effectively interrogated the temporally defined knockdown was so high that the flies did not live long enough to accurately score 24-hrm. This is attributable to genetic background differences between the lines going into the RNAi crosses and the DGRP lines. These issues highlight the differences in genetic architecture of subtly different DDT-related traits. While the knockdown trait may give us insight into the effect of perturbations of DDT on insect physiology, and the function of previously uncharacterized genes, one has to ask how relevant this variation is to the field exposure to DDT. We were unable to detect any signature of selection around variants in this gene, or indeed any high-frequency derived knockdown-associated variants, that would indicate microevolutionary change in the recent history of the Raleigh population of D. melanogaster – the source of the DGRP. However, this may reflect an issue of power to detect subtle changes in frequency of preexisting variants rather than an irrelevance of the phenotype to survivorship in the field. We also note that in early DDT selection experiments in D. melanogaster by Sokal and Hunter (1954), pupation height was correlated with DDT resistance and that Riedl mapped a pupal height QTL to a narrow region encompassing .

Triallelic variation at Cyp6w1 separates DDT-resistant allele from a selective sweep signal

Like , was previously identified as a DDT resistance candidate by Pedra , who found it had high transcriptional abundance in 91-R relative to the susceptible Canton-S. Again, inspection of the Steele genomes revealed that 91-R also carries the DDT resistance-associated Cyp6w1_ALA allele, whereas 91-C carries the Cyp6w1_GLY allele. While our transgenic manipulation of this gene shows that increased expression of three different alleles of can indeed lead to increased resistance to DDT, it is amino acid variants that appear to explain the difference in resistance between DGRP lines. The associated variant is a rare triallelic site where each state encodes a different amino acid (Valine GTG, Alanine GCG, or Glycine GGG). No sites are in LD with it that are also significantly associated with 24hrm, and the sharp and immediate decrease in EHH for both the Val370 and Ala370 variants makes it highly unlikely that this variant is a proxy for an undetected, linked regulatory change polymorphism. Amino acid site 370 of Cyp6w1 is predicted to occur in Substrate Recognition Domain SRS5 (Gotoh 1992; Zawaira ), of an enzyme in a gene family frequently associated with insecticide resistance. It is therefore completely credible that it acts to metabolize DDT, perhaps in a similar way to its paralog Cyp6g1 (Joußen ; Hoi ).

Inspection of CYP6W1 site 370 in related species suggests that Val370 is the ancestral state of this site; the Ala370 and Gly370 substitutions are therefore both derived. The Ala370 variant is strongly associated (P = 5 × 10−10) with reduced 24hrm in our GWAS, and this is supported by our transgenic work, which found overexpression of Cyp6w1_ALA increased DDT resistance relative to Cyp6w1_VAL and Cyp6w1_GLY. Calculation of iHS and nSL did identify a significant signal of selection at Cyp6w1_ALA, but it was associated with the other derived allele, Gly370. This allele is not significantly associated with 24-hrm in our GWAS, and transgenic overexpression of Cyp6w1_GLY does confer resistance to DDT relative to Cyp6w1_VAL. There is no evidence, therefore, that DDT is responsible for the selective footprint identified at , and we propose that the Gly370 mutation increases fitness under a different selective pressure, possibly by increasing the affinity of to another insecticide. The extreme population differentiation observed (Table 1) is also consistent with multiple selective agents acting on this site.

GWAS-associated variants in 91-R and 91-C

Interestingly, resistance-associated variants in and from our GWAS were present in 91-R but not 91-C. To see if any of our other DGRP resistance candidates were differentiated between 91-R and 91-C, we manually searched the genome sequences of these lines generated by Steele for each of our GWAS candidates (P < 1 × 10−5; File S2). Fifty distinct GWAS candidate variants varied between 91-R and 91-C. Of these, 17 (including the two SNPs in and the triallele) varied in the expected direction, with the “resistant” allele in 91-R and the “susceptible” allele in 91-C (Table 2). Thus, these variants may have been segregating in the population used to found 91-R and 91-C and potentially contributed to the strong selection response observed in 91-R.

Table 2

GWAS resistance variants enriched in 91-R but not 91-C

Phenotype	Pipeline	Location [V6]	P-value	Gene	Site class	91-C	91-R
24 hrm	PLINK	2R:6,174,944	1.72E−09	Cyp6w1	Nonsynonymous	T	C
4 hrkd	DGRP2, PLINK	2R:19,169,429	4.28E−07, 7.58E−07	CG10737	Intron	T	G
4 hrkd	DGRP2	2L:14,188,350	1.64E−06	smi35A	Intron	G	G/C
4 hrkd	PLINK	3L:4,059,748	1.75E−06	—	Intergenic	T	C
4 hrkd	PLINK	2L:5,600,056	1.97E−06	—	Intergenic	C	A/C
4 hrkd	PLINK	3L:12,404,415	3.00E−06	CG32103	Synonymous	C	C/T
4 hrkd	DGRP2	3R:30,968,896	3.70E−06	—	Intergenic	A	ATA
4 hrkd	DGRP2, PLINK	2R:19,169,399	5.45E−06; 7.62E−06	CG10737	Intron	T	C
24 hrm	PLINK	2L:10,546,377	5.68E−06	Trim9	Intron	A	A/G
4 hrkd	DGRP2	3L:7,102,962	6.79E−06	form3	Intron	C	C/T
4 hrkd	PLINK	3L:6,935,976	7.88E−06	—	Intergenic	A	A/G
24 hrm	PLINK	3L:15,854,279	8.76E−06	pHCl	Intron	G	A
4 hrkd	DGRP2	2L:8,232,517	9.49E−06	Pvr	Intron	C/A	C
4 hrkd	DGRP2	2L:4,086,206	9.57E−06	ed	Intron	T	G
4 hrkd	DGRP2	2L:19,231,491	1.66E−05	—	Intergenic	A	G
4 hrkd	DGRP2, PLINK	3L:7,495,781	2.02E−05; 7.38E−06	Cyp4d8	Intron	G	T
4 hrkd	DGRP2	3R:25,716,964	2.15E−05	LpR2	Intron	A	C

Other DDT resistance genes

Three substitutions in increasing its DDT metabolism capacity have been described in RalDDTR (Amichot ), a laboratory strain selected for DDT resistance over 50 generations (Cùany ). The first two substitutions (R335S and L336V) are absent from the DGRP, as well as population samples from Africa Drosophila Population Genomics Project; (DPGP; Pool ), France and Portugal in Europe, Florida and Maine in the USA, and Queensland and Tasmania in Australia (Pool ; Martins ; Reinhardt ; Bergman and Haddrill 2015). The third substitution, V476L, is found in around a quarter of DGRP and DPGP sequences (24.1% of DGRP, 26.5% of screened DGRP, 32.2% of DPGP), but this variant does not affect DDT metabolism by itself (Amichot ), and was not associated with either 4 hrkd or 24 hrm at P < 1 × 10−5. Similarly, the resistance-linked promoter variant described by Wan is present in the DGRP (68.8% of DGRP, 68.6% of screened DGRP), but was not associated with either of our phenotypes at P < 1 × 10−5.

Another P450, Cyp12d1, is overexpressed in the DDT-resistant Wisc-1 strain relative to Canton-S (Pedra ) and has been shown to increase DDT resistance when overexpressed transgenically (Daborn ), as well as reducing resistance when knocked down with RNAi (Gellatly ). We did not find associations with Cyp12d1 in any of our four GWAS, nor did we find strong correlations between our phenotypes and Cyp12d1 copy number variation or transcription level (File S3). Steele likewise found no differentiation in Cyp12d1 coding regions between 91-R and 91-C, and our inspection of these genomes found no differentiation in copy number between the two strains.

Cyp12d1 has been shown to be inducible by DDT (Brandt ; Festucci-Buselli ; Willoughby ), as well as other xenobiotics (Willoughby ; Misra ). Our results do not rule out the involvement of this gene in DDT resistance in the DGRP, but show that genomic variation at Cyp12d1 and its transcription level without induction do not correlate with the DDT resistance phenotypes measured in this study (File S5).

In the case of , neither the kdr nor super-kdr resistance mutations are present in the DGRP, nor any of the analogous sites conferring resistance in mutagenized lines characterized by Pittendrigh or Lindsay .

Conclusions

There are three main lines of evidence used in the literature to associate genes with DDT resistance in D. melanogaster: differentiation between resistant and susceptible laboratory strains, association of natural alleles with resistance in outbred populations, and transgenic manipulation in a controlled genetic background. Until now, only has bridged the divides between these approaches. With this study, and can be added to this list. These results alone support a polygenic model and are reinforced by the possibility that other genes from our GWAS, or other studies, also contribute to DDT resistance but have yet to be verified. We have focused on two specific DDT-related traits; 4hrkd and 24hrm, both in adult females with a single dose of DDT, and they give only partially correlated results. Which of these has greater relevance in the field, however, is not clear. Furthermore, there are many alternate parameter values and indeed parameters (e.g., temperature; Fournier-Level ) that may make laboratory-based DDT assays more accurately reflect the variation most relevant to field survivorship. Perhaps one of the most compelling diagnostics for a relevant assay would be if the loci uncovered showed the molecular signs of selection. Previous studies have revealed that the major DDT resistance locus of some populations, , shows the hallmarks of a selective sweep. We did not, however, observe signs of selection around our DGRP GWAS candidates. While we did detect signs of a sweep at , it was associated with a second derived variant, not with the variant associated with DDT resistance. Given that is also capable of providing resistance to other insecticides, the failure to find sweeps at DDT-specific resistance loci strengthens the possibility that other insecticides have contributed to driving the selective sweeps in field populations of D. melanogaster. Finally, it is also a possibility that statistics such as iHS and nSL are insensitive to the slight changes in allele frequency that might occur if a trait is largely polygenic.

Supplementary Material

Supplemental material is available online at www.genetics.org/lookup/suppl/doi:10.1534/genetics.117.300310/-/DC1.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.