Global Transcriptional Profiling of Diapause and Climatic Adaptation in Drosophila melanogaster.

Wild populations of the model organism Drosophila melanogaster experience highly heterogeneous environments over broad geographical ranges as well as over seasonal and annual timescales. Diapause is a primary adaptation to environmental heterogeneity, and in D. melanogaster the propensity to enter diapause varies predictably with latitude and season. Here we performed global transcriptomic profiling of naturally occurring variation in diapause expression elicited by short day photoperiod and moderately low temperature in two tissue types associated with neuroendocrine and endocrine signaling, heads, and ovaries. We show that diapause in D. melanogaster is an actively regulated phenotype at the transcriptional level, suggesting that diapause is not a simple physiological or reproductive quiescence. Differentially expressed genes and pathways are highly distinct in heads and ovaries, demonstrating that the diapause response is not uniform throughout the soma and suggesting that it may be comprised of functional modules associated with specific tissues. Genes downregulated in heads of diapausing flies are significantly enriched for clinally varying single nucleotide polymorphism (SNPs) and seasonally oscillating SNPs, consistent with the hypothesis that diapause is a driving phenotype of climatic adaptation. We also show that chromosome location-based coregulation of gene expression is present in the transcriptional regulation of diapause. Taken together, these results demonstrate that diapause is a complex phenotype actively regulated in multiple tissues, and support the hypothesis that natural variation in diapause propensity underlies adaptation to spatially and temporally varying selective pressures.

Introduction

Understanding the mechanisms by which populations adapt to distinct environments has remained a major goal of evolutionary and ecological genetics. A number of environmental parameters vary predictably with latitude, and the examination of natural populations across latitudinal gradients has been widely used to elucidate fundamental aspects of spatially variable selection (Blanckenhorn and Fairbairn 1995; Clapham et al. 1998; Huey et al. 2000; Loeschcke et al. 2000; García Gil et al. 2003; Stinchcombe et al. 2004; Zhang et al. 2008). Similarly, seasonality creates predictable temporal variation in a suite of environmental parameters, and changes in natural populations for various fitness-related traits among seasons demonstrate evolutionary responses to seasonal shifts in selective pressures (Hairston and Dillon 1990; Grant PR and Grant BR 2002; Korves et al. 2007; Brown et al. 2013; Tarwater and Beissinger 2013).

The model organism Drosophila melanogaster originated in tropical areas of Africa and has subsequently spread into temperate regions on multiple continents (David and Capy 1988). Natural populations now experience highly heterogeneous environments over broad geographical ranges, thus providing an excellent system to study the evolutionary responses to environmental variability. A number of latitudinal clines have been described in D. melanogaster, both at the phenotypic (Capy et al. 1993; James et al. 1997; Azevedo et al. 1998; Karan et al. 1998; Robinson et al. 2000; Schmidt, Matzkin, et al. 2005) and genetic levels (Berry and Kreitman 1993; Verrelli and Eanes 2001; Frydenberg et al. 2003; Sezgin et al. 2004; McKechnie et al. 2010); clinal patterns have been documented on most continents, including Australia (Mitrovski and Hoffmann 2001; Hoffmann et al. 2002; Frydenberg et al. 2003; Umina et al. 2005; Paaby et al. 2010), the Indian subcontinent (Rajpurohit, Parkash, and Ramniwas 2008; Rajpurohit et al. 2008; Rajpurohit and Nedved 2013), North America (Schmidt, Matzkin, et al. 2005; Schmidt, Paaby, et al. 2005; Schmidt et al. 2008; Paaby et al. 2010), South America (Folguera et al. 2008; Goenaga et al. 2013), and Africa (Fabian et al. 2015). The direct link between clinal genetic markers and fitness-related phenotypes (Schmidt et al. 2008; Paaby et al. 2010, 2014; Lee et al. 2013), as well as repeated clinal variation at multiple continents (Turner et al. 2008; Reinhardt et al. 2014) suggest that despite demographic processes (Bergland et al. 2015; Kao et al. 2015), local adaptation to spatially varying selective pressures has contributed to the observed phenotypic and genetic clines in D. melanogaster. Repeated and predictable changes in phenotype (Schmidt and Conde 2006; Behrman et al. 2015) and allele frequency (Cogni et al. 2013, 2014) over seasons have also been found in D. melanogaster. Genomic screen has identified hundreds of single nucleotide polymorphisms (SNPs) that oscillate predictably in frequency over multiple years (Bergland et al. 2014), indicating the widespread operation of balancing selection over seasonal and annual timescales as well as across spatial gradients.

The selective regimes associated with high and low latitudes are in many ways similar to those associated with winter and summer seasonality: Both high latitude and winter seasons are associated with lower temperature and scarce resources, while both low latitude and summer seasons are associated with elevated temperature and abundant resources. A natural question is whether the observed latitudinal clines and seasonal oscillations are generated by the same selection pressures, and if so, what are the common underlying mechanisms of adaptation.

The adult ovarian diapause in D. melanogaster may be one phenotype underlying adaptation to the cyclical seasonal and geographic gradients in climatic parameters and associated selection regimes. Diapause is induced by shortened photoperiod and moderately low temperature, and results in reproductive quiescence, extreme life span extension (more than 5-fold), negligible senescence, increased lipid storage, and increased stress tolerance (Saunders et al. 1989; Saunders and Gilbert 1990; Tatar, Chien, et al. 2001; Schmidt, Matzkin, et al. 2005; Schmidt, Paaby, et al. 2005). Notably, diapause expression is highly variable within and among D. melanogaster populations, and the variation has a significant genetic component: When exposed to the diapause-inducing conditions, some genotypes tend to become reproductively quiescent whereas others tend to proceed with vitellogenesis and reproductive development (Williams and Sokolowski 1993). The propensity to enter reproductive diapause exhibits a strong latitudinal cline in eastern North American populations: Greater than 90% in temperate areas and ∼30% in neotropical habitats (Schmidt, Matzkin, et al. 2005). Diapause incidence also varies with season: With high propensity to enter diapause (∼90%) in spring and lower incidence (40–50%) in fall (Schmidt and Conde 2006). Similarly, genetic polymorphism in candidate diapause genes also varies clinally (Schmidt et al. 2008; Fabian et al. 2012; Bergland et al. 2014) and seasonally (Cogni et al. 2013). Furthermore, diapause incidence is shown to be a powerful predictor of a suite of other traits associated with organismal fitness: Even without exposure to diapause-inducing conditions, genotypes that show a high diapause rate are constitutively longer lived, less fecund, and more stress resistant than the lines that show a low diapause incidence (Schmidt, Paaby, et al. 2005; Schmidt and Paaby 2008). These correlated traits have also been shown to be clinal (Hoffmann et al. 2002; Schmidt, Paaby, et al. 2005) as well as seasonal.

The clinal and seasonal patterns of diapause incidence and its association with other life history traits may represent a fundamental trade-off between somatic maintenance and reproduction: Under stressful and unfavorable conditions, hardiness is favored although associated with lower fecundity; while under benign conditions, high reproductive rate is favored at the cost of aspects of survival. It remains unclear what traits are directly under selection to generate the observed latitudinal clines and seasonal oscillations in phenotype and allele frequencies. Given the physiological consequences of diapause induction in D. melanogaster, the clinal and seasonal variation in diapause incidence, and the genetic correlations between diapause propensity and other fitness-related traits, we hypothesize that diapause is one major determinant of adaptation to spatial and temporal environmental heterogeneity in this taxon. If the natural variation of diapause does play critical roles in adaptation to environmental heterogeneity, genes differentially regulated as a function of the diapause phenotype are likely under spatially and/or temporally varying selectively pressures, thus genetic polymorphisms on these genes or in vicinity of these genes are likely to show clinal and seasonal patterns as they may have distinct cis-regulatory functions that regulates gene expression in diapausing nondiapausing individuals.

Unlike many other arthropod systems in which diapause is constitutive and invariant within populations under stressful environments (Tauber et al. 1986; Denlinger 2002; Ragland et al. 2010; Poelchau et al. 2013a), the intraspecific variation of diapause propensity in D. melanogaster provides an excellent opportunity to identify the genes and pathways associated with the observed genetic variance for diapause expression, and to identify the downstream targets that may ultimately result in the observed phenotypes, via comparisons between diapausing and nondiapausing genotypes exposed to identical environmental parameters. Three candidate genes for diapause in D. melanogaster have been previously identified: Dp110 (Williams et al. 2006), timeless (tim) (Tauber et al. 2007), and couch potato (cpo) (Schmidt et al. 2008). Similarly, two major pathways have been implicated in diapause expression: The circadian clock (Tauber et al. 2007; Pegoraro et al. 2014) and insulin/insulin-like growth factor signaling (IIS) that plays a central role in the neuroendocrine regulation (Tatar, Kopelman, et al. 2001; Williams et al. 2006; Schmidt 2011). However, a comprehensive genomic and transcriptomic analysis of diapause in Drosophila is lacking (but see Baker and Russell 2009, which utilizes diapause as a way to synchronize developmental stages to elucidate fundamental aspects of egg development).

Maternal photoperiod has been shown to influence whether offspring express diapause in other taxa (Henrich and Denlinger 1982; Saunders 1987; Rockey et al. 1989; Oku et al. 2003), and epigenetic processes are known to contribute to the control of diapause (Reynolds et al. 2013; Yocum et al. 2015). Accordingly, we hypothesized that higher level regulatory mechanisms such as local sharing of regulatory elements and chromatin structure may also be involved in the regulation of diapause in D. melanogaster, and seek to test if this is reflected in the differential expression profile of diapausing and nondiapausing individuals.

Here, we performed high-coverage, genome-wide transcriptomic profiling of heads and ovaries, two tissue types associated with neuroendocrine and endocrine signaling, in diapausing and nondiapausing individuals. To our knowledge, this is the first comprehensive expression profile of diapause in an organism that has naturally different perception of the same environmental cues. We use these data to address several fundamental questions: 1) What are the genes and transcripts that are differentially expressed (DE) as a function of the diapause phenotype? 2) Are the genes previously identified as being associated with variance in diapause propensity differentially regulated in diapausing and nondiapausing genotypes? 3) Is there evidence of chromatin-based gene regulation associated with the diapause phenotype? 4) Do the genes that are DE as a function of diapause also segregate for SNPs that vary predictably in frequency with latitude and season? And 5) Are the observed patterns parallel between heads and ovaries or are they tissue specific? Our data demonstrate that diapause is an actively regulated phenotype at the transcriptional level, suggesting that it is not a simple dormancy characterized by physiological quiescence. All genes previously shown to affect diapause have DE isoform(s) in at least one of the two tissues, confirming their association with the phenotype. Patterns of differential expression are spatially clustered across the major chromosomes, suggesting that chromatin-level regulation may contribute to the regulation and response of diapause. Genes downregulated in heads of diapausing flies are significantly enriched for clinal and seasonal SNPs, consistent with the hypothesis that diapause is a driving phenotype of climatic adaptation. Taken together, these results support the idea that diapause is a fundamental and complex phenotype that involves multiple levels of regulation, and underlies adaptation to spatially and temporally varying selective pressures.

Results

Differentially Expressed Genes

Our data demonstrate that diapause is actively regulated at the transcriptional level in both heads and ovaries in adult D. melanogaster (fig. 1). The numbers of genes having relatively higher and lower expression levels in diapausing individuals are approximately the same in both head and ovary, indicating that diapause is not a simple dampening of biological processes, but represents an alternative physiological state associated with the concomitant upregulation and downregulation of many genes.

F

Log2-fold change of expression levels versus mean normalized RPKM values in head and ovary. Genes that are significantly differentially expressed (FDR < 0.05) in D and ND are shown in red. The expression levels are normalized by library sizes. Genes with relatively higher expression levels in diapausing samples as compared with nondiapausing samples have positive log2-fold change values, and genes with lower expression levels in D have negative log2-fold change values. In the head, 6,237 genes fall above the line of y = 0 and 5,967 genes fall below the line. In the ovary, 5,272 genes fall above the line and 5,584 genes fall below.

A large number of genes are DE between diapausing (D) and nondiapausing (ND) individuals and this is true both in head and ovary (∼9% in all tested genes in head and ∼11% in all tested genes in ovary; table 1). Because D. melanogaster has only been in temperate environments for a relatively short period of time (David and Capy 1988), one hypothesis is that as a potentially recent adaptation to the novel temperate environment, diapause in this species might only involve a small number of genes. However, our results suggest that as in other insects (Ragland et al. 2010; Poelchau et al. 2013a, 2013b), the diapause phenotype in D. melanogaster is associated with differential expression of hundreds of genes.

Table 1.

Summary of the Gene-Level Tests of Differential Expression.

	Head	Ovary
Candidate genes tested for DEa	12,204	10,856
Differentially expressed genes between D and ND (FDR < 0.05)	1,094 (8.96%)	1,173 (10.81%)
Upregulated genes in D (FDR < 0.05)	510	426
Downregulated genes in D (FDR < 0.05)	584	747

Note.—Number in parenthesis is the percentage of DE genes in all tested candidate genes.

aTo eliminate noise and improve model fitting, only genes with RPKM > 10 in either D or ND sample are tested for differential expression.

The expression profiles of D and ND are highly distinct between head and ovary, and DE genes in the two tissues are not the same (supplementary table 1, Supplementary Material online). The reads per kilobase per million reads (RPKM) fold change of D over ND in the head is not correlated to that in the ovary (P value of Pearson’s product-moment correlation test = 0.9817, correlation = −0.000222). Only 165 genes are DE in both head and ovary, of which 41 show opposite regulation in the two tissues; by chance alone, we would expect to see ∼120 genes overlapping in the DE lists in the two tissues. The inconsistency of transcriptional regulation in head and ovary as a function of the diapause phenotype advocates analysis of different tissue and body parts in the dissection of mechanisms of diapause.

Because the most obvious phenotypic difference between D and ND is the development and maturation of the ovary, and genes involved in vitellogenesis and oogenesis are at much higher expression levels in developed ovaries compared with development-arrested ovaries (Baker and Russell 2009), it is possible that the drastic difference of these genes in D and ND ovaries alone can affect the detection of other DE genes. However, after removing reads mapped to these genes and rerunnning the differential expression analysis, the estimation of expression levels of the other genes and the differential expression calling was not affected. We therefore only report the results with all genes included in the analysis.

Because diapause incidence is clinal (Schmidt, Matzkin, et al. 2005), and distributions of some of the cosmopolitan inversions are also clinal (Knibb 1982; Kapun et al. 2014), one hypothesis is that the expression of diapause is associated with specific inversion status. We assessed the frequencies of the seven major cosmopolitan inversions (In(2L)t, In(2R)Ns, In(3L)P, In(3R)C, In(3R)K, In(3R)Mo, In(3R)Payne) and found that In(2R)Ns, In(3R)C, In(3R)K, and In(3R)Mo are at low frequencies in our experimental population. The frequencies of In(2L)t, In(2R)Ns, In(3R)C, In(3R)K, and In(3R)Mo are not different between diapausing and nondiapausing samples (supplementary table S1, Supplementary Material online). The regions containing the inversion-specific SNPs of In(3L)P and In(3R)Payne are not transcribed so we are not able to assess their frequencies from our data. However, according to previous results (Kapun et al. 2014), they are at low frequencies in the northern populations of North America. In addition, the DE genes are found across all major chromosomes and are not clustered around inversion break points. Based on these results, we conclude that inversion status is not driving the variation of diapause propensity, and does not have pronounced effects on diapause-associated patterns of gene expression in northern populations of North America.

Differentially Expressed Isoforms

To fully characterize the differences in transcriptional profiles of diapausing and nondiapausing Drosophila, we examined patterns of transcription at the gene level incorporating all isoforms. As a complementary analysis, we also examined differential expression of each individual isoform. This is especially interesting in the expression profiles of diapause, as candidate diapause genes may affect the phenotype in transcript-specific manners (Tauber et al. 2007; Schmidt et al. 2008).

Our results show that there are a large number of isoforms DE between D and ND in both the head and ovary. Compared with the gene-level DE detection, the isoform-level tests identified many more genes with at least 1 DE isoform (table 2). As expected, the majority of genes (more than 60%) identified from the gene-level tests have at least 1 DE isoform (fig. 2). However, there are also some genes that do not have any significantly DE isoforms, but when all isoforms are considered together, the genes are DE. There are also hundreds of genes carrying isoforms that have opposite patterns of regulation in D and ND (table 2), and these antagonistic effects may cancel out and make the gene not DE in the gene-level tests. Our results clearly demonstrate that the examination of gene-level DE associated with diapause fails to capture fundamental aspects of the transcriptional complexity.

Table 2.

Summary of the Isoform-Level Tests of Differential Expression.

	Head	Ovary
Candidate isoforms tested for DEa	21,087	18,956
DE isoforms between D and ND (FDR < 0.05)	3,762 (2,608)a	3,603 (2,491)
Upregulated isoforms in D (FDR < 0.05)	1,598 (1,384)	1,784 (1,407)
Downregulated isoforms in D (FDR < 0.05)	2,164 (1,705)	1,819 (1,522)
Genes not DE at gene level, but have DE isoforms	1,939	1,641
Genes with both up- and downregulated isoforms	481	438

Note.—Number in parenthesis is the number of genes corresponding to the isoforms listed before.

aTo eliminate noise and improve model fitting, only isoforms with RPKM > 10 in either D or ND sample are tested for differential expression.

F

Venn diagrams of numbers of genes identified from the gene-level and isoform-level differential expression tests.

Functional Categories of Differentially Expressed Genes

To place the comparisons in gene expression into a biologically meaningful context, we identified enriched functional categories using gene ontology (GO) enrichment analysis as well as the Kyoto encyclopedia of genes and genomes (KEGG) pathways analysis. To distinguish between functional categories of genes that are up- or downregulated as a function of the diapause phenotype, we performed the enrichment tests separately for each class of genes in both heads and ovaries. We see nonoverlapping GO terms (supplementary table S2, Supplementary Material online) and KEGG pathways (table 3) in different tissue types for up- and downregulated genes.

Table 3.

Enriched KEGG Pathways.

Tissue	Pathway	FDR	No. DE in Category	Total Category
Head	Upregulated in D
	Ribosome biogenesis in eukaryotes	5.37 × 10−27	32	71
	RNA polymerase	2.27 × 10−7	9	24
	Ribosome	4.65 × 10−5	12	81
	Pyrimidine metabolism	5.01 × 10−5	12	73
	Downregulated in D
	Metabolic pathways	8.98 × 10−9	74	876
	Glycine, serine, and threonine metabolism	2.47 × 10−7	10	27
	Lysosome	4.69 × 10−5	13	78
	Other glycan degradation	7.32 × 10−5	6	19
	Fatty acid degradation	3.82 ×10−4	7	28
	One carbon pool by folate	1.21 × 10−3	4	11
	Starch and sucrose metabolism	1.39 × 10−3	9	54
	Phenylalanine metabolism	1.51 × 10−3	4	10
	Galactose metabolism	2.03 × 10−3	6	26
	Gycerolipid metabolism	3.61 × 10−3	7	41
Ovary	Upregulated in D
	Phototransduction	3.34 × 10−7	7	25
	Downregulated in D
	No enrichment

Several metabolic pathways are downregulated in diapause heads according to KEGG enrichment, including carbohydrate (starch, sucrose, galactose) as well as amino acid (glycine, serine threonine, phenylalanine) metabolism. This indicates that flies in diapause may generate and/or consume less energy. There is also a downregulation of “fatty acid degradation” and “other glycan degradation” in diapause heads, suggesting a shift toward storage in diapausing flies. This is consistent with physiological characterizations (Tauber et al. 1986) as well as transcriptional profiling of diapause in other taxa (Ragland et al. 2010; Poelchau et al. 2013a).

Candidate Genes and Pathways

Although none of the previously identified candidate genes (cpo, tim, and Dp110) show differential expression at the gene level (table 4), all of them have at least 1 DE isoforms in at least one tissue (table 5). Two isoforms of cpo show opposite regulation in the ovaries: The cpo-RO transcript is upregulated in diapausing ovaries, whereas the cpo-RS transcript is downregulated in this same tissue. The two isoforms of tim, tim-RL, and tim-RM are differentially regulated in opposite directions in both heads and ovaries, and the regulation of each isoform shows antagonistic directions in different parts of the body. This suggests that the effects of tim on diapause are transcript and tissue specific. Taken together, genome-wide transcriptional analysis of naturally occurring variation in diapause induction demonstrates that the three previously known candidate genes underlying diapause propensity are DE; furthermore, the data are of sufficient resolution to detect that they are specific isoforms that are associated with the observed phenotypic variation.

Table 4.

Gene-Level Expression of Candidate Diapause Genes (expression levels are normalized and in units of RPKM).

	Head				Ovary
Gene	D	ND	Adjusted P value	DE	D	ND	Adjusted P value	DE
cpo	46,708	41,759	0.996	False	2,674	1,567	0.377	False
tim	87,867	8,842	0.997	False	5,731	3,611	0.714	False
Dp110	1,987	2,119	0.991	False	4,820	3,855	0.974	False

Table 5.

Transcript-Level Expression of Candidate Diapause Genes (expression levels are normalized and in units of RPKM).

		Head				Ovary
Gene	Transcript	D	ND	Adjusted P value	DE	D	ND	Adjusted P value	DE
cpo	cpo-RN	167	238	0.498	False	367	202	0.785	False
	cpo-RO	5,559	3,835	0.594	False	370	0	0.000	True
	cpo-RP	6,125	5,429	0.992	False	596	332	0.815	False
	cpo-RQ	378	295	0.925	False	18	5	0.566	False
	cpo-RR	Low	Low	NA	NA	Low	Low	NA	NA
	cpo-RS	2,091	3,038	0.081	False	0	198	0.000	True
	cpo-RT	12,130	10,747	0.995	False	112	67	0.865	False
	cpo-RU	Low	Low	NA	NA	Low	Low	NA	NA
	cpo-RV	21,426	17,070	0.987	False	1,382	669	0.579	False
	cpo-RW	96	0	0.000	True	Low	Low	NA	NA
	cpo-RX	Low	Low	NA	NA	Low	Low	NA	NA
	cpo-RY	Low	Low	NA	NA	Low	Low	NA	NA
tim	tim-RB	Low	Low	NA	NA	Low	Low	NA	NA
	tim-RL	22	0	0.000	True	0	35	0.000	True
	tim-RM	351	548	0.027	True	149	38	0.002	True
	tim-RN	0	32	0.000	True	31	7	0.162	False
	tim-RO	625	483	0.912	False	1,224	827	0.951	False
	tim-RP	Low	Low	NA	NA	Low	Low	NA	NA
	tim-RR	8,026	7,547	0.995	False	4,693	2,487	0.841	False
	tim-RS	Low	Low	NA	NA	Low	Low	NA	NA
Dp110	Pi3K92E-RA	1,193	1,539	0.693	False	3,835	3,623	0.981	False
	Pi3K92E-RB	848	525	0.058	False	1,293	0	0.000	True
	Pi3K92E-RC	Low	Low	NA	NA	Low	Low	NA	NA

Note.—NA, not applicable.

The IIS pathway regulates dauer formation in Caenorhabditis elegans (Kimura et al. 1997; Apfeld and Kenyon 1998) and dauer formation is in many ways similar with diapause in insects (Tatar and Yin 2001; Flatt et al. 2013; Sim and Denlinger 2013). Manipulations of different components of the IIS pathway in D. melanogaster phenocopy many aspects of diapause such as arrested cell cycle (LaFever and Drummond-Barbosa 2005; Hsu et al. 2008), extended lifespan (Clancy et al. 2001; Tatar, Kopelman, et al. 2001; Giannakou et al. 2004; Hwangbo et al. 2004; Lee et al. 2008; Grönke et al. 2010), and increased fat storage (Böhni et al. 1999). Therefore, the expression of genes in the insulin-signaling pathway in diapausing and nondiapausing flies is of particular interest.

We examined all genes belonging to the GO category “insulin receptor signaling pathway” (GO: 0008286) to determine if any differential regulation is present in either heads or ovaries as a function of diapause phenotype. None of these genes are DE in the heads. However, multiple components of the signaling pathway exhibit differential expression in the ovaries: The gene melt, which is a modulator of the IIS pathway, is upregulated in diapausing flies (supplementary table S4, Supplementary Material online). The insulin-like peptide Ilp-4 is upregulated, and Ilp7 and Ilp8 are downregulated, in ovaries of diapausing flies (supplementary table S4, Supplementary Material online). Because insulin-like peptides are known to control the germ cell cycle (LaFever and Drummond-Barbosa 2005), this result suggests that these peptides are related to the ovarian developmental arrest. These DE components of the IIS pathway may represent a cascade of signaling responses to upstream signals, although the exact mechanism of the endocrine control remains unclear. As with the evaluation of the candidate genes for diapause, an examination of specific transcripts demonstrates more complex patterns of expression for genes of this pathway in both heads and ovaries. The genes dock, foxo, InR, Pdk1, Pten, and RpL8 have DE transcripts in the heads; the genes chico, dock, melt, Pdk1, Pi3K21B, and Pi3K92E (Dp110) have DE transcripts in the ovaries (supplementary table S4, Supplementary Material online). Of these genes, chico affects body size and lipid storage (Böhni et al. 1999); the expression of chico and foxo is related to lifespan (Clancy et al. 2001; Giannakou et al. 2004; Hwangbo et al. 2004); natural variants of InR are associated with various aspects of life history evolution (Paaby et al. 2014); and Pi3K92E is a candidate diapause gene. Our data suggest that the IIS genes identified as DE in our transcriptional profile are excellent candidates for the regulation of diapause, and that their effects may be transcript specific.

Chromosome Location-Dependent Expression

To test if there is any higher order regulation of expression in diapause, we examined the spatial autocorrelation of log2-fold changes of expression levels in D versus ND along all major chromosomes, to test whether the fold changes tend to be clustered in space based on chromosome locations. The values of Moran’s I, which is a measure of the direction and level of spatial autocorrelation, are located on the far right side of the null distributions generated by randomizing the locations of genes within the chromosome; the sole exception is head data on chromosome X (fig. 3). The pattern of spatial autocorrelation is not driven by local coexpression of tandemly arrayed duplicated genes, as they were excluded from the analysis. This result suggests that chromosome location-based coregulation of gene expression is present in the transcriptional regulation of diapause.

F

The null distributions (shown as curves) and real (shown as vertical lines) Moran’s I of log2-fold change of D over ND on the major chromosome arms in head and ovary. A positive Moran’s I indicates positive spatial autocorrelation. The null distribution was generated by scrambling the starting sites of genes within each chromosome 1,000 times.

Clinal and Seasonal Overlap

To test the hypothesis that genes differentially regulated in diapausing and nondiapausing Drosophila are likely to be under spatially and/or temporally varying selective pressure, we identified the overlaps of DE genes with clinal and seasonal genes according to Bergland et al. (2014). We found hundreds of DE genes harbor clinally varying polymorphisms, and dozens of DE genes harbor seasonally varying polymorphisms (summarized in table 6, listed in supplementary table S5, Supplementary Material online). The enrichment of seasonal/clinal genes in DE gene sets was tested with a block-bootstrapping method as described in Materials and Methods section.

Table 6.

Summary of the Overlap of Differentially Expressed Genes and Clinal/Seasonal Genes

	Head	Ovary
Clinal genes tested for differential expressiona	7,468	6,565
Clinal genes upregulated in diapause	289	279
Clinal genes downregulated in diapause	357	395
Seasonal genes tested for differential expressiona	1,689	1,469
Seasonal genes upregulated in diapause	54	87
Seasonal genes downregulated in diapause	98	75

aTo eliminate noise and improve model fitting, only genes with RPKM > 10 in either D or ND sample are tested for differential expression.

There is a significant enrichment of clinal genes in genes that are downregulated in diapausing heads (P = 0.022) (fig. 4). The 357 clinal genes that are downregulated in diapausing heads are primarily located throughout chromosomes 2 and 3, and the locations of these genes are not associated with the cosmopolitan inversion breakpoints. These genes are enriched for GO terms “proteolysis,” “oxidation–reduction process,” and KEGG “metabolic pathways.” The other classes of DE genes are not enriched for clinal polymorphisms.

F

Clinal and seasonal enrichment of differentially expressed genes. Enrichment (log2 odds ratio) of differentially expressed genes among genes that harbor, or are near clinally varying (left) and seasonally oscillating (right) polymorphisms. Error bars represent 95% confidence intervals based on 500 block bootstrap resampling.

Similarly, we also observe a significant enrichment of seasonal genes in genes that are downregulated in diapausing heads (P = 0.01), but not in other classes of DE genes (fig. 4). As with the clinal overlap gene list, the 98 seasonal genes that are downregulated in diapausing heads are primarily located on chromosomes 2 and 3, and the locations of these genes are not associated with the cosmopolitan inversions. Unlike the clinal overlap gene list, however, no GO or pathway enrichment was identified for this set of genes. Many of these genes’ functions are completely unknown either from experimental evidence or from prediction based on sequence similarities. Interestingly, 29 of the 98 downregulated genes have seasonal SNP(s) in vicinity of the known cis-eQTL peaks of these genes (supplementary table S5, Supplementary Material online) (cis-eQTL regions contributing to regulatory variation in D. melanogaster female head were obtained from King et al. 2014). This suggests that these seasonally fluctuating SNPs, or sequence variants linked to them, may be under seasonally oscillating selection because they are associated with different gene expression levels that are associated with the diapause and nondiapause phenotypes. These results prompt future research on function of these candidate genes, as well as the interplay among segregating variation, expression level, phenotype, and dynamics of selection.

Discussion

Diapause is the most intensively studied adaptation to seasonality in arthropods, and is a complex physiological syndrome that in adults involves sequestered reproduction, extended lifespan, altered metabolism, and increased stress tolerance. Global transcriptional profiles of diapause have been performed on whole bodies of the flesh fly Sarcophaga crassipalpis (Ragland et al. 2010) and the Asian tiger mosquito Aedes albopictus (Poelchau et al. 2013a, 2013b). These studies identified metabolic pathway transitions, cell cycle arrest, and altered endocrine regulation as fundamental features of the diapause program. Both species have constitutive diapause, therefore the difference between the diapause and nondiapause phenotypic states were inferred from comparisons between insects exposed to different environments as well as insects at different diapause stages. In D. melanogaster, the expression of glucagon- and insulin-like peptide genes was shown to increase during diapause, and the expression of several target genes of these peptides is altered by the diapause phenotype (Kubrak et al. 2014). Transcriptomic profiling of female abdomens has identified lipid metabolism, insulin signaling, and structural constituents of chitin-based cuticle being altered by the transition from diapause-inducing conditions to diapause-terminating conditions (Baker and Russell 2009).

Here we provide the first global transcriptomic profiling of naturally occurring variation in diapause expression that is elicited under the exact same environmental parameters. Like the aforementioned studies, we find that diapause is an actively regulated transcriptional state where almost as many genes are upregulated as downregulated. We have also identified that aspects of metabolism are altered as a function of diapause. Components of the IIS pathway were identified as DE in our study, consistent with previous investigations in dauer formation in C. elegans (Tatar and Yin 2001; Sim and Denlinger 2013) and diapause in mosquitos (Sim and Denlinger 2008; Sim et al. 2015). Unlike the aforementioned studies, our analysis focused on how differential effective perception of low temperature and a short photoperiod affected patterns of transcription in the adult head and ovary, tissues/structures associated with juvenile hormone and ecdysteroid signaling rather than systemic metabolism (e.g., fat body, flight muscle, digestive system). Although both the adult head and ovary may be involved in both environmental perception and the downstream transcriptional regulation of this perception (Liu et al. 1988; Plautz et al. 1997), diapause expression in Drosophila is regulated by measuring the length of night (Meuti and Denlinger 2013; Saunders 2013) and may involve an inhibition of juvenile hormone release from the corpus allatum (Saunders and Gilbert 1990; Yamamoto et al. 2013). Further investigations have also implicated ecdysteroids in the regulation of diapause in this taxon (Richard et al. 1998). Thus, transcriptional profiles of flies differentially responding to the same environmental parameters eliminate confounding factors such as diet, age, developmental state, temperature, and photoperiod. Confining the transcriptional profile to heads and ovaries reduces the influence of systemic metabolic processes in other parts of the body. This type of global transcriptional screen may reveal genes and pathways differentially responding to specific environmental cues, as well as those under distinct neuroendocrine and endocrine control. The fat body is responsible for a suite of metabolic and biosynthetic activities, therefore expression patterns in the fat body in diapausing and nondiapausing individuals may reflect distinct downstream metabolic changes as a consequence of the differential environmental perception. However this was not tracked in this study.

Biological Processes Revealed by GO and KEGG Enrichment Tests

Among genes that are upregulated in diapausing head compared with nondiapausing head, GO terms such as “ribosome biogenesis,” “ribonucleoprotein complex biogenesis,” “rRNA (ribosome RNA) processing,” and “rRNA metabolic process” are among the most significantly enriched biological processes. Ribosome biogenesis is a major metabolic activity and accounts for large proportion of energy consumption in the cell (Wullschleger et al. 2006). Therefore, at first glance it is paradoxical that genes related to ribosome biogenesis are upregulated in diapausing heads, where energy consumption is presumably reduced. However, in the flesh fly S. crassipalpis, the ribosomal protein P0 was found to be upregulated in diapausing pupae brains (Flannagan et al. 1998), and the upregulation of P0 is tightly linked to the downregulation of O2 consumption (Craig and Denlinger 2000). Ribosome biogenesis proteins were thought to contribute to the cellular maintenance of the insect during the diapause period (Denlinger 2002). Although the homologous gene in D. melanogaster, RpLP0, is not significantly upregulated in diapause heads, the upregulation of these functional categories suggests that low metabolic rate and similar response to hypoxia occur in diapausing heads in D. melanogaster.

Vitellogenesis and egg formation related terms are enriched in genes that are downregulated in diapausing head. Although vitellogenesis and egg formation obviously does not occur in the head, these genes are likely pleiotropic and may have important functions outside the ovary. In fact, the expression levels of the yolk proteins Yp1, Yp2, and Yp3 are dozens of times higher in heads than in ovaries, and experimental evidence has shown that these genes are involved in neurogenesis (Neumüller et al. 2011). The consistency of the shut down of some vitellogensis and egg formation related genes in diapausing heads and ovaries indicates some shared mechanisms of gene regulation in different parts of the body.

The phototransduction pathway is upregulated in diapausing ovaries. Phototransduction genes are expected to play a role in the regulation of diapause as they are involved in light sensing, and the perceived change in photoperiod is an anticipatory environmental cue used in both the initiation and termination of diapause (Denlinger 2002; Pegoraro et al. 2014). The conundrum is why the upregulation of phototransduction occurs in the ovary instead of head. It was previously shown that circadian genes are expressed in multiple tissues (Liu et al. 1988) and circadian oscillators are present throughout the body of D. melanogaster (Plautz et al. 1997). It may be that the ovaries sense light directly, and that increased light sensing in diapausing ovaries affects the expression of diapause.

Cell Cycle Arrest in the Ovaries

GO terms such as “chromosome organization,” “nucleosome assembly,” “protein–DNA complex assembly,” and “chromatin assembly” are highly over represented in genes that are downregulated in diapausing ovaries. A closer examination reveals that these GO enrichments are driven by the shut down of a group of ∼90 histones located on the right end of chromosome 2L from cytogenetic location 39D3 to 39E1. Canonical histones are greatly upregulated during S phase of the cell cycle for the assembly of newly replicated chromatin (Osley 1991), but the two genes H3.3A and H3.3B that encode the variant H3.3 are expressed at a relatively constant level throughout the cell cycle (Ahmad and Henikoff 2002). In the ovary data, there is a bulk loss of canonical histone transcripts in diapausing individuals but H3.3 is not DE between diapausing and nondiapausing individuals: This suggests arrested cell division (somatic and germline) during diapause. In fact, GO terms such as “cell cycle,” “cell division,” “meiotic cell cycle,” and “mitotic cell cycle” are also enriched in genes that are downregulated in diapausing ovaries.

However, many of these histone encoding genes are characterized by a high degree of sequence similarity, and thus it is possible that the mapping software is not able to accurately differentiate between them. If true, assigned expression values may have been generated for similar histones by dividing reads equally among them; thus, the observed differential expression of these functional classes may be driven by a subset of the associated genes. Despite that, the observed data for head and ovary are still quite distinct: The canonical histones are shut down in diapausing ovaries but not in diapausing heads, suggesting that the cell cycle arrest is specific to the ovaries.

Chromosome Location Dependent Expression Regulation

Neighboring genes in the D. melanogaster genome can exhibit similar expression patterns across different experimental conditions, and such similarity is not explained by gene function or homology (Boutanaev et al. 2002; Spellman and Rubin 2002; Kalmykova et al. 2005; Levine et al. 2010). Although the mechanism of this neighborhood effect is not clear, we speculated that it can be explained by regulation at the level of chromatin structure, where chromosome regions containing genes open for translational regulation increase the accessibility of the cis-regulatory elements of their neighboring genes. Alternatively, local sharing of cis- and trans-regulatory elements could activate the target genes as well as its adjacent genes (Oliver et al. 2002; Michalak 2008).

We have also observed a high level of spatial autocorrelation of transcriptional regulation as a function of the diapause phenotype. Whether the coregulation of genes involved is constitutive or specific to the diapause phenotype is beyond the scope of our data and analyses. However, we can make two speculations: 1) Higher order regulation of gene expression, such as regulatory element sharing and chromatin remodeling, is present in diapause; and 2) because the transcriptional regulation is not always perfectly precise to individual genes, it is possible that some genes are DE as a result of “transcriptional hitch-hiking,” where the neighboring genes are the targets of transcriptional regulation (Oliver et al. 2002). Regardless, the spatial autocorrelation of transcriptional profiles does not conflict with our previous conclusion that diapause involves the differential expression of many genes, as DE genes are found across all chromosomes.

Seasonal and Clinal Overlap

Evidence suggests that the ancestral source of D. melanogaster is the subtropical Miombo and Mopane woodland of Zambia and Zimbabwe, where the genetic diversity is the highest (Pool et al. 2012). Because there is wet–dry seasonality in this region, it is possible that D. melanogaster has an ancestral seasonal syndrome or response that is associated with the wet/dry seasonality. Migration out of subtropical Africa to the high latitude temperate regions in North America, Europe, Australia, and the Indian subcontinent with pronounced winters has resulted in the exposure of D. melanogaster to novel environmental conditions, as well as selective pressures, associated with temperate climates. The physiological and genetic mechanism of adaptation to temperate climate is not clear but increased diapause propensity is likely involved (Fabian et al. 2015), and this may in part be achieved by selection on gene expression levels.

Our global transcriptional profiling of diapause has suggested that, as an adaptation to seasonality and temperate climate, diapause in D. melanogaster is a globally regulated syndrome. The genes that are DE as a function of the naturally occurring variation in diapause phenotype represent a nonrandom set, as genes that are downregulated in diapausing heads are more likely than expected by chance to show both seasonal and clinal patterns of adaptation. Contrary to the traditional hypothesis that genes that are shut down during diapause may simply represent a suspension of development or reduced metabolic activity (Denlinger 2002), our data demonstrate that it is the transcriptionally suppressed rather than transcriptionally activated genes that exhibit a robust signal of climatic adaptation. We suggest that subsequent functional analysis of this set of genes could provide further insight into the mechanisms of climatic adaptation in Drosophila and other taxa.

Materials and Methods

Fly Populations

Drosophila melanogaster was collected as isofemale lines by direct aspiration on wind fallen fruit from four natural populations of northern regions of the east coast of the United States: Rocky Ridge Orchard in Bowdoin, ME (44.03°N, 69.95°W), Champlain Orchards in Shoreham, VT (43.86°N, 73.35°W), Westward Orchards in Harvard, MA (42.49°N, 71.56°W), and Lyman Orchards in Middlefield, CT (41.50°N, 72.71°W). All populations were sampled in October 2009. Fifty isofemale lines from each orchard were used to construct an outbred population. Ten females and ten males from a single age cohort for each isofemale line were released into laboratory population cages and were allowed to interbreed for six generations. The outbred cage was constructed less than 6 months after the flies were collected from the wild. Four cages were created from the F1 of the initial cohort, and the population was maintained as four experimental replicates.

Diapause Assay, cDNA Library Construction, and Sequencing

Vials containing standard cornmeal-molasses medium were placed in the population cage and flies were allowed to lay eggs on the surface for approximately 4 h. These eggs were reared at low density at 25 °C, 12 h light: 12 h dark. Females were collected within 2 h posteclosion and placed at 11 °C at a photoperiod of 9 h light: 15 h dark. After being exposed to the diapause-inducing conditions of low temperature and short days for 4 weeks, all experimental females were dissected and the developmental status of the ovaries was assessed (King 1970). A female was scored as D if the most advanced oocyte was arrested before vitellogenesis (before stage 8); a female was scored as ND if vitellogenin was observed in either ovary (stage 8 or later).

The heads and ovaries of either diapausing or nondiapausing females were pooled and separately flash frozen in liquid nitrogen, yielding 4 samples: D_head, ND_head, D_ovary, and ND_ovary. A total of 92 flies went into the pooled diapause samples (D_head and D_ovary), and a total of 106 flies went into the pooled nondiapause samples (ND_head and ND_ovary). cDNA libraries were prepared as previously described (Elliott et al. 2013), and sequenced on the Illumina Hiseq 2000 platform using standard protocols for 100-bp single-end read sequencing.

Read Processing, Alignment, Counts Estimation, and Differential Expression

As some of the fragments in the libraries were shorter than 100 bp, the 3′ adapters were sequenced in part or full in some of the reads. Therefore raw reads were first trimmed using the cutadapt program (Martin 2011). Trimmed reads were mapped to the D. melanogaster reference transcriptome and expression levels of genes and isoforms were estimated following the RNA-Seq by Expectation-Maximization (RSEM) pipeline version 1.2.8 (Li and Dewey 2011). The reference transcriptome was constructed with the D. melanogaster reference genome version 5.54 and the Ensembl annotation version 5.74. The mean and standard deviation of read length in each library was calculated using fastqstats of the ea-utils toolset (Aronesty 2013), and provided to the RSEM pipeline.

The empirical Bayesian approach-based tool EBSeq (version 1.5.3) (Leng et al. 2013) was implemented to call for genes and isoforms that have differential expression levels between D_head and ND_head, D_ovary and ND_ovary. The RPKM values given by RSEM were used as input. Library sizes were normalized using the median normalization approach provided by EBSeq. Genes and isoforms that have adjusted P-values (false discovery rate [FDR]) smaller than 0.05 were considered DE as a function of the diapause phenotype. Due to the similarities of the isoforms, it is difficult to assign short reads unambiguously to individual isoforms (Roberts and Pachter 2011). Therefore we confine our examination of the isoform-level expression profiles to a general analysis, and do not make further set-based inferences.

Because the most obvious phenotypic difference between D and ND is the development and maturation of the ovary, and genes involved in vitellogenesis and oogenesis are at much higher expression levels in developed ovaries compared with development-arrested ovaries (Baker and Russell 2009), it is possible that the drastic difference of these genes in D and ND ovaries alone can affect the detection of other DE genes. To address this issue, we excluded the genes under the GO term “vitellogenesis” (GO:0007296) and “oogenesis” (GO: 0048477) and reran the differential expression analysis on the ovary data. The estimation of the expression levels for the other genes is not affected by the removal of vitellogenesis and oogenesis genes, and the list of DE genes is not affected. Therefore, these genes were kept in the analysis.

To estimate the inversion frequencies in the samples, we called SNPs from the RNAseq data using the SAMtools mpileup version 1.1 (Li et al. 2009), and assessed the frequencies of major cosmopolitan inversions using the diagnostic SNPs provided by Kapun et al. (2014).

Functional Category and Pathway Enrichment

GO and KEGG pathway enrichment analysis was performed using the GOseq package (Young et al. 2010) that accounts for transcript length bias associated with RNA-seq data. Only genes that had RPKM values greater than 1 in at least one of the D or ND samples in the respective tissues were used as the background gene list for that tissue. Genes that were of significantly high and low expression in diapausing head/ovary samples were separately tested to identify which categories of genes are enhanced and suppressed as a function of diapause. The P values were corrected using the Benjamini and Hochberg procedure (Benjamini and Hochberg 1995), and the FDR threshold was set to be 0.05.

Location-Dependent Expression

To investigate whether genes located closely on the chromosome are similarly regulated between D and ND, we examined the log2-fold changes of expression levels in D versus ND along each major chromosome arm, and calculated the levels of spatial autocorrelation of the log2-fold changes using Moran’s I (Moran 1950). Null distributions of spatial autocorrelation for each chromosome arm were generated by randomly scrambling the locations of genes within the chromosome. Tandemly arrayed duplicated genes were removed from the analysis due to the difficulty of unambiguously mapping reads to these genes, as well as the potential for local coexpression. The list of tandemly arrayed duplicated genes was obtained from Quijano et al. (2008). Moran’s I calculation was implemented with Moran.I function of the R package “ape” (Paradis et al. 2004) version 3.1.

Seasonal and Clinal Overlap with Differentially Expressed Genes

Seasonal and clinal genes were identified based on the results of Bergland et al. (2014), where generalized linear models with binomial error structure were used to call SNPs that oscillate with season or change with latitude. SNPs with a seasonal q value of 0.3 or less were considered to be seasonal, and those with a clinal q value of 0.01 or less were considered to be clinal (as per Bergland et al. 2014). Genes with at least one seasonal/clinal SNP on them, or within 5 kb upstream or downstream of an identified SNP were considered seasonal/clinal genes, respectively.

We then tested if DE genes are enriched in the seasonal/clinal set. Because DE genes and seasonal/clinal SNPs are spatially clustered, to eliminate the possibility that the signal of enrichment is driven by a few clusters of DE genes that happen to harbor clusters of seasonal/clinal SNPs, we implemented a block-bootstrap procedure: We generated 500 subsets of DE genes where at most 1 DE gene was sampled from each 100 kb consecutive interval of the genome. For each DE gene, a gene that was not DE, but of similar length and normalized mean RPKM value, was established as a control. Five hundred sets of control genes were generated for each DE list to ensure the reproducibility of the results. The numbers of genes that are seasonal/clinal were numerated in both the DE and control gene sets, and odds ratios of DE genes being seasonal/clinal relative to the control genes are calculated. Estimates of the mean log2 odds ratios and the 95% confidence intervals were taken across the 500 block bootstraps. Up- and downregulated genes in either tissue type were then tested for enrichment separately.

Supplementary Material

Supplementary tables S1–S5 are available at Molecular Biology and Evolution online (http://www.mbe.oxfordjournals.org/).