Changes in Neuronal Signaling and Cell Stress Response Pathways are Associated with a Multigenic Response of Drosophila melanogaster to DDT Selection.

The adaptation of insect populations to insecticidal control is a continual threat to human health and sustainable agricultural practices, but many complex genomic mechanisms involved in this adaption remain poorly understood. This study applied a systems approach to investigate the interconnections between structural and functional variance in response to dichlorodiphenyltrichloroethane (DDT) within the Drosophila melanogaster strain 91-R. Directional selection in 6 selective sweeps coincided with constitutive gene expression differences in DDT resistant flies, including the most highly upregulated transcript, Unc-115 b, which plays a central role in axon guidance, and the most highly downregulated transcript, the angiopoietin-like CG31832, which is involved in directing vascular branching and dendrite outgrowth but likely may be under trans-regulatory control. Direct functions and protein-protein interactions mediated by differentially expressed transcripts control changes in cell migration, signal transduction, and gene regulatory cascades that impact the nervous system. Although changes to cellular stress response pathways involve 8 different cytochrome P450s, stress response, and apoptosis is controlled by a multifacetted regulatory mechanism. These data demonstrate that DDT selection in 91-R may have resulted in genome-wide adaptations that impacts genetic and signal transduction pathways that converge to modify stress response, cell survival, and neurological functions. This study implicates the involvement of a multigenic mechanism in the adaptation to a chemical insecticide, which impact interconnected regulatory cascades. We propose that DDT selection within 91-R might act systemically, wherein pathway interactions function to reinforce the epistatic effects of individual adaptive changes on an additive or nonadditive basis.

Introduction

Although adaptation by insect species to survive increasing exposures to insecticidal compounds poses a widespread threat to agriculture and human health (Mallet 1989; Ranson and Lissenden 2016), the systemic or genome-wide impacts of such adaptation are yet not fully understood. The classic paradigm ascribes arthropod resistance to insecticidal neurotoxicants either to the effects of one or perhaps several genes that facilitate the metabolic breakdown of xenobiotic compounds (Ffrench-Constant et al. 2004) or to modification of receptor proteins that inhibit toxin binding and their subsequent neurogenic effects (Ffrench-Constant et al. 1993; Ffrench-Constant 2013).

Indeed, low- to moderate-levels of DDT resistance in Drosophila melanogaster have been associated with single locus mutations in the sodium channel gene, para (Pittendrigh et al. 1997; Martin et al. 2000) in laboratory generated para mutant populations, as well as with metabolic resistance conferred by the upregulation of cytochrome P450s (Brandt et al. 2002; Daborn et al. 2007). In contrast, high levels of DDT resistance in the D. melanogaster strain 91-R is associated with multiple changes in physiology (Strycharz et al. 2013), with expression levels of ABC transporters (Strycharz et al. 2013; Seong et al. 2016), and with genes involved in xenobiotic detoxification, ion transport, transcription, and signal transduction (Pedra et al. 2004; Qiu et al. 2013). The necessity of candidate ABC transporters and P450s for DDT detoxification in stress response were functionally validated subsequently via RNAi knockdown (Gellatly et al. 2015); similarly, expression of several detoxification enzymes are controlled by the Nrf2/Keap2 pathway (Misra et al. 2013).

In general, many complex traits result from selection at multiple independent genetic loci (Najarro et al. 2015), and this likely extends to high levels of DDT resistance in 91-R as well. Previously, structural signatures of adaptation to DDT within the genome of 91-R were shown by near allelic fixation in 13 different genomic regions (Steele et al. 2015), with protein coding regions under the influence of directional selection (Steele et al. 2014). Specifically, the effects of DDT selection on 13 genome intervals on chromosomes 2L, 2R, and 3R in 91-R were inferred by reductions in nucleotide diversity and corresponding estimates of directional selection using Tajima’s D when compared with the susceptible control 91-C. Additional functional analysis of the ABC transporter, MDR49, within a selective sweep of the 91-R genome suggested a likely role in conferring DDT resistance (Seong et al. 2016). Despite these advances, the impact of these structural changes on functional variation between the 91-R and its corresponding DDT-susceptible strain 91-C remain not fully understood, particularly with respect to their relative contributions to DDT resistance levels.

The control of cell survival and the response to stress conditions involves intertwined sensing and response pathways, which include both general (constitutive) and stimulus-specific (inducible) interactions to modulate gene expression (Spriggs et al. 2010) leading to changes in cytoskeletal reorganization (Galbraith et al. 1998; Tzima et al. 2002), the accumulation of stress response proteins (Guertin et al. 2010), and changes within signal transduction pathways (Gehart et al. 2010). These responses often converge to effectively increase resilience to oxidative stress (Wang et al. 2003), regulate cell homeostasis (Slack et al. 2010), maintain the integrity of cellular DNA, lipids, and proteins via molecular chaperones (Tower 2011), and repair enzymes that help to increase fly lifespan (Shaposhnikov et al. 2015). For example, changes in pathways that modulate the expression of genes involved in detoxification and increased lifespan, such as the nuclear hormone receptor (NHR) insulin/insulin-like growth factor-like (IGF) signaling (IIS) pathway, may increase the capacity of an organism to resist cytotoxic effects of xenobiotics (Slack et al. 2011). Where the downstream IIS target the forkhead box-containing protein subfamily O, FOXO, transcription factors modulate the expression of genes involved in cycle control, metabolism, apoptosis, and cellular stress responses (Slack et al. 2011). More recently, the Drosophila NHR, DHR96, was shown to be a crucial factor for increasing resistance to DDT; while often co-occurring with IIS in natural conditions, it is functionally independent of IIS (Afschar et al. 2016). Thus, a multitude of metabolic responses and fixed adaptive mutations have been shown to contribute to DDT resistance in 91-R, but the complex epistatic interactions and modifications to inclusive pathway components in relation to adaptations to xenobiotic resistance remain unknown. Moreover, DDT target site mutations and changes in neurologic function have not been discovered to date in 91-R.

To partially address this knowledge gap, this study applied a systems approach to concurrently assess genome-wide structural and functional variation in the DDT resistant D. melanogaster strain 91-R. Specifically, constitutive alteration of RNA-seq estimated transcript levels in 91-R compared with its corresponding DDT-susceptible 91-C counterpart were assessed in the context of regulatory pathways and any causal adaptive mutations underlying these changes.

Materials and Methods

Topical DDT Bioassay

To confirm resistance in 91-R, we used a previously described bioassay approach with some modifications (Strycharz et al. 2013). Females (1–5 days old, mated) were used as the experimental group to assess adult mortality in 91-R and 91-C using an endpoint measure within a DDT topical application bioassay. Mortality rates were then analyzed within strain by the probit analysis (XLSTAT 2008, Addinsoft). The 50% lethal dose (LD50), 95% confidence intervals, and slopes and intercepts for dose-response curves were generated.

RNA-Seq and Estimation of Differential Gene Expression

The DDT resistant D. melanogaster strain, 91-R, and corresponding non-DDT selected control strain, 91-C, were established previously (Merrell and Underhill 1956; Merrell 1960, 1965) and reared in separate colonies on a commercially available medium (Jazz-Mix Drosophila Food, Fischer Scientific) at 25 °C in plastic bottles with transfer to new bottles about every 2 weeks. The 91-R strain has been continually selected by maintaining the flies in a colony bottle with the presence of a 150-mg DDT/filter paper disk while 91-C was maintained without any exposure to DDT. A total of 1,000 5-day-old D. melanogaster flies (500 male and 500 female), not exposed to DDT within that generation, were pooled for each of three biological replicates per strain. Total RNA was extracted from each pool (n = 6) using the Qiagen RNeasy Maxi Kit according to the manufacturer’s instructions (Qiagen, Valencia, CA). Bidirectional Illumina RNA-seq libraries were constructed from each pool, and 160-bp paired-end (PE) sequence read data generated on an Illumina HiSeq 2500 at the University of Illinois Urbana-Champaign, W. M. Keck Center for Comparative and Functional Genomics (one lane per library; three lanes run for each 91-R and 91-C library). Six sequence runs were submitted to the National Center for Biotechnology Information (NCBI) sequence read archive (SRA; accession numbers SRX2611754-SRX2611759). Raw PE fastq-formatted read data were imported into the CLC Genomics Workbench software version 8.5 (Qiagen), each then individually trimmed to remove Illumina adapter sequences, low Phred quality regions (q < 20), short reads (<50 bp), and long oligo Ts (>50). The Drosophila melanogaster Release 6 reference genome assembly including annotated genes and transcripts was downloaded from Flybase database (http://flybase.org; last accessed December 06, 2017) (Dos Santos et al. 2015) and imported into the CLC Genomics Workbench RNA-seq Analysis Tool where trimmed RNA-seq read data from each library were independently aligned with parameters: minimum length fraction 0.9, minimum similarity fraction 0.8, insertion/deletion cost 3, and mismatch cost 3. Counts for the number of reads that aligned to each gene model were output in tab-delimited format.

Differences in gene expression (transcript abundance) were compared between 91-R and 91-C strains using the number of reads per kilobase per million (RPKM) of RNA-seq reads mapped against the annotated gene models in the Drosophila genome assembly v.6.07. Output was log2 transformed and normalized using the median method (Allison et al. 2006). Statistical significance of any differences in normalized read counts (expression levels) between strains were conducted using the “Exact Test” between two groups within the Differential Gene Expression tool of CLC Genomic Workbench (Robinson and Smyth 2008), as implemented in the EdgeR Bioconductor package (Robinson et al. 2010) to iteratively model dispersion for each gene used as a weighted average for subsequent comparisons (Robinson and Smyth 2007). Significances of any differences in RNA-seq reads mapped to gene models between 91-C and 91-R were corrected through multiple testing (Benjamini and Hochberg 1995). A low-stringency significance threshold of a corrected ≤0.01 FDR (false discovery rate) and log2 (fold change) ≥|1.0|, and a high-stringency threshold of a corrected ≤0.001 FDR and log2 (fold change) ≥|2.0|, were implemented when determining differential transcript expression.

Gene Expression by Reverse Transcriptase-Quantitative PCR

Reverse transcriptase-quantitative PCR (RT-qPCR) was carried out on selected transcripts for validation of RNA-seq estimated differences between the 91-R and 91-C strains. Both the RT-qPCR method and first-strand cDNA synthesis were similar to that previously described by Sun et al. (2015). Potential genomic DNA contamination was removed by treatment with RNase-free DNase (Qiagen). Relative expression was calculated using the 2−ΔΔ method (Livak and Schmittgen 2001) with rp49 as the reference gene (Ponton et al. 2011). Multiple linear regressions were carried out using XLSTAT (Addinsoft) to compare relative expression based on RT-qPCR and RNA-seq for each gene.

To analyze differential gene expression at various developmental stages of the 91-R and 91-C strain, total RNA samples were prepared using RNeasy Mini Kit (Qiagen) from a pool of 20 eggs, an individual larva, an individual pupa, and an individual virgin adult female. At least three independent biological samples were prepared for each. The specific primer pairs for RT-qPCR were used (supplementary table S1, Supplementary Material online) with statistical analysis by XLSTAT (Addinsoft) using a two-way ANOVA followed by Tukey’s multiple comparison tests for significant differences (P <0.05) in gene expression during developmental stages.

Genome Positioning of Differentially Expressed Transcripts

Thirteen regions in the genome of 91-R were previously defined as having undergone a selective sweep that reduced nucleotide diversities following DDT selection for survival (Selective sweeps 1–13; Steele et al. 2015). The genome intervals spanned by genes encoding differentially expressed transcripts were retrieved from Flybase.org and compare with the set of 100,000-bp intervals corresponding to the 13 selective sweeps. The gene symbols for differentially expressed genes (DEGs) within selective sweeps were used to query the collection of curated transcription factor (TF)-gene interactions maintained in the REDFly database (http://redfly.ccr.buffalo.edu; last accessed December 06, 2017) and experimentally supported TF-binding domain retrieved from modENCODE (http://www.modencode.org/; last accessed December 06, 2017). Corresponding upstream regulatory regions within extended gene regions were obtained from FlyBase.org, to which genomic read data were mapped and SNPs detected as described previously (Steele et al. 2014).

Ontologies, Pathways, and Interactions among Differentially Expressed Genes

For all transcripts surpassing the high-stringency threshold, information regarding protein–protein interaction from yeast-two-hybrid data were obtained from Flymine.org (Gene -> Interactions) and defined as first-order protein–protein interactions. Gene descriptions and protein domain information were retrieved from Flybase database (http://flybase.org; last accessed December 06, 2017) and Gene Ontology (GO) terms were retrieved from Flymine database (http://flymine.org; last accessed December 06, 2017) (Gene -> GO term queries). Additionally, the Flymine database (http://flymine.org; last accessed December 06, 2017) of phenotypes resulting from RNAi knockdown experiments (Gene -> RNAi phenotypes) were retrieved and clustered based on shared phenotypic term. The biochemical pathway identifiers in which gene products reside were analogously obtained from Flymine.org (Gene -> pathway) and clustered by shared pathway. Functional interaction evidence (text mining, experimental, database, and coexpression data) was obtained from string-db version 10 (http://version10.string-db.org/; last accessed December 06, 2017) (Szklarczyk et al. 2015) for all derived gene products that surpassed the high-stringency threshold.

All of the above searches were analogously performed as well for yeast-2-hybrid partners for each DEG identified from interaction data (http://flymine.org/; last accessed December 06, 2017); Gene -> Interactions). Relevant pathway descriptions and pathways diagrams were obtained from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg/pathway.html; last accessed December 06, 2017) Interactions among GO terms and RNAi phenotypes shared among yeast-2-hybrid partners for DEGs were visualized using Cytoscape 3.4.0 (Shannon et al. 2003), defining GO or RNAi phenotypes as the Target, DEG, and the Source, and yeast-2-hybrid partners as the Interaction. Regulation of each DEG (up- or down-) was used as a Source Attribute and mapped onto resulting networks.

Results

Mortalities among 91-R adults after topical exposures of DDT were significantly lower compared with 91-C when measured at a 24-h bioassay end-point. Specifically, the estimated LD50 values across replicates of 91-R had a mean of 3,746.8 ng (2,577.2–5,202.2 ng, 95% CL) per individual versus 35.0 ng (31.9–38.4 ng, 95% CL) per individual for 91-C, and a resultant >107.1-fold resistance to DDT among 91-R compared with 91-C (supplementary table S2, Supplementary Material online).

Overall, ≥82.68% RNA-seq read data from replicates of 91-R (n = 3) and 91-C (n = 3; supplementary table S3, Supplementary Material online), from which 14,112,989 (8.9%) of the reads were trimmed, were subsequently aligned to the D. melanogaster genome v.6.07. Estimates of differential expression based on the normalized counts of RNA-seq read data that mapped to gene models predicted a total of 137 differentially expressed transcripts at the low-stringency, with 79 up- and 58 downregulated in 91-R compared with 91-C (supplementary table S4, Supplementary Material online) and; 69 transcripts for the high-stringency, with 42 up- and 27 downregulated (tables 1 and 2, respectively). Estimates of fold-change between 11 transcripts from 91-R and 91-C were highly correlated in the RT-qPCR and RNA-seq methods (R2 = 0.88, F = 102.2, P < 0.001; supplementary fig. S1, Supplementary Material online).

Table 1

Annotations for Upregulated Transcripts in the DDT Resistant Strain 91-R

Gene Symbol	Functional Annotation	Log2 Fold Changea	FDRb
Unc-115b	Uncoordinated 115b	5.89	3.31E-06
Cyp4p2	Cytochrome P450-4p2	5.64	5.03E-04
Cyp6a8	Cytochrome P450-6a8	5.59	7.72E-92
Cyp4p1	Cytochrome P450-4p1	5.02	1.19E-27
Cyp6g1	Cytochrome P450-6g1	4.08	1.14E-23
Cyp6w1	Cytochrome P450-6w1	3.56	6.73E-94
Cyp6g2	Cytochrome P450-6g2	2.33	7.06E-03
CG14820	Peptidase_M14	3.82	7.05E-07
Qtzl	Quetzalcoatl	3.71	2.21E-04
CG7542	Trypsin	3.65	7.89E-45
CG31683	Unknown	3.58	4.30E-18
CG4250	Zn finger protein	3.54	6.20E-45
Jon65Aii	Jonah 65Aii	3.53	6.75E-23
Jon99Ci	Jonah 99Cii	2.28	1.10E-09
CG18787	Nucleoporin-like protein 2	3.52	1.23E-07
Mal-B1	Maltase-B1	3.44	7.64E-03
CG13658	Alpha-amylase	3.28	2.37E-06
CG5770	Peptidase_C11	3.04	3.71E-06
CG3819	Endonuclease_NS	2.96	4.55E-09
Npc2e	Niemann-Pick type C-2e	2.93	6.73E-03
CG33301	EcKinase	2.91	4.48E-17
CG43055	C-lectin domains-containing	2.86	1.61E-09
CG31703	RNA polymerase-associated protein RTF1	2.74	3.03E-03
Spn43Aa	Serpin 43Aa	2.73	6.17E-03
CG4210	Histone N-acetyltransferases	2.54	1.60E-04
Sfp24Ba	Seminal fluid protein 24Ba	2.50	1.76E-14
CG6908	EcKinase	2.20	1.99E-20
Lectin-galC1	Galactose-specific C-type lectin	2.16	5.23E-15
Sodh-2	Sorbitol dehydrogenase-2	2.14	2.85E-03
CG31809	Alcohol dehydrogenase	2.06	1.18E-04
CR33319	Unknown	4.84	2.16E-19
CR42653	Unknown	4.72	5.03E-04
CR44030	Unknown	4.65	1.80E-03
CR44308	Unknown	4.64	1.25E-04
CG33966	Unknown	4.34	1.80E-03
CG42521	Unknown	4.18	3.72E-05
CG31683	Unknown	3.58	4.30E-18
CR45191	Unknown	3.24	2.30E-04
CG14191	Unknown	3.23	3.91E-03
CR43239	Unknown	2.74	6.08E-07
CG6834	Unknown	2.71	7.23E-05
CG15293	Unknown	2.35	3.64E-66
CG33666	Unknown	2.25	2.98E-05

Fold change was calculated as log2 91-C/91-R.

FDR: False discovery rate. Differentially expressed genes were identified at the high-stringency threshold [FDR < 0.001 and log2(fold change) ≥|2.0|] of 91-C/91-R.

Table 2

Annotations for Downregulated Transcripts in the DDT Resistant Strain 91-R

Gene	Functional Annotation	Log2 Fold Changea	FDRb
CG31832	Tyrosine kinase activator and a negative regulator of apoptosis	−6.76	7.33E-14
Ser12	Serine protease 12	−6.62	7.87E-09
CG30025	Serine protease	−5.13	2.79E-13
Hsp70Ba	Heat shock protein 70	−5.11	8.93E-03
PPO1	Prophenoloxidase 1	−4.35	1.17E-07
CG17234	Serine protease	−4.48	1.87E-03
CG31274	Leucine-rich_repeat	−3.57	1.23E-07
CG4927	Serine protease	−3.47	4.33E-03
CG30054	Guanine nucleotide-binding protein G(q) subunit alpha	−3.42	1.31E-05
PPO2	Prophenoloxidase 2	−3.31	1.87E-07
Hsp70Bb	Heat shock protein 70Bb	−3.28	8.56E-03
Jon99Fii	Jonah 99Fii	−3.24	1.05E-91
Cyp4e3	Cytochrome P450-4e3	−3.19	6.72E-18
CG3397	Aldo/keto reductaseNADP-dependent oxidoreductase	−3.06	2.93E-08
TotA	Turandot A	−2.83	1.44E-30
CG4563	AMP-dependent synthetase/ligase	−2.67	9.26E-04
CG2064	Alcohol dehydrogenase	−2.58	5.74E-12
CG31091	Abhydro lipase	−2.36	4.29E-04
TotC	Turandot C	−2.33	4.17E-11
CG43061	Pathogenesis	−2.31	4.58E-06
Idgf1	Imaginal disc growth factor 1	−2.26	7.47E-07
CG4757	Acetylcholinesterase	−2.19	7.78E-08
GILT2	Gamma interferon inducible lysosomal thiol reductase 2	−2.10	7.83E-33
CG6277	Lipase	−2.05	1.19E-11
CR33487	Unknown	−5.51	1.18E-05
CR43263	Unknown	−5.31	6.28E-03
CR43158	Unknown	−5.16	6.24E-06

Fold change was calculated as log2 91-C/91-R.

FDR: False discovery rate. Differentially expressed genes were identified at the high-stringency threshold [FDR < 0.001 and log2(fold change) ≥|2.0|] of 91-C/91-R.

Stage-Specific Changes in Gene Expression

RT-qPCR measurements of a subset of transcripts differentially expressed across 91-R and 91-C egg, larva, pupa, and adult stages indicated that changes are mostly constitutive across those growth stages. Results showed that uncoordinated 115b (Unc-115 b) and jonah (Jon) 65Aii were constitutively more highly expressed in 91-R compared with 91-C across all stages (fig. 1), whereas imaginal disc growth factor 1 (Idgf1) was not differentially regulated in any stage (fig. 2). Several transcripts showed greater stage specificity. For example, the C-type lectin domain-encoding CG43055 was predominantly expressed in the adult stage for both strains, but 91-R had a ∼15-fold higher level compared with 91-C. In contrast, expression of the other C-type lectin, Lectin-galC1, occurred during larval and pupal stages (fig. 1). Similarly, Jon99Cii was predominantly expressed in the 91-R larval and adult stages (fig. 1).

. 1.

—Relative expression levels of several upregulated genes in different developmental stages of 91-R and 91-C strains at the high-stringency threshold [FDR < 0.001 and log2(fold change) ≥|2.0|], determined by RT-qPCR. The mRNA level in different developmental stages was normalized using rp49 as a reference gene. The mean expression in each developmental stage is shown as a fold change compared with the lowest mean expression, which has been ascribed an arbitrary value of 1. The vertical bars indicate standard error of the mean (SEM) (n = 3). Different letters on the bars indicate that the means are significantly different throughout the different developmental stages.

. 2.

—Relative expression levels of several downregulated genes in different developmental stages of 91-R and 91-C strain at the high-stringency threshold [FDR < 0.001 and log2(fold change) ≥|2.0|], determined by RT-qPCR. The mRNA level in different developmental stages was normalized using rp49 as a reference gene. The mean expression in each developmental stage is shown as a fold change compared with the lowest mean expression, which has been ascribed an arbitrary value of 1. The vertical bars indicate standard error of the mean (SEM) (n = 3). Different letters on the bars indicate that the means are significantly different throughout the different development stages.

Interestingly, the serine protease inhibitor, Spn43Aa, was expressed predominantly in the pupal stage for both strains, but at ∼4.02-fold higher levels in adults from 91-R compared with 91-C (fig. 1). Analogously, while Jon99Fii was highly expressed in larval stage for both strains, the transcript was ∼1.9-fold more abundant in 91-R larvae (fig. 2). Conversely, while the prophenoloxidases, PPO1 and PPO2, were predominantly expressed during the larval stage for both strains, 91-C showed a ∼1.6- and ∼2.25-fold higher expression as compared with 91-R (fig. 2). In addition, changes in differential regulation between growth stages were shown for maltase-B1 (MalB1), with a higher expression during the pupal stage for 91-C compared with 91-R (fig. 1). Lastly, while of approximately equal expression level for 91-R and 91-C adults alike, the niemann-pick type C-2e (Npc2e) peptidoglycan recognition protein was the only transcript that partially conflicted with RNA-seq. In contrast, Npc2e was most highly expressed in 91-C at the larval and pupal stages. Only Unc-115 b, Jon65Aii, and Idgf1 were expressed at high levels in eggs.

About 6 of the 69 DEGs based on RNA-seq analyses were colocated with genome regions of 91-R previously placed within selective sweeps (e.g., reduced nucleotide variation; Steele et al. 2015). Four of these transcripts were upregulated in 91-R, including the most highly upregulated transcript, Unc-115 b, located within sweep 13 on chromosome 3R (table 3). The three remaining upregulated transcripts, CG31683, CR43239, and Lectin-galC1, resided within a successive series of regions with reduced nucleotide diversity along 2L. Two transcripts with reduced relative expression in 91-R, CR43263, and PPO2, were colocated with selective sweeps on 2L (sweep 1) and 2R (sweep 11), respectively. This is based on a question if the 6 DEGs located within 178,075 bp (size of the 13 selective sweeps; Steele et al. 2014; P = 3.369×10−5) is more rare of an event as compared with 69 total DEGs being randomly distributed across the 143.9 Mb of the entire genome (4.795×10−7). To do this, we calculated the Z-score for two population proportions, which provides Z-score of 0.0441 and a nonsignificant P value of 0.48405.

Table 3

Differentially Expressed Transcripts Encoded within 91-R Genome Regions Putatively under the Influence of Selective Sweeps in Response to DDT Selection

Gene	Functional Annotation	Log2 Fold Changea	FDRb	Genome Position	Selective Sweep IDc
Unc-115b	Cytoskeletal linker protein that acts in axon guidance	5.89	3.31E-06	3R: 9,687,060–9,692,775	13
CG31683	Lipid metabolic process	3.58	4.30E-18	2L: 20,447,734–20,449,964	6
CR43239	Unknown	2.74	6.08E-07	2L: 17,353,730–17,354,668	4
Lectin-galC1	Immune pathway	2.16	5.23E-15	2L: 19,417,447–19,418,274	5
CR43263	Unknown	−5.31	0.006	2L: 1,216,854–1,217,785	1
PPO2	Immune pathway	−3.31	1.87E-07	2R: 9,042,261–9,044,709	11

Fold change was calculated as log2 91-C/91-R.

FDR: False discovery rate. Differentially expressed genes were identified on the basis of FDR < 0.001.

Thirteen genomic intervals in Drosophila melanogaster strain 91-R identified regions under the influence of a selective sweep (ID 1–13), where regions show reduced nucleotide diversity in 91-R compared with 91-C and putatively under the influence of directional selection in response to differential survivorship when exposed to DDT (Steele et al. 2015).

Promoters of 6 differentially expressed genes have 36 previously identified TF-gene interactions involving 25 different TFs (droidb.org database Feb 2017; 1.57 ± 0.79 promoters per TF). Four of these TFs, chromato (Chro), centrosomal protein 190kD (Cp190), CTCF, and trithorax (trx), all bind at upregulated CG31683, Lectin-galC, and Unc-115 b (supplementary table S5, Supplementary Material online). Protein–protein interactions between boundary element-associated factor of 32kD (BEAF-32), Chro, and Cp190 are supported by yeast-2-hybrid experiments reported at Flybase.org (not shown). Furthermore, these four TFs, along with senseless (sens) and BEAF-32, bind the promoters of both Unc-115 b and its nearest paralog, Unc-115ba (supplementary table S5, Supplementary Material online). The promoter binding domains in modENCODE provide evidence that these TFs are mostly positioned upstream of Unc-115 b splice variants, except for variant Unc-115 b-RE (fig. 3).

. 3.

—Splice variants for uncoordinated-115b (unc-115 b; reverse strand of 3R) and corresponding putative upstream transcription factor binding domains; cAMP responsive element binding protein 2 (CRE-BP1), DNA polymerase II (DNA Pol II), TATA binding protein (TBP), BEAF-32 (boundary element-associated factor-32 kDa), CP190 (centrosomal protein 190 kDa), CTCF, Chro (Chromator), Chriz, and Sens (senseless). Flanking CG32939 defines the putative gene boundary.

Mapping of pooled genomic-seq reads (Steele et al. 2014) identified three mutation sites that were within CHiP-seq defined TF binding regions, and these mutations were also in or flanking consensus TF binding domains upstream of Unc-115 b. Specifically, the cAMP responsive element binding protein 1 (CRE-BP1) consensus binding sequence contained an indel (CTC insertion) and a substitution mutation within different reads from 91-R. None of these mutations were present among reads from 91-C. Additionally, a substitution mutation flanking a BEAF-32 recognition sequence (CGATA) was at higher frequency in 91-R (0.60) compared with 91-C (0.08).

A total of 230 terms for molecular function (F; n = 97; 40 up- and 57 downregulated), biological process (P; n = 92; 59 up- and 33 downregulated), and cellular component categories (C; n = 41; 22 up- and 19 downregulated; supplementary fig. S2, Supplementary Material online) were assigned to DEGs. GO category C revealed extracellular region/space was highly represented among gene products from both up- and downregulated DEGs, and the GO category P showed oxidative–reductive processes were most prevalent among upregulated transcripts. GO category F showed heme and iron binding, as well as electron carrier activity, were most represented among upregulated transcripts, whereas serine-type endopeptidase activity was prevalent among downregulated transcripts (supplementary table S6, Supplementary Material online).

Reduction–oxidation (redox) activity was predicted among six upregulated cytochrome P450 monooxygenases (P450s; GO: 0016705) and sorbitol dehydrogenase 2 (sodh-2; GO: 0055114) and analogously among downregulated gamma-interferon-inducible lysosomal thiol reductase 2, GILT2, and PPO 1 and 2. The second-most highly downregulated gene, serine protease 12 (Ser12), is annotated as an extracellular serine-like protease (GO: 0004185) involved in proteolysis (GO: 0006508) and wound healing (GO: 0042060). Downregulated transcripts encoding serine proteases (GO: 0004252) include CG17234, CG30025, CG4927, and Jon99Fii. Of these, CG17234 is involved in sensory perception of pain (GO: 0019233).

In addition to downregulated PPO1 and 2, turandot A and C (TotA and TotC), heat shock proteins (Hsp70a and Hsp70b), and CG43061 respond to cellular stress or pathogen attack. Furthermore, DEGs may be involved in binding to actin (Unc-115 b), myosin (CG6834), galactose (Lectin-galC1 and CG43055), and lipids (CG31683 and Npc2e). The most highly upregulated transcript, Unc-115 b, has actin-binding function (GO: 0003779) and involvement in cytoskeletal organization (GO: 0007010). The most highly downregulated transcript, CG31832, is predicted to function as a transmembrane tyrosine kinase receptor (GO: 0030297) that negatively regulates apoptosis (GO: 0043066).

About 6 (14.3%) and 5 (11.9%) of 42 upregulated transcripts in 91-R, respectively, encode cytochrome P450s and serine protease/serine protease-inhibitor-like domains (supplementary table S7, Supplementary Material online). Protein kinase domains are present in CG6908, CG13658, CG6834, and CG33301; lectin C-type domains were predicted in derived protein products for Lectin-galC1 and CG43055. The most highly upregulated Unc-115 b transcript predicts abLIM-type and villin headpiece domains.

Among the 34 pathways assigned to 19 of the 69 DEGs, 4—21.1%, the downregulated Jon99Fii and upregulated Jon99Ci, Jon65Aii, and CG7542—are involved in the activation of extracellular matrix metalloproteinases (MMPs; supplementary table S8, Supplementary Material online). All other pathways were represented by ≤2 protein coding genes. Five transcripts are nonprotein coding and received no GO or pathway annotations. Three protein coding genes (CG14191, CG33666, and CR45191) are of unknown function.

Protein–protein interactions derived from prior two-hybrid screens within Flymine identified 254 experimental lines of evidence for 201 unique interactions involving 33 DEGs (supplementary table S9, Supplementary Material online). High connectivity was shown via network analyses between downregulated PPO1 and PPO2 and upregulated CG14820 through four different interactions, with CG14820 interacting with the phosphatidylinositol-4-phosphate 5 kinase (PIP5K), skittles (sktl; supplementary fig. S3, Supplementary Material online). The downregulated Hsp70Ba and Hsp70Bb share 22 partners including the Jun N-terminal kinase (bsk). Additionally, the most highly upregulated transcript in 91-R, Unc-115 b, physically interacts with the DNA mismatch repair proteins, Mlh1 and Pms2, and the tetratricopeptide repeat domain-containing Pex5.

While no known interactions have been reported for the most highly downregulated transcript, CG31832, the second most highly downregulated Ser12 interacts with thioredoxin peroxidase 2 (Jafrac2), with yeast-2-hybrid partner interactions showing an abundance of GOs related to stress response and control of apoptosis (fig. 4). Analogously, 8 up- and 6 downregulated transcripts (20.3% of DEGs) have protein network interactions that functionally impact neuron development or function (supplementary fig. S4, Supplementary Material online). Protein–protein interactions for 12 of 69 DEGs (17.4%) tether them to the cytoskeleton or extracellular matrix (supplementary fig. S5, Supplementary Material online).

. 4.

—Functional interaction network connecting the downregulated serine protease, Ser12, to cellular pathways involved in apoptosis and response to stress. Protein–protein interaction and activation of the thioredoxin peroxidase 2, Jafrac2, by Ser12 directly signals the activation of apoptotic cysteine-type endopeptidase, as well as potential regulatory cascades involving prophenol oxidase 2 (PPO2). Jafrac2 is also involved in cellular response to oxidative stress conditions, and regulated cell survival via the Janus activated kinase/signal transducer and activator of transcription (JAK/STAT) pathway. Up- and downregulated transcripts were indicated in green and orange boxes, respectively. Putative networks defined between DEGs and gene ontologies (biological processes; gray nodes) were connected by DEG yeast-2-hybrid protein–protein interaction partners.

Flymine Gene ->RNAi experiments identified 143 and 1,338 previously defined phenotypes, respectively, associations with the knockdown of 40 DEGs and 160 of their 205 yeast-2-hybrid partners (supplementary table S10, Supplementary Material online). Among the DEGs and partners, 44 (30.8%) and 354 (26.5%), respectively, were lethal (supplementary fig. S6, Supplementary Material online). Among the remaining phenotypes associated with DEGs, changes in the phosphorylation state of extracellular signal-regulated kinase (ERK) was most prevalent (n = 12 phenotypes) following knockdown of nine different transcripts (CG4250, CG4563, CG4927, CG7542, Cyp4e3, Hsp70Ba, Idgf1, Mal-B1, and Ser12). This was followed by 11 RNAi phenotypes that resulted in alteration of wing development and flight capacity through the knockdown of CG5770, CG6277, CG6834, Cyp6g1, CG6908, CG14191, and CG18787. In addition, bristle morphology or development was impacted in eight different cases when seven different transcripts were knocked down (CG14191, CG4757, CG5770, CG6277, CG6834, Sodh-2, and TotA). Other phenotypes involved changes in Ca2+ transport and altered pathogen defenses, heat avoidance behaviors, and neural defects.

Discussion

In this study, we identified constitutive differentially expressed genes in DDT resistant 91-R that are associated with neural cell development, signal transduction, and gene regulatory cascades that impact the nervous system. Previously, significant differences in expression were documented for several cytochrome P450s (Cyp6a2, Cyp6g1, and Cyp9c1) between DDT resistant 91-R and the susceptible strain Canton-S (Qiu et al. 2013), as well as Cyp6g1 between 91-R and Wisconsin strain (Pedra et al. 2004). Moreover, use of UAS-RNAi transgenic lines of D. melanogaster has validated the role of several cuticular proteins (Cyp4g1 and Lcp1) and membrane transport genes (e.g., ABC transporters; MDR50, MDR65, and MRP1) involved in DDT resistance for 91-R strain (Strycharz et al. 2013). Despite previous studies pointing out MDR50, MDR65, and MRP1 or Cyp4g1 were more highly expressed in 91-R strain as compared with susceptible strain Canton-S (Gellatly et al. 2015), those genes were not differentially transcribed in the comparison of 91-R with 91-C. The fact that 91-R is 1,500-fold more resistant to DDT when compared with Canton-S (Strycharz et al. 2013) suggests that the possibility that 91-C was already slightly resistant to DDT prior to the laboratory selection of 91-R strain. This possibility is supported by previous data indicating that 91-C strain is slightly more tolerant to DDT (7.5-fold) compared with the Canton-S strain (Festucci-Buselli et al. 2005).

Adaptations in the Migration of Neural Cells

Transcripts for uncoordinated-115b, Unc-115 b, and angiopoietin-like CG31832 are, respectively, the most constitutively up- and downregulated transcripts from the DDT resistant strain 91-R, with both functioning to direct neuron cell movements. Specifically, Unc-115 b is an actin-binding and LIM domain containing protein orthologous to human abLIM (Lundquist et al. 1998). It comprises part of the netrin signaling pathway that conveys guidance signals to advancing distal dendritic growth cones and determines direction of finger-like filopodia and laminopodia outgrowths into the extracellular matrix (ECM; Yang and Lundquist 2005). These nerve outgrowths involve the coordinated modulation of underlying cytoskeletal structures in response to environmental stimuli through redundant Rac GTPases MIG-2/RhoG and CED-10/Rac1(Lundquist et al. 2001), with Unc-115 b as a downstream effector anchored to cytoskeletal F-actin (Struckhoff and Lundquist 2003). Opposing cascades of cone protrusion and retraction are initiated through interaction of cell surface receptors UNC-5 and UNC-40 (Frazzled, fra, or Deletion in Colorectal Cancer, DCC) with the ECM-embedded laminin-like protein, UNC-6 (Netrin-A or -B proteins; Chan et al. 1996; Leonardo et al. 1997). Coordinated movements toward UNC-6 are mediated by UNC-40 homodimers, whereas repulsive signals from UNC-6 at growth cones are mediated by interaction with UNC-5: UNC-40 heterodimers (Hong et al. 1999; MacNeil et al. 2009; Norris and Lundquist 2011). The intracellular repulsive signaling cascades mediated by UNC-5: UNC-40 involve the activation of CED-10/Rac1 and MIG-2/RhoG via the guanidine nucleotide exchange factor (GEF), UNC-73/Trio, at the UNC-40 scaffold (Forsthoefel et al. 2005; Norris et al. 2014). In contrast, the cytoplasmic domains (P1 to P3) of UNC-40 homodimers transduce attractive intracellular signals in response to UNC-6 by interactions between Rho-GTPase CDC-42 and the T-cell lymphoma Invasion and Metastasis Factor 1 (TIAM1)/Drosophila Still Life (sif) GEF (Demarco et al. 2012).

The Unc-115 b CDS is within selective sweep 13 (Steele et al. 2015). The nearly 6-fold increase in transcript levels in 91-R may be influenced by cis-regulatory mutations in, or flanking, TF binding sequences (fig. 3). The exact role of Unc-115b/abLIM in DDT resistance remains unclear. Regardless, functional data suggest that increased cellular concentrations of Unc-115 b could enhance the growth, stability, or signaling to actin-based dendritic extensions. In addition to modulation of cytoskeletal dynamics, abLIMs participate in intracellular signal transduction events and impact transcription (Barrientos et al. 2007), suggesting that upregulation of Unc-115 b may regulate transcription factor activity and be a causal factor that influences global gene expression changes predicted between 91-R and 91-C.

Annotations for the most highly downregulated transcript, CG31832, indicate functioning as a transmembrane receptor protein tyrosine kinase activator and negative regulator of apoptosis. The encoded fibrinogen C-terminal domain suggests possible roles in blood clotting and platelet (cell) aggregation as shown in its nearest ortholog, the agonistic human vascular endothelial growth factor (VEGF) family member angiopoietin 4 (ANGPT4) involved in angiogenesis (Valenzuela et al. 1999). The antagonistic D. melanogaster angiopoietin, Pvf1, inhibits apoptosis in stem cell lines (Xing et al. 2015), whereas the bovine angiopoietin ortholog, Ang-1, also inhibits apoptosis and promotes cell survival in epithelial cells by activation of the Akt pathway (Papapetropoulos et al. 2000); Ang-2 promotes apoptosis via inhibition of Ang-1 in mouse glomeruli (Davis et al. 2007).

Angiopoietins cross-talk with Janus kinase/signal transducers and activators of transcription (JAK/STAT) and mitogen-activated protein kinase (MAPK) pathways to promote cell survival (Rychli et al. 2010). Interestingly, Netrins might regulate aspects of tumor angiogenesis (Klagsbrun and Eichmann 2005), such that analogies between guidance mechanisms in axon and capillary network formation have been made since these involve many of the same receptors and signaling ligands (Autiero et al. 2005; Freitas et al. 2008). The downregulation of the angiopoietin CG31832 in 91-R suggests, analogous to Unc-115 b, that DDT resistance in 91-R might involve changes in neurological cell growth and maintenance. Although angiopoietin orthologs can be agonistic or antagonistic with regard to cell outgrowth and apoptosis, the downregulation of CG31832 might be predicted to enhance axon growth or neuronal resistance to apoptosis, but additional experiments are required to investigate this.

The ECM is composed of a complex matrix of proteins and glycans involved in structural maintenance and cell signaling and is a crucial determinant of neural growth, development, and function (Broadie et al. 2011). DEGs acting within the cell migration process are observed in the upregulation of the serine proteases Jon99Ci, Jon65Aii, and CG7542, which function to cleave leader peptides from zymogens of MMPs. This activation of MMPs, likely by c-Jun N-terminal kinases (JNKs) or phosphatidylinositol 3 kinase (PI3K; Cheng et al. 2012), leads to the degradation of ECM components and subsequently allows for neuronal cell migration and extension of lamellipodia (Zimmermann and Dours-Zimmermann 2008), but also could modulate various growth factor signals from the ECM and thereby have a role in a range of cell processes including immune response and neuron outgrowth (Loffek et al. 2011). Additionally, MMP expression is directed by various growth factors and cytokines (Mauviel 1993) which also lead to cell responses via activation of intracellular signals from Janus protein-tyrosine kinases/transducers and activators of transcription (JAKs/STATs) classes of receptors (Ihle 1996), suggesting the coordination between ECM modulation and cell migration.

Adaptations in Neuron Cell Function

Neurotoxic effects of acute DDT exposure include characteristic twitches induced by an inhibition of the sodium channels para and DSC1 gate function (Rinkevich et al. 2015), the former of which is required for Drosophila mating behavior (Lilly et al. 1994). Furthermore, the chronic systemic effects of DDT and its metabolites are known also to impact brain development and function (Iwaniuk et al. 2006), to cause neurological defects (Grandjean and Landrigan 2014), and to influence the onset of neurological disorders (Li et al. 2015). The significantly downregulated transcript CG3397 in 91-R has previously been shown to be upregulated within the brain of four fly lines with period (per), clock (Clk), timeless (tim), and cycle (cyc) null mutations (Lin et al. 2002), suggesting it may possibly play a role in brain function. Alternatively, since yeast-2-hybrid evidence indicates that CG3397 interacts with thioredoxin-like CG11790, which may in turn mediate redox balance in response to DDT detoxification induced oxidative stress (Li et al. 2008), but further stuides are required to further investigate the possible impacts on DDT resistance in 91-R.

The downregulated serine protease CG17234 and lipase CG6277 are involved in the sensory perception of pain (Neely et al. 2010) and neurogenesis. The CG6277 encoded phospholipase A1 (PLA1)-like domain and phosphatidylcholine 1-acylhydrolase activity (GO: 0008970) catalyze release of fatty acid groups from the 1-position of phosphatidylcholine, while RNAi-mediated knockdown causes a significant reduction in neuroblast proliferation (Neumuller et al. 2011). PLA1 has limited involvement in lysophosphatidylcholine production but is a potent regulator of phosphatidic acid levels in the testes and brain (Higgs and Glomset 1994).

Phosphatidic acid (PA) is an intermediate in the biosynthesis of phospholipids and glycolipids and is implicated in cell signaling through the binding or sequestering of proteins to membranes. For example, PA facilitates Ras activation by recruiting Son of Sevenless (SOS) to plasma membranes (Zhao et al. 2007), and similarly recruits Raf-1, which in turn induces the mitogen activated protein kinase/extracellular signal-regulating kinase (MAPK/ERK) pathway (Rizzo et al. 2000). Additionally, PA is required for activity by the homeotic sensor mammalian target of rapamycin (mTOR), which regulates progression into the cell division cycle (Foster 2013). Thus, cellular levels of PA are crucial for regulating cell survival (growth and migration) and apoptosis along with modulation of cell stress response and cytoskeletal organization (Wang et al. 2006). While the possible role of CG6277 in DDT resistance remains undetermined, significant reduction in PLA1 activity in 91-R may reduce PA levels in the brain, thereby attenuating or modulating the dynamics of cell growth, maintenance, and apoptosis.

Interestingly, along with CG17234, the knockdown of the downregulated extracellular serine protease CG30025 has been shown to increase the sensory perception of pain from heat exposure (Neely et al. 2010). Protease-activated receptors (PARs) are G protein coupled receptors (GPCRs) expressed in tissues, including neurons, that function in the transduction of cytoplasmic signaling from extracellular proteases (Russo et al. 2009b). This occurs via specific protease-dependent N-terminal cleavages that irreversibly activate PARs (Finigan et al. 2005; Russo et al. 2009a) and thus evoke different intracellular responses, including not only the transmission of pain sensation (Bunnett 2006) but possibly also the mediation of inflammatory neuronal cell death (pyroptosis; Li et al. 2016). These lines of evidence suggest that repression of CG17234 and CG30025 might be advantageous in 91-R by reducing PAR activation and thereby hindering signal transduction of signals involved in apoptosis or the perception of pain and inflammation.

Our analysis indicates that the protein products of 8 up- and 6 downregulated transcripts in 91-R are involved in known protein–protein interactions, where the partner proteins function in neural development or function (supplementary fig. S4, Supplementary Material online). Within these, the upregulated extracellular localized CG15293 and downregulated endomembrane system localized reductase CG2064 were the most highly centralized with ≥4 protein interactions that impact ≥9 neurological functions. Little is known regarding the function of these two genes, but transcription of CG15293 is increased in response to sigma virus infection (Carpenter et al. 2009), and CG2064 is upregulated in responses to phospholipase C2 knockdown-mediated cellular repression of ERK response to oxidative stress (Liu et al. 2007).

Additionally, physical interaction of the downregulated CG2064 with Rab2 and Rab5 GTPases—which mediate intracellular vesicle transport (Stenmark 2009), as well as the clathrin heavy chain (Chc), Golgi vesicle transport protein, and small synaptic vesicle protein, synaptobrevin (Syb)—might be hypothesized to impact neurotransmitter release. This suggests, consonant with previous work, that CG2064 is likely intertwined within neural function and cellular stress response pathways (Liu et al. 2007).

Changes in Cell Stress Response

Cellular reactions to stress conditions involve a delicate balance between cell survival (proliferation) and programmed death (apoptotic) pathways. Redox reactions are essential for maintaining cell homeostasis by modulating electron-transfer reactions that occur during normal metabolic processes and in response to environmental stress.

91-R exhibits constitutive upregulation of 6 cytochrome P450 monooxygenases, with only one, Cyp4e3, significantly downregulated. While Drosophila Cyp6g1, Cyp6g2, Cyp6t3, Cyp6a2, Cyp6a8, Cyp6a19, Cyp6w1, and Cyp6a23 have been validated as DDT or imidacloprid metabolizers (Daborn et al. 2007; Morra et al. 2010; Harrop et al. 2014; Le Goff and Hilliou 2017), Cyp4p2 and Cyp4e3 have been found to be highly upregulated in Drosophila resistant to imidacloprid (Kalajdzic et al. 2012) and permethrin (Terhzaz et al. 2015). Concurring with this prior evidence, our data found a constitutive increase in Cyp4 and Cyp6 family transcript levels in 91-R with the co-occurring Cyp4e3 downregulation analogous to other insecticide-resistant insects (Yang and Liu 2011). Furthermore, DDT selection could induce elevated levels of reactive oxygen species (ROS), resulting in mitochondria-mediated apoptotic changes (Jin et al. 2014). Since P450s respond to various endogenous and exogenous compounds (Riddick 2004; Davies et al. 2006), changes in expression may enhance general stress responses, if not also remediation of specific DDT secondary metabolites.

The second most downregulated transcript, Ser12, physically interacts with the thioredoxin peroxidase 2, (Jafrac 2; Guruharsha et al. 2011). Thioredoxin peroxidases (Tpx) catalyze antioxidant reactions to neutralize ROS but also function in the regulation of apoptosis. The thioredoxin (Trx) redox system is based on the reversible oxidation of Trx cysteine disulfide bonds (S2) to cysteine-thiol (SH2) via concomitant reduction of equimolar equivalents of S2, which is catalyzed by the NADPH-dependent thioredoxin reductase (TrxR). Trx remediates accumulated ROS through the reduction of superoxide dismutase (sod) generated H2O2, often through cytosolic thioredoxin peroxidase (Tpx) or other oxidized (damaged) proteins. The Trx pathway is activated by exposure to hypoxic conditions or ROS to remediate conditions that lead to cell damage via oxidation of lipids, nucleotides, and proteins and to produce signals that increase the rate of cell growth and elevate cellular resistance to apoptosis. Furthermore, cell damage and subsequent generation of ROS act as second messengers within signal transduction cascades (Thannickal and Fanburg 2000), including the induction of apoptosis (Sarafian and Bredesen 1994). The endoplasmic-reticulum localized Jafrac2 is redistributed to the cytoplasm during apoptosis, where it interacts with the Drosophila inhibitor of apoptosis proteins (IAPs), DIAP1 and 2 (thread; Tenev et al. 2002). These DIAP family members possess a Baculovirus IAP Repeat (BIR) domain that acts to bind and sequester proteolytically active caspases that would otherwise initiate apoptosis in an unbound state (Crook et al. 1993). The RING finder domain of DIAP1, moreover, regulates the ubiquination, and possible endosomal degradation, of DRONC. Therefore, DIAP1 binding and inactivation of the three Drosophila caspases drICE (Kaiser et al. 1998), DCP-1 (Hawkins et al. 1999), and DRONC (Meier et al. 2000; Quinn et al. 2000) suppresses apoptosis. This release of Jafrac2, however, is contingent upon N-terminal cleavage exposing an IAP-binding motif (IBM) to allow subsequent competitive displacement of the active form of the caspase DRONC from DIAP1. Based upon this observed colocalization, we propose that the serine protease Ser12 may constitute an endoplasmic reticulum protease able to N-terminal cleave Jafrac2, which upon subsequent downregulation in 91-R may hinder transport of Jafrac2 to the cytoplasm, thereby suppressing the onset or progression of apoptotic events.

Serine proteases are a diverse family of proteolytic enzymes involved in a wide array of extracellular processes, including digestion, immune response, degradation of extracellular matrix components, and inflammation, along with intracellular protein cycling, apoptosis, and mediation of signal transduction events (Weiss 1989; Puente et al. 2005; Ramachandran and Hollenberg 2009). Four serine-like proteases, Jon99Fii, CG4927, CG17234, CG30025, and Ser12 are downregulated in 91-R. Previous research showed that Ser12 mutation can alter epithelial wound healing capacity (Campos et al. 2010) and is analogous to GO annotations and RNAi knockdown for PPO 1 and the glycoside hydrolase/chitinase II and imaginal disc growth factor 1, Idfg1 (Pesch et al. 2016). Systemic wound response is orchestrated by the serine protease-mediated conversion of PPO to phenoloxidase (PO) (Bidla et al. 2009), which in turn converts monophenols to o-quinones during melanin nodule and scab formation (Galko and Krasnow 2004), with the hemolymph hayan protease (CG6361) cleaving PPO to the activated PO (Nam et al. 2012). Additionally, c-JNK is activated in response to cell wounding (Rämet et al. 2002) and triggers cell migration during wound closure (Lesch et al. 2010).

Since PO-catalyzed melanization reactions result in the production of ROS (Nappi and Vass 1993), an integrated system of cellular response has been described. Specifically, serine proteases were shown to trigger JNK-dependent activation of redox pathways in D. melanogaster (Nam et al. 2012), which may serve to protect cells from oxidative damage during wound healing. Results from this study disclose a constitutive downregulation of PPO1 and 2, as well as an upregulation of serine protease inhibitors, serpin Spn43Aa and the Kunitz trypsin-inhibitor domain-containing, sfp24Ba, in 91-R. These suggest that adaptations resulting in decreased endogenous ROS production in response to (xenobiotic-induce) cell damage may counterbalance the excess produced when detoxifying high levels of xenobiotics like DDT. Since hemolymph ROS are systemic, organism-wide upregulation of these protective redox mechanisms may be selectively advantageous.

Three serine-like proteases, CG7542, Jon99Ci, and Jon65Aii, are upregulated in 91-R, with the latter two involved in Toll-signaling. The Drosophila Toll pathway helps to modulate cellular immune response as well as several developmental events (Valanne et al. 2011) and is regulated via serpin family protease inhibitors that also directly control melanization via phenol oxidases (Ligoxygakis et al. 2002). Products from the upregulated CG43055 and Lectin-galC1 (located in selective sweep 5) transcripts have C-lectin domains that define extracellular, carbohydrate-binding, Ca2+-dependent cell surface adhesion proteins involved in agglutination and hemocyte-mediation mechanism during bacterial infection response (Tanji et al. 2006). C-type lectins initiate immune response via a diverse set of signaling pathways including through PPO1 and PPO2 (Geijtenbeek and Gringhuis 2009).

Conclusions

The regulation of gene expression is an important mechanism by which organisms respond and adapt to their environment (Carroll 2000). In this study, we identified a total of 69 constitutive DEGs in 91-R as compared with 91-C strains, and wherein identified DEGs associated with neural cell development, signal transduction, and gene regulatory cascades. When comparing the current functional data to prior structure mutation data, only 6 out of 69 DEGs were physically linked to 1 of the 13 different genome regions showing evidence of strong directional selection between strains (selective sweeps; Steele et al. 2015). Therefore, it remains difficult to make concrete associations between structural and functional changes in 91-R in relation to the evolution of DDT resistance. For example, even though mutations within 91-R were predicted in the CRE-BP1 and flanking the BEAF-32 consensus TF binding sequence upstream of the most highly upregulated transcript, Unc-115 b, the fact that mutations remaining segregating remains difficult to resolve for subsequent associations with DDT resistance when dominance of the alleles or other factors influencing expression in the genome architecture of 91-R are unknown. Furthermore, in spite of evidence that selective sweeps have likely occurred in response to directional selection of DDT exposure in 91-R, the specious associations with ancestral neutral mutations embedded within these sweeps, as well as across the entire genome, that have been impacted by random genetic drift since strains with a common genetic background ∼60 years ago cannot be ruled out. Regardless, the majority of environment-dependent changes in gene expression are predicted to be caused by trans- and cis-acting regulatory mutations (Grundberg et al. 2011; Fear et al. 2016), suggesting that regulatory complexity for gene expression may be important to modulating adaptive divergence. Further studies are undoubtedly required in order to directly correlate structural mutations with complex changes in systemic regulation of gene expression.

Given the above caveats, this current study does provide novel evidence that alteration of neurological function or neurogenesis may be associated with DDT resistance in 91-R. Specifically, this is shown for the most highly upregulated transcript in 91-R, the actin-binding Unc-115 b, that functions in axon guidance and maintenance, and the most downregulated transcript, the angiopoietin-like CG31832, which has a putative vascular and neural guidance function. One might suggest that the most dramatic changes in constitutive gene expression, both impacting the regulation of neural cell migration or maintenance (survival/resistance to apoptosis), may not be coincidental, and possibly be contributing factors in the adaption of 91-R to high DDT exposure levels. Furthermore, the second most downregulated transcript, Ser12, comprises a proteolytic putative activator of the thioredoxin 2, Jafrac2. This reduced expression of Ser12 might be hypothesized to reduce to proteolytic cleavage of Jafrac2 and subsequently inhibit cell progression into apoptosis. While the current study provides enticing insights into the evolution of DDT resistance in 91-R, additional functional evidence is required to establish causation with respect to axon growth and cellular resistance to apoptosis pathways. Additional experiments utilizing analysis of reselected recombinant inbred lines (RILs), and reverse genetics will likely be required to identify the causal basis of DDT resistance in 91-R, thus results here must be interpreted with caution.

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.