Transcriptome profiling reveals mechanisms for the evolution of insect seasonality

Defined as a state of environmentally induced developmental arrest, dormancy is a physiologically dynamic developmental trajectory and is a nearly ubiquitous life-history strategy for animals and plants to survive inhospitable climatic conditions in temperate and polar environments. Diapause, a widespread form of insect dormancy, is also viewed as the primary mechanism through which the annual rhythm of major life history phases such as reproduction, growth, development and migration are synchronized with seasonally varying biotic and abiotic requirements (Danilevskii, 1965; Tauber and Tauber, 1981; Danks, 1987; Leopold and Lang, 1996). Diapause occurs at various life-stages in insects (egg, larval, pupal, adult) and proceeds through pre-diapause, diapause and post-diapause developmental phases (Koštál, 2006). Diapause itself can be further divided into an initiation phase in which direct development ceases, a maintenance phase in which developmental arrest continues and a termination phase in which a stimulus triggers an active potential for development although development may not be overtly visible (Tauber et al., 1986; Koštál, 2006).

Species and populations can reveal subtle shifts in the timing of the unique stages of diapause that can substantially alter annual life-history patterns. Whereas the onset of diapause commonly initiates seasonal migration or polyphenism at appropriate times, the termination of diapause synchronizes active growth and reproduction with favorable conditions (Tauber et al., 1986). Examples of rapid shifts in the timing of the onset or termination of diapause are common, with many believed to be adaptive responses to a range of contemporary environmental perturbations, including changing climates, anthropogenic impacts and species introductions (Bradshaw and Holzapfel, 2001, 2006, 2008; Bradshaw et al., 2004; Filchak et al., 2000; Mathias et al., 2005; Schmidt et al., 2005; Gomi et al., 2007).

To understand how organisms adapt and persist in novel seasonal environments like those experienced during global warming, we must understand the molecular bases for shifts in diapause timing. In several study systems, different populations or species show shifts to their seasonal dormancy timing (Pickett and Neary, 1940; Tauber and Tauber, 1976; Bradshaw and Lounibos, 1977; Glover et al., 1992). However, complications in these systems arising from genetic architecture (e.g. polygenic, epistasis) or genomic architecture (e.g. chromosomal rearrangements), have confounded the identification of genes underlying variation in diapause timing (Tauber et al., 1977; Feder et al., 2002; Bradshaw et al., 2005; Mathias et al., 2007; Emerson et al., 2010; Wadsworth et al., 2015). Despite these difficulties, genes involved in the transitions of diapause phases have been identified. Cell cycling, DNA replication/transcription and stress response genes have been proposed to be important for diapause initiation in moths and mosquitos (Bao and Xu, 2011; Poelchau et al., 2011). Additionally, genes in the Wnt and TOR signaling pathways have been nominated as strong regulatory candidates during diapause termination in Rhagoletis flies (Ragland et al., 2011). Perhaps the best single candidate thus far identified is for the egg diapause of the silkmoth (Bombyx mori), in which transgenerational diapause induction involves the temperature-sensitive Transient receptor potential A1 (Trpa1) gene (Sato et al., 2014). Nevertheless, in many of these systems it is still unclear how the regulation of diapause genes and pathways are altered to produce ecologically relevant variation in diapause timing.

A useful model for studying the molecular basis of shifts in seasonality is the European corn borer moth (Ostrinia nubilalis). Introduced to North American from Europe in the early 20th century, multiple corn borer ecotypes exist that differ by 20–30 days in the time needed for overwintering larvae in diapause to reinitiate development in spring (Glover et al., 1992; Dopman et al., 2005). Over a half century of research has connected variation in diapause timing to latitudinal differences in generation number (Beck and Apple, 1961; Palmer et al., 1985; Levy et al., 2015). More recently, differences in diapause timing have been emphasized because they contribute to asynchronous adult mating flights of two-generation E-strain borers and one-generation Z-strain borers in New York (Dopman et al., 2010). Thus, divergence in the timing of diapause appears relevant for both successful colonization and ecological speciation.

Two observations provide important insight about possible loci responsible for diapause timing change in European corn borers. First, the major genetic factor Pdd (post-diapause development) has been shown to map to a region of the Z (sex) chromosome that may be rearranged between corn borer strains (Glover et al., 1992; Dopman et al., 2004, 2005; Wadsworth et al., 2015). Second, high temperatures and/or long-day photoperiod cues are known to stimulate diapause termination in corn borers (McLeod and Beck, 1963). Thus, pathways and genes on the sex chromosome that are sensitive to light and temperature represent strong candidate genes for Pdd.

In a prior study, metabolic rate trajectories were used to identify a shift in the timing of the end of the diapause termination phase as critical for explaining earlier emergence in the E-strain and later emergence in the Z-strain of the European corn borer moth (Fig. 1A) (Wadsworth et al., 2013). Within a week of exposure to diapause-breaking (long-day) conditions, E-strain moths rapidly progress through the diapause termination phase and enter post-diapause development, as indicated by increased respiratory metabolism. In contrast, Z-strain borers exhibit continued suppression of metabolic rate for nearly a month longer, suggesting that timing divergence stems from an extended diapause termination phase. Here, we sample from the developmental time course of diapause in both strains and we profile regulatory and amino acid changes coincident with changes in timing. We then address two main questions. First, what candidate genes might underlie Pdd and therefore divergence in diapause timing? Second, what downstream pathways might direct insects to an earlier or later end of diapause?

Diapause was induced in larvae of the European corn borer moth (Ostrinia nubilalis Hübner 1825) using short-day conditions (12 h:12 h light:dark at 23°C). After 35 days of exposure to these conditions, diapausing borers remain at the 5th instar whereas direct developing insects will have already pupated or enclosed as adults (Glover et al., 1992; Dopman et al., 2005; Wadsworth et al., 2013). Diapause was terminated in 36-day-old larvae using long-day conditions (16:8 LD at 26°C) (Glover et al., 1992; Dopman et al., 2005; Wadsworth et al., 2013). Alternative developmental genetic pathways underlying differences in diapause termination timing seem to be triggered by Pdd within the first 7 days after insects experience diapause-breaking cues, as indicated by elevated respiratory metabolism in the E-strain but not the Z-strain (Wadsworth et al., 2013). Therefore, we collected tissues at three time points from both strains within this interval (Fig. 1B): (1) in diapause maintenance prior to a long-day (diapause breaking) cue; (2) 1 day after transfer to long-day conditions when strains are physiologically indistinguishable but have entered the diapause termination phase; and (3) on day 7 after transfer when the Z-strain remains in diapause termination (previously referred to as ‘diapause maintenance’ in Wadsworth et al., 2013) and most E-strain moths have entered post-diapause development. Head tissues were sampled to capture changes in major neuroendocrine glands and the primary endogenous clock (Tauber et al., 1986; Koštál, 2011). Five individuals were pooled for each of five replicates per strain per time point. All samples were collected at the same time of day (14:00 h–16:00 h) to avoid spurious effects of circadian rhythmicity. Tissues were ground in liquid nitrogen and total RNA was prepared for multiplexed Illumina sequencing using the Qiagen RNeasy kit (Qiagen Inc., Germantown, MD, USA) followed by the TruSeq RNA Prep kit (Illumina Inc., San Diego, CA, USA). Six libraries were multiplexed per lane and subjected to single-end 50 bp sequencing on an Illumina HiSeq 2000 (Illumina) at Tufts University Core Facility.

Low-quality reads and adapter sequences were removed using Trimmomatic v0.27 (Lohse et al., 2012). Retained reads were ≥36 bp after eliminating reads with an average quality of <15 and trimming off ends with a quality of <5 in leading and trailing bases. Mitochondrial (NC_003367.1) and ribosomal contaminants (all Ostrinia ribosomal sequences in NCBI) were removed using Bowtie2 v2.1.0 (Langmead and Salzberg, 2012). Library quality was subsequently assessed using FastQC v0.10.1 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Prior to de novo assembly, duplicate reads were collapsed to reduce library complexity using FASTX-Toolkit v0.0.13 (Gordon and Hannon, 2010). Trinity was chosen for de novo assembly because it has been shown to recover more full-length genes than other single-kmer de novo assembly programs (Grabherr et al., 2011; Zhao et al., 2011). It was run at default parameters (kmer=25, minimum contig length=48) on all libraries and the longest transcript per transcript cluster was used for subsequent analyses. Annotation of genes was by TBLASTX search to FlyBase and the NCBI non-redundant (nr) database and retention of transcripts with E-values ≤10⁻⁵. Putative locations for corn borer transcripts was determined using TBLASTX (E-value ≤10⁻⁵) and the genome of B. mori (KAIKObase v3.2.2) (Shimomura et al., 2009), which shows conserved macrosynteny with O. nubilalis (Kroemer et al., 2011). The Pdd QTL is known to be between Ket and Ldh (Dopman et al., 2005) or between 6.5 and 17.3 Mb on the B. mori Z chromosome. We used Bowtie 2 v2.1.0 with ‘end-to-end’ and ‘very-sensitive’ options, which resulted in a higher percentage of uniquely mapped reads than any other combination of preset options. Count data were generated using the single best alignment for each read. To help mitigate G-C bias during PCR amplification, within and between sample data were corrected for G-C content using EDAseq (Fig. S2) (Risso et al., 2011).

Annotation of genes was conducted by TBLASTX search to both the FlyBase and NCBI nr databases, retaining annotations with E-values ≤10⁻⁵ as putative orthologs. GOstats v2.12 conditional hypergeometric test (Falcon and Gentleman, 2007) was used to test for ‘Biological Process’, ‘Molecular Function’ and ‘KEGG pathway’ enrichment using the ontology terms from the Bioconductor package org.Dm.eg.db v2.10.1 with the P-value cut-off of ≤0.001. Short Time-series Expression Miner v1.3.8 (Ernst et al., 2005) was used to search for groups of genes that showed coordinated expression change from diapause-maintaining to diapause-breaking conditions. Normalized expression data were randomly coupled from diapause maintenance to later developmental time points and then log₂ normalized to simulate microarray M-values. Multiple transcripts with duplicate annotations were averaged. STEM Miner determines all possible iterations gene expression change through time, with the diapause maintenance time point normalized to 0, increasing or decreasing up to a user defined step size (=2). Each gene in the dataset was then assigned to a given model based on its correlation coefficient, which was defined as 1. P-values were determined by calculating the observed number of genes assigned to each profile compared with the expected number, which is calculated by permuting the fit of each gene to the model 1000 times (P<0.01). Significantly enriched gene ontology terms for each time-series were determined using Bonferroni corrected P-values (≤0.0001).

Differential expression analyses used edgeR (Robinson and Oshlack, 2010), which employs an empirical Bayes method to estimate gene-specific biological variation. A multi-factor analysis was used to account for both day and strain factors. Per-gene dispersions were estimated using the Cox–Reid profile-adjusted method and a negative binomial generalized log-linear model was fit to identify significantly differentially expressed transcripts (Figs S3 and S4). Subsequently, gene lists were corrected for multiple testing using a false discovery rate (FDR) (cut-off of ≤0.01). Candidates within the predicted QTL interval were selected if they showed no differential expression between strains during diapause maintenance, were differentially expressed between strains on day 1 or 7, or were differentially expressed within strains through time.

We identified fixed single nucleotide polymorphisms (SNPs) between strains by pooling all libraries by strain, mapping reads and estimating SNP frequencies using Popoolation2 v1.201 (Kofler et al., 2011). This method has been shown to be effective at recovering the majority of true allele frequencies in populations (Konczal et al., 2014). A minimum coverage of 8 reads per population and a pool size of 75 were used to calculate F_ST  (a measure of genetic variation between populations, with a value of 1 indicating complete divergence and a value of 0 indicating panmixis) for every SNP. To determine if fixed SNPs coded for an amino acid change in the subset of transcripts whose putative locations were within the QTL, ESTScan was used to predict the reading frame (Iseli et al., 1999).

Libraries were archived in GenBank (Bioproject ID: PRJNA294976; for accession numbers see Table S4).

For Z- and E-strain corn borers, we sampled three time points along the photoperiodically initiated and terminated diapause developmental pathway (Fig. 1B). Larval diapause in corn borers is initiated under short-day photoperiod, even under temperatures that are permissive for development (McLeod and Beck, 1963; Glover et al., 1992). Therefore, the first sampled time point occurred when both strains were in the diapause maintenance phase under short-day diapause inducing conditions (12 h:12 h light:dark). Exposure to long-day conditions stimulates entrance of corn borer larvae into the diapause termination phase and leads to downstream pupal development (McLeod and Beck, 1963; Glover et al., 1992). Therefore, the second sampled time point occurred on day 1 after long-day exposure (16 h:8 h light:dark), but before a physiological change in respiratory metabolism in either strain (hereafter ‘day 1’). The final time point was on day 7 after long-day exposure and after physiological divergence between strains (hereafter ‘day 7’). Within this sampling scheme, we expected changes on day 1 to be more reflective of strain-specific differences in upstream causal factors. Day 7, by contrast, might be more informative of downstream developmental consequences of Pdd such as the transition from diapause termination to post-diapause development (as indicated by a metabolic uptick, Fig. 1A). Nevertheless, within-strain variation in the length of the diapause termination phase exists (Wadsworth et al., 2013) and therefore patterns of gene expression on day 1 and day 7 probably represent a mixture of cause and effect.

As the most likely tissues involved in photoperiodic responses have been hypothesized to be located in the brain, eyes or major neuroendocrine glands (Koštál, 2011), at each of the three time points we sampled larval head tissues (N=5 individuals for n=5 replicates per strain). A total of ∼562 million 50 bp single-end reads were generated across 30 libraries. Each library had, on average, 19.46±8.99 million reads. After removal of adapters, mtDNA and rRNA contaminants, average library size decreased to 18.38±8.24 million reads. Our de novo assembly consisted of 47,565 transcripts (mean length=686, N50=1025). Of these, 14,177 were annotated via TBLASTX to FlyBase and 23,077 were annotated via TBLASTX to the NCBI nr database.

There were a total of 15,208 transcripts that were significantly differentially expressed in at least one comparison involving strains or time points (Table 1). We removed some transcripts from our analysis to minimize the number of genes that were differentially expressed between strains because of population divergence (strains diverged approximately 100,000 years ago, Malausa et al., 2007) and not because of associations with differences in the timing of diapause termination. Because E- and Z-strain corn borers are physiologically indistinguishable under short-day diapause-inducing conditions (Wadsworth et al., 2013), we considered differentially expressed transcripts between strains at this time point as a control for regulatory divergence unrelated to photoperiodic response or divergence in length of the diapause termination phase. We therefore removed these transcripts from differential expression comparisons between strains at other time points (day 1 and day 7) to help control for false positives.

To nominate transcripts as candidates for the 30 day shift in the end of diapause, we initially considered genomic position. Pdd maps to a region of the European corn borer Z chromosome between Ket and Ldh (Dopman et al., 2005). Therefore, an initial candidate list was nominated by TBLASTX hits of corn borer transcripts to the macrosyntenic interval of the B. mori Z chromosome containing these genes (6.5–17.4 Mb) (Kroemer et al., 2011). Of the corn borer transcripts that mapped to this interval, two forms of evidence were evaluated. First, the causal factor(s) might be regulatory. Therefore, patterns of differential expression either within strain from diapause maintenance to a later developmental time point, or between strains on day 1 or day 7, were considered. Second, as differences in diapause timing could be due to structural mutation(s), we considered the presence of fixed amino-acid-changing substitutions between strains.

Of the 563 (∼1%) transcripts with putative locations between Ket and Ldh, a total of 48 were nominated as candidates for Pdd. Of these, 41 genes were nominated based on patterns of significant differential expression between strains, or significant differential expression within strains across the developmental time course (Fig. 2; Table 2). Six broad patterns of expression differences (A–F) were observed between strains: (A) five genes were upregulated in the Z-strain on day 1 relative to the E-strain; (B) three genes were upregulated in the Z-strain on day 7 relative to the E-strain; (C) three genes were upregulated in the Z-strain on day 1 and day 7 relative to the E-strain; (D) ten genes were upregulated in the E-strain on day 1 relative to the Z-strain; (E) eight genes were upregulated in the E-strain on day 7 relative to the Z-strain; and (F) four genes were upregulated in the E-strain on days 1 and 7 relative to the Z-strain. In addition, within the E-strain we observed two further expression patterns across the time course: (G) four downregulated genes from diapause maintenance to day 1; and (H) one upregulated gene from diapause maintenance to day 7. Of the transcripts nominated by regulatory patterns, all had some associated annotation information (Table 2).

We next considered candidates in the Ket–Ldh interval that harbored fixed amino acid changes between Z- and E-strains. Of the total number of variable SNPs between strains (N=1,347,090), 250 were located in European corn borer transcripts with significant TBLASTX hits to the Ket–Ldh interval in B. mori. These were distributed among 114 unique transcripts, 11 of which contained at least one predicted amino acid change that was fixed between strains. Eight of the 11 transcripts had annotation information (Table 3). Four amino acid changes were switches in polarity, two were switches in polarity and charge and one was a switch in charge.

To characterize regulatory mechanisms that might represent consequences of differences at Pdd, we focused on differential expression of transcripts between strains after exposure to diapause-breaking cues (Table 1). For gene ontology (GO) and KEGG (Kyoto encyclopedia of genes and genomes) enrichment (P-value cut-off ≤0.001) of differentially expressed gene lists, on day 1 there were 14 enriched ‘biological process’ GO terms, four ‘molecular function’ terms and five KEGG pathways (Fig. 3; Table S1). Two terms associated with Wnt signaling were also enriched (P≤0.05). Between strains on day 7 there were 24 enriched ‘biological process’ terms, nine ‘molecular function’ terms and one KEGG pathway (Fig. 3; Table S1).

To identify pathways under the control of similar regulatory mechanisms, for each strain we evaluated correlated changes in gene expression across the three time points. Normalized expression data from the diapause maintenance time point were randomly coupled to later developmental time points and then log₂ normalized to simulate microarray M-values. Short time-series expression miner (STEM) (Ernst et al., 2005) was then used to discover five significant gene expression profiles in the E-strain and two in the Z-strain (P≤0.01) (Fig. 4). Pattern A had 105 transcripts assigned to it, pattern B had 109, pattern C had 108, pattern D had 67, pattern E had 64, pattern F had 38 and patterns G had 32. The most extreme expression changes were seen from diapause maintenance to day 1. Groups of genes with similar expression patterns across the time course in the E-strain were enriched for GO categories including development, metabolism, protein catabolism and oxidoreductase activity (Table S2). In contrast, there were no significantly enriched GO terms for the two temporal expression profiles in the Z-strain. Time series expression profiles were also plotted for individual transcripts in pathways of interest. Pathways included cell cycling, ecdysone, circadian rhythms, heat shock, Wnt signaling and insulin signaling (Fig. 5; Fig. S1). At least one transcript within each of these pathways showed significant differential expression between the strains (Table S3).

In insects, divergence in seasonal life-cycle timing may commonly involve changes in the duration of the diapause termination phase in spring. Z- and E-strain European corn borers differ by ∼30 days in the length of the diapause termination phase as indicated by a shift in the timing of metabolic release coincident with the end of diapause (Wadsworth et al., 2013) (Fig. 1A). Although strains are physiologically indistinguishable in early stages of diapause termination, RNA sequencing shows this to be a period of intense transcriptional activity in the E-strain, as indicated by a large number of differentially expressed transcripts from diapause maintenance to day 1 (Table 1). Conversely, the Z-strain shows negligible change from diapause maintenance to day 1 or 7 (Table 1). Thus, the major genetic cause of altered seasonal diapause timing, Pdd, appears to act as a genetic switch operating shortly after the beginning of the growing season in the E-strain, allowing for rapid development, while being delayed in the Z-strain.

Specific gene candidates for shifts in insect seasonal timing have been difficult to accumulate via QTL mapping because of complex genetic bases or genomic architectures between populations (Tauber et al., 1977; Feder et al., 2002; Bradshaw et al., 2005; Mathias et al., 2007; Wadsworth et al., 2015). In such systems, transcriptome profiling can be a useful alternative to QTL mapping to nominate candidate genes. For example, in the pitcher-plant mosquito (Wyeomyia smithii), northern and southern populations differ in the critical day length required for diapause termination, but the genetic basis is complex, involving many loci and epistasis (Bradshaw et al., 2005; Mathias et al., 2007). Using microarrays, a total of 29 genes were nominated as candidates with one, a cuticular protein with cornea-specific expression, falling within a previously defined QTL for altered termination response (Emerson et al., 2010).

In European corn borers, early results implicated a major mendelian factor underlying differences in diapause termination timing (Glover et al., 1992; Dopman et al., 2004, 2005). However, recent findings show that Pdd and roughly 20% of the sex chromosome may be caught up in a large inversion (Wadsworth et al., 2015). As inversions can cause large swaths of genes to be inherited together and thus act as ‘supergenes’, there is potential for Pdd to be composed of multiple linked loci. Of the genes that are predicted to be within the chromosomal interval containing Pdd, 48 have fixed amino acid changes between strains or experience a change in expression through diapause termination (Tables 2 and 3). Candidate gene lists did not overlap between W. smithii and European corn borers (Emerson et al., 2010). Within our nominated gene lists, transcripts involved in insulin signaling and circadian rhythms are some of the most promising candidates for Pdd because both pathways are sensitive to light and temperature (Suttie et al., 1991; Chan et al., 1999; Shingleton et al., 2005; Taylor et al., 2005; Luckenbach et al., 2007; Sim and Denlinger, 2009), environmental factors that are known to stimulate diapause termination in corn borers (McLeod and Beck, 1963; Glover et al., 1992). These two pathways are also proposed regulators of dormancy or diapause in other species (Ikeno et al., 2010; Hahn and Denlinger, 2011). Of the 48 candidates, we focus on four with direct involvement in these pathways.

The γ subunit (Snf4Aγ, 15.6 Mb in B. mori) of AMP-activated protein kinase (AMPK) is a highly conserved gene that makes up the energy-sensing component of the AMPK protein (Kahn et al., 2005). AMPK is a member of the insulin, target of rapamycin (TOR) and cell death pathways (Kahn et al., 2005; Lippai et al., 2008), all of which are known to be important in the progression of development and molting in Lepidoptera (Lippai et al., 2008; Gu et al., 2009; Smith et al., 2014). In our transcriptome, Snf4Aγ was significantly downregulated in the later-emerging Z-strain compared with the earlier-emerging E-strain (FDR<0.001 on day 7) (Table 2; Fig. 2). In Drosophila, loss of Snf4Aγ expression has been shown to inhibit the ability of 20-hydroxyecdysone (20HE) to promote cell death, a major component in metamorphosis in insects (Lippai et al., 2008). Thus, perhaps suppression of Snf4Aγ in Z-strain moths promotes a state of delayed development through inhibition of cell death and turnover, processes that occur ∼20–30 days earlier in the E-strain.

Two members of the insulin growth factor (IGF) pathway were identified as candidate genes. The first, IGF-II mRNA binding protein (Imp, 11.4 Mb in B. mori), is associated with growth and development and, like the insulin signaling pathway, the IGF pathway commonly shows sensitivity to both temperature and photoperiod (Suttie et al., 1991; Dahl et al., 1997; Luckenbach et al., 2007). Imp was found to be downregulated (day 1 and 7) in the later-emerging Z-strain larvae (FDR<0.001) (Table 2; Fig. 2). This result is notable because suppression of the IMP homolog in mice prevents IMP from binding to insulin growth factors, resulting in a reduction in growth (Hansen et al., 2004). In both mice and nematodes, disruption of other genes within the IGF pathway result in phenotypes reminiscent of longer dormancy such as growth reduction and extend lifespan (Liu et al., 1993; Kenyon, 2010). Hence, downregulation of Imp in corn borers may promote later emergence from diapause in the Z-strain through decreased IGF signaling.

The second IGF candidate, Receptor for activated C kinase 1  (Rack1, 10.9 Mb in B. mori), interacts with the IGF pathway through the insulin-like growth factor I receptor (IGF-IR). In corn borers, Rack1 was upregulated in the later-emerging Z-strain on day 7 (FDR<0.001) (Table 2; Fig. 2). Rack1 shows a similar association with diapause maintenance in the Asian tiger mosquito (Aedes albopictus), where it was one of 314 transcripts to be upregulated upon exposure to short-day conditions and entrance into diapause (Poelchau et al., 2011). Overexpression of RACK1 protein has also been shown to lead to reduced cell growth in mammals (Hermanto et al., 2002). With the upregulation of Rack1 in the Z-strain and known upregulation in other species delaying growth or promoting diapause, this gene is an interesting candidate for Pdd.

A final gene candidate for Pdd is the clock gene Period (Per, 13.0 Mb in B. mori). Per shows a predicted glycine to arginine amino acid change between Z- and E-strain insects (Table 3) and is modestly upregulated in the E-strain on both day 1 and day 7 (FDR=0.58 and 0.89, respectively) (Fig. 5A,B). Evidence is accumulating that clock genes like Per are important for photoperiodic responses and diapause transitions. Within the Z-strain of corn borers, a different amino acid substitution in Per varies with diapause emergence timing and latitudinal changes in generation number (Levy et al., 2015). In D. melanogaster, geographic variation in a threonine-glycine repeat within Per associates with diapause incidence (Kyriacou et al., 2008). Laboratory studies also support a role for Per in diapause. RNA interference of Per in the bean bug (Riptortus pedestris) disrupts both circadian rhythms and photoperiodic response (Ikeno et al., 2010), and artificial selection for enhanced diapause in grey flesh fly (Sarcophaga bullata) is associated with a deletion at the Per locus that is 33 amino acids long (Han and Denlinger, 2009).

Our results show that the diapause termination phase and the transition into post-diapause development involves stress response, hormone processing, cell cycling, developmental patterning and metabolic pathways (Figs 3 and 4). Similar results were found in flies (Emerson et al., 2010; Poelchau et al., 2011; Ragland et al., 2010, 2011). This suggests a conserved transcriptional basis for the end of the diapause developmental pathway, despite ∼200 million years of divergence (Pringle et al., 2007) and distinct diapausing life-stages (larval versus pupal) in different environments (host plant versus soil).

Previous results in flies suggest that a common characteristic of diapause is cell-cycle arrest in the G0/G1 or G2 phases, with transition into cell cycle progression being indicative of the end of diapause (Koštál et al., 2009; Ragland et al., 2010, 2011). Gene expression profiles in corn borers support a transition into active cell cycling upon exposure to long-day cues in the E-strain but not the Z-strain (Fig. 5C,D), suggesting E-strain borers are rapidly exiting diapause. For example, proliferating cell nuclear antigen (Pcna), a DNA-replication-related protein that is active in the S phase of the cell cycle, is upregulated during the rapid diapause termination phase in E-strain European corn borers and also diapause termination of Rhagoletis flies (Ragland et al., 2011), but it is downregulated in Z-strain borers (FDR<0.00001 on day 1). Differentially expressed genes between Z- and E-strains were highly enriched for cell cycling GO terms (Fig. 3; Table S1). Additionally, genes involved in the progression of the cell cycle were upregulated in the E-strain, whereas they were downregulated in Z-strain [e.g. CycA, CycB (FDR<0.001 days 1 and 7); Cdk4, Polo (FDR<0.01 on day 1)]. Overall, patterns suggest that active cell cycling occurs within days of exposure to long-day conditions in the E-strain allowing for release from diapause-induced metabolic depression and progression into post-diapause development, whereas cellular development remains suppressed for ∼30 additional days in the Z-strain.

Transcripts involved in development showed congruent patterns to those involved in mitotic progression. Wnt signaling pathways promote proper developmental patterning in insects (e.g. metamorphosis) through cellular communication, resumption of the cell cycle and also may regulate growth through the canonical growth-regulatory pathways (Logan and Nusse, 2004; Gokhale and Shingleton, 2015). The Wnt signaling and canonical Wnt signaling pathways showed enrichment of GO terms on day 1 between strains (Fig. 3; Table S1). Many transcripts involved in these pathways were downregulated in the Z-strain and upregulated in the E-strain [e.g. Fz3, Arm, Smo, Wg and Pan (FDR<0.01 on day 1)] (Fig. S1A,B; Table S3). Coupled with the regulation of members of the cell cycling pathway, upregulation of members of Wnt signaling pathways in the E-strain indicate a reversal of arrest on both the cellular and developmental levels.

Studies in flies have documented shifts from anaerobic to aerobic metabolism during the diapause termination phase as the hypoxic stress of winter is lifted (Emerson et al., 2010; Ragland et al., 2010, 2011). In moths, we find evidence for an equivalent shift in the E-strain but not in the Z-strain. Members of the anaerobic pyruvate metabolic pathway experienced a ∼2-fold  decrease in expression within the E-strain from diapause maintenance to day 1 or day 7 [e.g. Aldh, CG6084, CG9629, Men, CG31674, Hex-t2 (FDR<0.001)] and the pathway was upregulated in the Z-strain compared with the E-strain during diapause termination [e.g. Aldh, CG9629, CG31075 (FDR<0.01 on day 1)]. These results suggest that a switch to aerobic metabolism is not associated with the diapause termination phase in moths per se, but is primarily associated with the metabolic uptick that marks the end of the diapause termination phase.

Diapause is commonly associated with cold-hardiness across many insect species, yet it is unclear if the cold-hardy state persists through all of the phases of diapause. Our results suggest that in moths, the reduction in cold-hardiness out of winter is primarily a function of a change in photoperiod (short to long days) and entrance into diapause termination phase rather than the end of diapause. One class of proteins associated with cold-hardiness and diapause are heat shock proteins (HSPs), which act as molecular chaperones (Rinehart et al., 2007; Koštál and Tollarová-Borovanská, 2009; Ragland et al., 2010, 2011). Proteins within the HSP70 family are essential to recovery after chill shock treatments with the silencing of these proteins increasing mortality in flesh flies and linden bugs (Rinehart et al., 2007; Koštál and Tollarová-Borovanská, 2009). Our results show significant downregulation of genes in this family within the E-strain, both in the diapause termination phase and at the end of diapause [e.g. Hsp70AB, Hsp70BC (FDR<0.01 on day 1 and 7)]. Although not significant, a similar pattern was seen for Hsp70BC within the Z-strain, even though they remained in the diapause termination phase for much longer (Fig. 5E,F). Consistent downregulation of HSPs in both strains prior to the end of diapause implies that cold-hardiness is a property of diapause maintenance in moths, rather than the entire diapause developmental pathway.

Once hormones are released from neurosecretory cells in the brain, downstream development and metamorphosis become inevitable. One important hormone, 20-hydroxyecdysone (20HE), plays a role in inducing insect molting and metamorphosis at appropriate developmental stages (Denlinger, 2002) and commonly breaks diapause in many insect species when injected, fed or topically applied (Žd'árek and Denlinger, 1975; Gadenne et al., 1990; Kidokoro et al., 2006; Yamamoto et al., 2008). Therefore, extension of the diapause termination phase should be correlated with suppressed transcription of ecdysone-related genes through decreased availability of 20HE or the inability to bind it. Indeed, Neverland (Nvd), Spook (Spo) and Phantom (Phm) are involved in 20HE production and were upregulated in the E-strain compared with the Z-strain [e.g. Nvd (FDR<0.01 day 1 and 7); Spo, Phm (FDR<0.0001 on day 1 and 7)] (Fig. S1C,D; Table S3). Loss of function of Nvd is known to cause developmental arrest and reduced growth in D. melanogaster through reduced ecdysteroid production (Yoshiyama et al., 2006). In the E-strain, 20HE receptors were also upregulated [e.g. EcR (FDR<0.01 on day 1); Ftz-f1 (FDR<0.01 on day 7)]. Surprisingly, one upregulated receptor gene, Ftz transcription factor 1 (Ftz-f1), controls stage-specific responses to 20HE (Broadus et al., 1999) and was located within the QTL for Pdd (Table 2).

It has long been hypothesized that the circadian clock is intimately connected to photoperiodic responses such as diapause (Bünning, 1936). Supporting this idea, many latitudinal clines show co-variation between clock genes and diapause response (Mathias et al., 2005; Schmidt et al., 2008; Tauber et al., 2007; Cogni et al., 2013; Levy et al., 2015). In addition to the possible involvement of Per, we found the circadian gene PAR-domain protein 1 (Pdp1) to be significantly differently expressed between the strains (e.g. FDR<0.01 on day 1) (Fig. 5A,B; Table S3). However, although predicted to be Z-linked, Pdp1 is not located in the QTL (21.9 Mb in B. mori). In Drosophila, Pdp1 encodes a transcription factor that promotes the transcription of the gene Clock (Clk) (Cyran et al., 2003). In the apple maggot fly, Pdp1 is significantly downregulated from diapause maintenance to the diapause termination phase (Ragland et al., 2011). Similarly, in the E-strain Pdp1 is downregulated from diapause maintenance to day 1 (Fig. 5A). Although the functional significance of this expression change is unclear, conserved expression of this gene between species may be indicative of an important role for Pdp1 in the transition to the end of diapause.

Although diapause is deployed at various life stages in insect taxa, studies in flies and moths now indicate deep conservation of transcriptional pathways during the phases of diapause (Fig. 3; Table S3) (Emerson et al., 2010; Ragland et al., 2010, 2011). However, there appears to be little evidence in these organisms for shared genetic control of natural variation in diapause timing (Tables 2 and 3) (Emerson et al., 2010). The situation might be different within the Lepidoptera, such as Ostrinia moths and swallowtail butterflies, in which variation in diapause termination timing is associated with a common region of the sex chromosome (Hagen and Scriber, 1989; Ording et al., 2010; Dopman et al., 2005; Wadsworth et al., 2015). Hence, independent physiological and genetic mechanisms appear to govern the evolution of diapause timing across distantly related taxa, whereas a common origin may characterize more closely related lineages.

Elucidating the regulatory mechanisms of diapause evolution has long been challenging because of a lack of anatomical landmarks and a complex genetic basis. Our results demonstrate the combined power of high-throughput transcriptome and physiological screens for understanding the evolution of diapause. We anticipate that expanding this framework to other diapause phases (i.e. diapause induction and maintenance), taxonomic groups and at various stages of evolutionary divergence will provide a foundation to understand how animals adapt to climate change and synchronize major life-history events with seasonally varying environments.

We would like to thank the DopChewLew laboratory group for early discussion of this project and Genevieve Kozak for providing feedback on early drafts of this manuscript.