key: cord-0002377-g5qkgwti authors: Zhao, Xiaqing; Bergland, Alan O.; Behrman, Emily L.; Gregory, Brian D.; Petrov, Dmitri A.; Schmidt, Paul S. title: Global Transcriptional Profiling of Diapause and Climatic Adaptation in Drosophila melanogaster date: 2015-11-13 journal: Molecular Biology and Evolution DOI: 10.1093/molbev/msv263 sha: 173558a2a7b31fd73e4d89593bde8896efd30a93 doc_id: 2377 cord_uid: g5qkgwti 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. 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 and Nedved 2013) , North America (Schmidt, Matzkin, et al. 2005; Schmidt, Paaby, et al. 2005; 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 Paaby et al. 2010 Paaby et al. , 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 (Cogni et al. , 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 , 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 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 fitnessrelated 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) . 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 708 Zhao et al. . doi:10.1093/molbev/msv263 MBE 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. 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. 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 Poelchau et al. , 2013b , the diapause phenotype in D. melanogaster is associated with differential expression of hundreds of genes. 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 diapauseassociated patterns of gene expression in northern populations of North America. 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 transcriptspecific manners (Tauber et al. 2007; . 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. 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. 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, To eliminate noise and improve model fitting, only genes with RPKM 4 10 in either D or ND sample are tested for differential expression. 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 ). 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, genomewide 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. 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 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€ onke et al. 2010) , and increased fat storage (B€ ohni 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€ ohni 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 To test if there is any higher order regulation of expression in diapause, we examined the spatial autocorrelation of log2fold 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. 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. 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 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. Diapause is the most intensively studied adaptation to seasonality in arthropods, and is a complex physiological Diapause and Climatic Adaptation . doi:10.1093/molbev/msv263 MBE 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 (Poelchau et al. , 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 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 * * * * * * * * * * * * * * * * * * To eliminate noise and improve model fitting, only genes with RPKM 4 10 in either D or ND sample are tested for differential expression. 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. 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 O 2 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€ uller 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. 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) . Diapause and Climatic Adaptation . doi:10.1093/molbev/msv263 MBE 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. 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. 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. As some of the fragments in the libraries were shorter than 100 bp, the 3 0 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 Pvalues (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) . 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. 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 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 log 2 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 tables S1-S5 are available at Molecular Biology and Evolution online (http://www.mbe.oxfordjournals.org/). The histone variant H3.3 marks active chromatin by replication-independent nucleosome assembly Cell nonautonomy of C. elegans daf-2 function in the regulation of diapause and life span Comparison of sequencing utility programs Latitudinal variation of wing: thorax size ratio and wing-aspect ratio in Drosophila melanogaster Gene expression during Drosophila melanogaster egg development before and after reproductive diapause Seasonal variation in life history traits in two Drosophila species Controlling the false discovery rate: a practical and powerful approach to multiple testing Genomic evidence of rapid and stable adaptive oscillations over seasonal time scales in Drosophila Secondary contact and local adaptation contribute to genome-wide patterns of clinal variation in Drosophila melanogaster Molecular analysis of an allozyme cline: alcohol dehydrogenase in Drosophila melanogaster on the east coast of North America Life history adaptation along a latitudinal cline in the water strider Aquarius remigis (Heteroptera: Gerridae) Autonomous control of cell and organ size by CHICO, a Drosophila homolog of vertebrate IRS1-4 Large clusters of co-expressed genes in the Drosophila genome Fluctuating viability selection on morphology of cliff swallows is driven by climate Phenotypic and genetic variability of morphometrical traits in natural populations of Drosophila melanogaster and D. simulans. I. Geographic variations Extension of life-span by loss of CHICO, a Drosophila insulin receptor substrate protein Latitudinal cline of requirement for far red light for the photoperiodic control of budset and extension growth in Picea abies (Norway spruce) The intensity of selection acting on couch potato gene-spatial-temporal variation in a diapause cline Variation in Drosophila melanogaster central metabolic genes appears driven by natural selection both within and between populations Sequence and transcription patterns of 60S ribosomal protein P0, a diapause-regulated AP endonuclease in the flesh fly, Sarcophaga crassipalpis Genetic variation of Drosophila melanogaster natural populations Regulation of diapause Analysis of the host transcriptome from demyelinating spinal cord of murine coronavirus-infected mice Genome-wide patterns of latitudinal differentiation among populations of Drosophila melanogaster from North America Spatially varying selection shapes life history clines among populations of Drosophila melanogaster from sub-Saharan Africa Diapause-specific gene expression in pupae of the flesh fly Sarcophaga crassipalpis Life-history evolution and the polyphenic regulation of somatic maintenance and survival Clinal variation in developmental time and viability, and the response to thermal treatments in two species of Drosophila DNA sequence variation and latitudinal associations in hsp23, hsp26 and hsp27 from natural populations of Drosophila melanogaster Nucleotide diversity at two phytochrome loci along a latitudinal cline in Pinus sylvestris Long-lived Drosophila with overexpressed dFOXO in adult fat body Latitudinal variation in starvation resistance is explained by lipid content in natural populations of Drosophila melanogaster Unpredictable evolution in a 30-year study of Darwin's finches Molecular evolution and functional characterization of Drosophila insulin-like peptides Fluctuating selection and response in a population of freshwater copepods A maternal effect that eliminates pupal diapause in progeny of the flesh fly, Sarcophaga bullata Opposing clines for high and low temperature resistance in Drosophila melanogaster Diet controls normal and tumorous germline stem cells via insulin-dependent andindependent mechanisms in Drosophila Rapid evolution of a geographic cline in size in an introduced fly Drosophila dFOXO controls lifespan and regulates insulin signaling in brain and fat body Genetic and environmental responses to temperature of Drosophila melanogaster from a latitudinal cline Regulated chromatin domain comprising cluster of co-expressed genes in Drosophila melanogaster Population genomic analysis uncovers African and European admixture in Drosophila melanogaster populations from the south-eastern United States and Caribbean Islands Inference of chromosomal inversion dynamics from Pool-Seq data in natural and laboratory populations of Drosophila melanogaster Desiccation and starvation tolerance of adult Drosophila: opposite latitudinal clines in natural populations of three different species daf-2, an insulin receptor-like gene that regulates longevity and diapause in Caenorhabditis elegans Genetic dissection of the Drosophila melanogaster female head transcriptome reveals widespread allelic heterogeneity Ovarian development in Drosophila melanogaster Chromosome inversion polymorphisms in Drosophila melanogaster II. Geographic clines and climatic associations in Australasia Fitness effects associated with the major flowering time gene FRIGIDA in Arabidopsis thaliana in the field The sleeping beauty: how reproductive diapause affects hormone signaling, metabolism, immune response and somatic maintenance in Drosophila melanogaster Direct control of germline stem cell division and cyst growth by neural insulin in Drosophila Drosophila short neuropeptide F signalling regulates growth by ERK-mediated insulin signalling Polymorphism in the neurofibromin gene, Nf1, is associated with antagonistic selection on wing size and development time in Drosophila melanogaster EBSeq: an empirical Bayes hierarchical model for inference in RNA-seq experiments Whole-genome expression plasticity across tropical and temperate Drosophila melanogaster populations from eastern Australia RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome Genome Project Data Processing Subgroup. 2009. The Sequence Alignment/Map format and SAMtools Spatial and temporal expression of the period gene in Drosophila melanogaster Variation in body size and life history traits in Drosophila aldrichi and D. buzzatii from a latitudinal cline in eastern Australia Cutadapt removes adapter sequences from highthroughput sequencing reads A clinally varying promoter polymorphism associated with adaptive variation in wing size in Drosophila Evolutionary links between circadian clocks and photoperiodic diapause in insects Coexpression, coregulation, and cofunctionality of neighboring genes in eukaryotic genomes Postponed reproduction as an adaptation to winter conditions in Drosophila melanogaster: evidence for clinal variation under semi-natural conditions Notes on continuous stochastic phenomena Genome-wide analysis of self-renewal in Drosophila neural stem cells by transcgenic RNAi Different maternal effects on diapause induction of tetranychid mites, Tetranychus urticae and T. kanzawai (Acari: Tetranychidae) The regulation of histone synthesis in the cell cycle A highly pleiotropic amino acid polymorphism in the Drosophila insulin receptor contributes to life-history adaptation Identification of a candidate adaptive polymorphism for Drosophila life history by parallel independent clines on two continents APE: Analyses of Phylogenetics and Evolution in R language Role for circadian clock genes in seasonal timing: testing the B€ unning hypothesis Independent photoreceptive circadian clocks throughout Drosophila Deep sequencing reveals complex mechanisms of diapause preparation in the invasive mosquito, Aedes albopictus RNA-Seq reveals early distinctions and late convergence of gene expression between diapause and quiescence in the Asian tiger mosquito, Aedes albopictus Population genomics of sub-Saharan Drosophila melanogaster: African diversity and non-African admixture Selective maintenance of Drosophila tandemly arranged duplicated genes during evolution Mechanisms of suspended animation are revealed by transcript profiling of diapause in the flesh fly Clinal variation in fitness related traits in tropical drosophilids of the Indian subcontinent Body melanization and its adaptive role in thermoregulation and tolerance against desiccating conditions in drosophilids Variations in body melanisation, ovariole number and fecundity in highland and lowland populations of Drosophila melanogaster from the Indian subcontinent Parallel geographic variation in Drosophila melanogaster Transcriptional evidence for small RNA regulation of pupal diapause in the flesh fly, Sarcophaga bullata Ecdysteroids regulate yolk protein uptake by Drosophila melanogaster oocytes RNA-Seq and find: entering the RNA deep field Starvation resistance and adult body composition in a latitudinal cline of Drosophila melanogaster A diapause maternal effect in the flesh fly, Sarcophaga bullata: transfer of information from mother to progeny Maternal influence on the incidence and duration of larval diapause in Calliphora vicina Insect photoperiodism: measuring the night Regulation of ovarian diapause in Drosophila melanogaster by photoperiod and moderately low temperature Induction of diapause in Drosophila melanogaster: photoperiodic regulation and the impact of arrhythmic clock mutations on time measurement Evolution and mechanisms of insect reproductive diapause: a plastic and pleiotropic life history syndrome Environmental heterogeneity and the maintenance of genetic variation for reproductive diapause in Drosophila melanogaster Geographic variation in diapause incidence, life-history traits, and climatic adaptation in Drosophila melanogaster Reproductive diapause and life-history clines in North American populations of Drosophila melanogaster Genetic variance for diapause expression and associated life histories in Drosophila melanogaster An amino acid polymorphism in the couch potato gene forms the basis for climatic adaptation in Drosophila melanogaster Single-locus latitudinal clines and their relationship to temperate adaptation in metabolic genes and derived alleles in Drosophila melanogaster Insulin signaling and FOXO regulate the overwintering diapause of the mosquito Culex pipens Insulin signaling and the regulation of insect diapause Identification of FOXO targets that generate diverse features of the diapause phenotype in the mosquito Culex pipiens Evidence for large domains of similarly expressed genes in the Drosophila genome A latitudinal cline in flowering time in Arabidopsis thaliana modulated by the flowering time gene FRIGIDA Opposing selection and environmental variation modify optimal timing of breeding Negligible senescence during reproductive dormancy in Drosophila melanogaster A mutant Drosophila insulin receptor homolog that extends life-span and impairs neuroendocrine function Slow aging during insect reproductive diapause: why butterflies, grasshoppers and flies are like worms Natural selection favors a newly derived timeless allele in Drosophila melanogaster Seasonal adaptations of insects Genomic analysis of adaptive differentiation in Drosophila melanogaster A rapid shift in a classic clinal pattern in Drosophila reflecting climate change Clinal variation for amino acid polymorphisms at the Pgm locus in Drosophila melanogaster Natural variation in Drosophila melanogaster diapause due to the insulin-regulated PI3-kinase Diapause in Drosophila melanogaster females: a genetic analysis TOR signaling in growth and metabolism Juvenile hormone regulation of Drosophila aging Key molecular processes of the diapause to postdiapause quiescence transition in the alfalfa leafcutting bee Megachile rotundata identified by comparative transcriptome analysis Gene ontology analysis for RNA-seq: accounting for selection bias Association of the circadian rhythmic expression of GmCRY1a with a latitudinal cline in photoperiodic flowering of soybean We would like to thank Mia Levine for helpful discussion on chromosome location-based gene regulation and cell cycle arrest. We would like to thank Peter Petraitis for discussion on statistical methods to test for spatial autocorrelation. Two anonymous reviewers have offered helpful suggestions that improved the quality of the paper. This work was supported by a National Science Foundation grant to P.S. (DEB 0921307) and a National Institute of Health grant to D.P. and P.S. (R01 GM100366).