Developmental trajectories of gene expression reveal candidates for diapause termination: a key life-history transition in the apple maggot fly Rhagoletis pomonella

Dormancy is a pervasive adaptive strategy in seasonal environments, suppressing development and metabolism during stressful or resource-limited portions of the year. Here we define dormancy as a state of developmental arrest that organisms enter in response to environmental cues before the onset of stress or resource limitation. The pre-programmed nature of dormancy responses sets them apart from environmentally induced quiescence, i.e. a period when development is arrested in direct response to poor environmental conditions. Programmed dormancy is an alternative and physiologically dynamic developmental pathway co-ordinated by a suite of neuroendocrine and metabolic controls (Denlinger, 2002; Hahn and Denlinger, 2011; MacRae, 2010).

Mechanisms that govern transitions into and out of dormancy are particularly crucial because these transitions synchronize dormant stages with harsh conditions and active life stages with permissive conditions. The importance of synchrony between seasonal events and animal and plant life histories is well illustrated by several recent examples of rapid, contemporary evolution of seasonal timing in migration, breeding and dormancy in response to short-term environmental fluctuations (Hairston and Walton, 1986), species introductions/invasions (Huang et al., 2010; Schmidt et al., 2005; Thomas et al., 2003) and earlier spring emergence timing associated with climate change (reviewed in Bradshaw and Holzapfel, 2006). Many of these studies focus on arthropod dormancy, or diapause, which commonly synchronizes life cycles of temperate insects. Forcing seasonal asynchrony between insect diapause and natural seasonality causes substantial fitness costs (Bradshaw et al., 2004), and similar selective pressures have driven the evolution of latitudinal clines in seasonal diapause timing (Bradshaw and Lounibos, 1977) and the recent evolution of diapause timing concordant with climate change (Bradshaw and Holzapfel, 2001). Evolutionary change in seasonal timing of diapause requires naturally segregating variation in the underlying regulatory mechanisms, and characterizing these physiological processes remains a central challenge.

The mechanisms controlling major life-history transitions such as diapause initiation and termination can act at several different time points along an ontogenetic trajectory. Defining these crucial time points is a necessary first step towards identifying and testing candidate genes underlying the evolution of diapause regulation. Insect diapause is a physiologically dynamic developmental process that can be conceptually divided into three major time periods, or phases, when diapause genes may act: (1) diapause induction and preparation; (2) diapause maintenance; and (3) diapause termination (Kostal, 2006). In species with a facultative diapause response, token environmental cues such as temperature, photoperiod or rainfall stimulate a switch to diapause preparatory development long before an organism actually enters the metabolically and developmentally suppressed diapause state. For example, cues programming diapause are often received in one stage of the life cycle while diapause itself occurs in a later stage (Tauber et al., 1986). In contrast, many univoltine insect species enter a developmentally programmed obligate diapause period each generation and do not require token stimuli for diapause induction and preparation. Once the animal has committed to diapause, the diapause maintenance stage is characterized by developmental arrest coupled with metabolic depression and increased stress resistance. Although animals are in developmental arrest, the diapause maintenance phase itself is physiologically dynamic and changes over time in response to environmental stimuli and internal cues (Denlinger, 2002; Hayward et al., 2005). Like initiation, termination of diapause can be pre-programmed or stimulated directly by an environmental cue, typically either temperature (e.g. increasing temperature after a required chilling period) or photoperiod. Termination may span days to weeks and culminates in either the resumption of active development (in permissive conditions) or the achieved potential to resume development when permissive conditions return (Hayward et al., 2005; Kostal, 2006). Diapause timing is determined by both the timing of diapause induction and the transition from diapause maintenance to diapause termination.

Our understanding of the molecular mechanisms of insect diapause has expanded rapidly over the last decade with the majority of work focusing on comparisons of diapausing and non-diapausing individuals, reflecting mechanisms important to maintaining diapause (reviewed by Denlinger, 2002; MacRae, 2010). Despite the importance to diapause timing, few studies have addressed the molecular basis of diapause induction or the diapause preparatory period (but see Kostal et al., 2009; Urbanski et al., 2010). The molecular mechanisms of diapause termination have received somewhat more scrutiny, often using chemical or environmental manipulations to immediately terminate diapause and initiate post-diapause development. Although specific mechanisms may vary among species, the neuroendocrine system is a crucial regulator of diapause termination across insects. Manipulations with hormones, hormone analogues or other pharmacological agents triggering endocrine activity provide insights into early events during diapause termination. For example, pupal diapause in the flesh fly, Sarcophaga crassipalpis, may be terminated by treatment with the steroid hormone ecdysone or pharmacologically with an array of organic solvents (e.g. hexane) that trigger ecdysteroid production and initiate adult morphogenesis (Denlinger, 2002). Numerous studies have used ecdysteroids or hexane to rapidly terminate diapause in flesh flies and have identified rapid, post-treatment changes in the expression of stress response proteins, metabolite pools, metabolic enzymes, proteins involved in cell cycling and several regulatory pathways including ERK/MAPK and insulin signaling (Fujiwara and Denlinger, 2007; Michaud and Denlinger, 2007; Ragland et al., 2010; Rinehart and Denlinger, 2000; Tammariello and Denlinger, 1998). This is an excellent approach for characterizing downstream gene action, but it cannot reveal processes that occur upstream of the ecdysteroid pulse that signals diapause termination. Further, pharmacological effects of treatments with exogenous hormones or other agents may obscure the pattern of molecular events normally occurring at diapause termination in non-manipulated individuals (Ragland et al., 2009). Therefore, identifying candidate genes that may shape the responses of diapause timing to natural selection will require an understanding of the molecular events that regulate diapause termination in non-manipulated individuals (e.g. Emerson et al., 2010).

Here we describe patterns of transcript abundance along an ontogenetic trajectory from diapause through to diapause termination and subsequent development in the apple maggot fly Rhagoletis pomonella (Walsh): a model species for the evolution of diapause timing, synchronization with host plants and ecological speciation. We take advantage of metabolic signatures of diapause measured using respirometry to sample the transcriptome along a natural trajectory of pupal diapause termination and post-diapause adult development. Our immediate experimental goals were to: (1) pinpoint the precise timing of rapid transcriptional changes that have been implicitly associated with diapause termination and the resumption of development (Denlinger, 2002); and (2) identify important signaling cascades associated with these rapid changes.

Rhagoletis pomonella is an ecological model system for rapid evolution and speciation via changes in diapause timing, providing the opportunity to examine how quickly adaptive shifts in diapause timing can occur. Briefly, this fly historically infested only the fruits of native hawthorns (Crataegus spp.) in North America. However, sometime in the mid-1800s the fly shifted onto a new host plant, the introduced domesticated apple (Malus domestica), and has subsequently formed a genetically distinct host race considered to be an early stage in the speciation process (Feder et al., 1988). Moreover, shifts in diapause regulation and subsequently the timing of adult emergence have created a major pre-zygotic reproductive barrier facilitating speciation (Feder and Filchak, 1999). Adult emergence is synchronized with the fruiting phenology of each host through differences in the timing of diapause termination and adult development. Apples fruit earlier in the summer than hawthorns, and the apple host race has evolved an earlier exit from pupal diapause. Differences in seasonal phenology between the fly host races are clearly observed in the field and in a common environment in the laboratory (Feder and Filchak, 1999; Filchak et al., 2000).

Beyond our immediate experimental goals, functional characterization of gene expression during diapause termination within the ancestral hawthorn host lays the foundation for comparative studies of how diapause mechanisms may be modified during evolutionary differentiation between the host races. Our long-term goal is to identify candidate genes and physiological pathways harboring naturally segregating variation for diapause timing in the field, ultimately providing insights into the evolutionary mechanisms of seasonal timing and speciation in Rhagoletis and other insects.

Rhagoletis (Diptera: Tephritidae) fruit flies infest a variety of North American host fruits from several different families (Berlocher and Bush, 1982). Host fruits are typically available for only a small portion of the year, each fly species has one generation per year, and adults typically rendezvous and mate on host fruits (Feder et al., 1994). Females oviposit into host fruits and larvae consume the fruit, burrow into the soil to pupate, and enter a deep pupal diapause that lasts until the following year. Diapause occurs in phanerocephalic pupae prior to the molt into pharate adults (pupal–adult apolysis) (Denlinger and Zdarek, 1994). Diapausing pupae typically require a chilling period (usually 0–4°C in the laboratory), and after a return to warmer conditions (usually ∼24°C in the laboratory) pupae spontaneously break diapause weeks to months later without any additional temperature or photoperiodic cues (Dambroski and Feder, 2007). Rhagoletis pomonella historically infested a number of native hawthorn species (Crataegus spp.), but a genetically distinct host race has evolved in the last 150 years to infest apples (M. domesticus) recently (∼400 years) introduced from Eurasia (Feder et al., 1988). Apple varieties commonly infested by R. pomonella fruit about three weeks earlier than hawthorns in the late summer or early autumn, and the apple host race has evolved earlier emergence times via earlier timing of diapause termination (Feder et al., 1993).

We sampled a midwestern population of the hawthorn host race to focus on ancestral mechanisms of diapause termination. We collected infested hawthorn fruit from a field site near Fennville, MI, USA, in September 2007. Fruit were transported to the University of Florida, FL, USA, placed in mesh baskets over plastic collecting trays, and held at 24°C, 14 h:10 h light:dark. Pupae were collected daily, held under the same temperature/photoperiodic conditions for 10 days, then transferred to a 4°C incubator simulating winter conditions. This rearing regime promotes field-relevant diapause responses and maximizes survival in the laboratory (Dambroski and Feder, 2007).

We have previously shown that metabolic rate, as estimated by CO₂ production, provides a precise and convenient indicator for the transition from diapause into post-diapause development in R. pomonella (Ragland et al., 2009). Thus, respirometry allows us to sample at clear developmental landmarks across the diapause termination transition and post-diapause development. Diapausing pupae are highly metabolically depressed (∼90% depression compared with the baseline of non-diapause development) and typically remain depressed for weeks to months after transfer to warmer temperatures before the occurrence of a rapid increase in metabolic rate that is clearly measurable over a 24 h period (Ragland et al., 2009). At and beyond this increase, pupae exhibit a characteristic, biphasic metabolic rate trajectory that culminates in adult eclosion (Fig. 1A). The sigmoid, initial rate increase is unique to diapause development and brings metabolic rate back up to the baseline of a U-shaped curve defining a direct (non-diapause) developmental trajectory (Fig. 1B).

A crucial goal of this study was to characterize the transcriptional events occurring between early diapause and late diapause prior to the initial metabolic increase and the 24–48 h period just after the metabolic increase. The initial increase in metabolic rate is consistent with an increase in either cell differentiation and proliferation or metabolic preparation for development. In addition, the molt from the pupal stage into the developing, pharate adult stage happens after the initial metabolic increase. Thus, we previously hypothesized that the initial increase in metabolic rate signals the end of diapause and the resumption of active development (Ragland et al., 2009). We designed a sampling scheme that centered on this transition sampling both before and after the initial increase in metabolic rate to: (1) pinpoint the resumption of development based on changes in the expression of genes implicated in morphogenesis; and (2) determine whether any changes in the expression in signaling pathways associated with cell quiescence and cell growth occur prior to the release from metabolic depression.

Forty diapausing pupae were moved to a 24°C, 14 h:10 h light:dark incubator after 20 weeks of simulated winter at 4°C. We sampled an initial set of pupae after 24 h, and after one week we began daily measurements of metabolic rates on the remaining pupae by estimating CO₂ production using a Sable Systems (Las Vegas, NV, USA) respirometry system and data collection software built around a Li-Cor 7000 CO₂ analyzer (Lincoln, NE, USA) (Ragland et al., 2009). Briefly, we sealed individual pupae in syringe respirometry chambers for 24 h, injected air samples into a flow-through system, and used bolus integration to calculate CO₂ production compared with empty control syringes (Lighton, 2008). We sampled a minimum of six individuals at each of six landmark time points (Fig. 1A) representing (1) early diapause prior to the initial metabolic rate increase (ED), (2) late diapause prior to the initial metabolic rate increase (LD), (3) 24 h (24H) and (4) 48 h (48H) after the initial metabolic increase, (5) the respiration rate plateau (PLAT), and (6) the final, exponential phase of metabolic increase indicating the later stages of adult morphogenesis (EXP). We used the following sampling criteria to define the six landmark time points: (1) 24 h after pupae were transferred from 4°C overwintering to 24°C (ED); (2) two weeks after transfer from 4°C to 24°C (LD; only pupae with stable CO₂ baselines and no metabolic increase over this two-week period were considered LD); (3) pupae whose metabolic rate increased 40% or more from a stable baseline in a 24 h period (24H); (4) 48 h following the initial metabolic rate increase (48H); (5) 5–7 days after metabolic rates stabilized following the initial increase (PLAT); and (6) 2–3 days after pupae had entered the final, exponential phase of metabolic increase (EXP; ∼100% or more increase in respiration rate compared with the plateau). Sampled pupae were replaced with fresh pupae from the 4°C incubator until adequate sampling from each time point had been achieved, and all samples were taken within a 60 day span. Respirometry was performed and samples were taken (as described below) following the same schedule every day (between 12:00 h and 14:00 h), limiting the variation caused by rhythmic patterns of daily expression observed in many genes (McDonald and Rosbach, 2001).

Sampled pupae were immediately homogenized in 1.2 ml Tri Reagent (Ambion, Austin, TX, USA) and frozen at –80°C. Within three months RNA was extracted using the Ambion Ribopure kit following the manufacturer's protocol and assessed for purity and quantity using a NanoDrop Spectrophotometer (NanoDrop Technologies, Inc., Thermo Scientific, Pittsburgh, PA, USA) and the Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA). Labeled amplified RNA (aRNA) from four individuals (out of the six sampled individuals) at each time point was generated using the Amino Allyl MessageAmp II Amplification kit (Ambion). We used 325 ng of total RNA for reverse transcription to synthesize first-strand cDNA, then generated labeled aRNA according to the manufacturer's protocol. Microarray analysis was performed using custom Agilent 4×44k cDNA arrays (described below). A total of 825 ng of labeled aRNA per sample was used for the hybridization, and hybridization and washing were performed according to the manufacturer's protocol using the Gene Expression Hybridization Kit and Gene Expression Wash Buffers (Agilent Technologies, Inc.). Arrays were scanned using a dual-laser DNA microarray scanner (Model G2505C, Agilent Technologies, Inc.). Data were extracted from the scanned images using Feature Extraction 10.1.1.1 software (Agilent Technologies, Inc.). We used a two-color hybridization design with dye swapping, using four hybridizations for each of the six direct comparisons of our six time points (Fig. 1A): ED vs LD; LD vs 24H; LD vs 48H; 24H vs 48H; 48H vs PLAT; and PLAT vs EXP. All data are deposited in the NCBI Gene Expression Omnibus (Accession GSE30785).

We designed 44k arrays from R. pomonella EST sequences generated from a previous 454 pyrosequencing run (Schwarz et al., 2009). We used Agilent eArray to design 60-mer probes for selected ESTs. We initially chose two sense strand probes and one antisense strand probe for each of 11,520 ESTs annotated to the National Center for Biotechnology Information NR database. After correcting for some initial errors in strand orientation and annotation, however, 6082 annotated sequences were each represented by a single sense strand probe. We also included forward and reverse probes from 1886 un-annotated ESTs, 1170 probes designed from diapause candidate genes from S. crassipalpis (Ragland et al., 2010), and 178 probes from Drosophila genes implicated in circadian rhythm, stress resistance and endocrine signaling. The non-R. pomonella probes generally hybridized poorly and are excluded from the results. Antisense probes were initially included on the array to establish a distribution for presence/absence calls (Ragland et al., 2010), but we instead used present/absent calls from the Agilent Feature Extraction software to filter probes. All antisense probes and probes whose green or red signals were not significantly different from background (Agilent Feature Extraction reference) in at least three of the four replicate hybridizations were filtered prior to linear model analysis.

Hybridization intensity data were corrected to Agilent internal controls, background-corrected and Loess normalized (Smyth and Speed, 2003), and analyzed using the Limma package in R (R Development Core Team, 2009), applying a linear model with treatment effects (i.e. sampled time points; four biological replicates per time point) and empirical Bayes estimates of standard errors (Smyth, 2004). We designed contrasts to estimate log₂-fold changes corresponding to each of the direct hybridization comparisons (i.e. ED vs LD, LD vs 24H, LD vs 48H, 24H vs 48H, 48H vs PLAT and PLAT vs EXP), and P-values were corrected for multiple testing using the Benjamini and Hochberg false discovery rate (FDR). Mean log₂-fold changes were calculated for each Rhagoletis gene by averaging over all replicate probes within all probe sets that annotated to the same Drosophila melanogaster gene (several ESTs on the array annotate to the same gene). We estimated a combined P-value for the mean fold change estimates using Fisher's method for combining multiple experimental replicates (Hess and Iyer, 2007).

Enrichment analysis was performed for each time point comparison to examine pathways differentially regulated along the transition from diapause to post-diapause development. To accomplish this we first performed across-array (scale) normalization to the red/green dye intensities. Next we selected probes annotated to D. melanogaster (Blastx search against Flybase protein database with e<10^–5; http://flybase.org/), removed probeset redundancy by selecting the replicate for a given probeset that showed the greatest variance in intensity among arrays for a given EST, and then used the same criterion to select a single probe from sets of ESTs that annotate to the same gene. We used the Gene Set Analysis (GSA) package in R (Efron and Tibshirani, 2007) to test for differential expression of all D. melanogaster-annotated Kyoto Encyclopaedia of Genes and Genomes pathways (KEGG; http://www.genome.jp/kegg/) using 1000 permutations to estimate FDR-corrected P-values. Finally, we employed short time-series analysis to identify groups of genes with similar temporal expression profiles. We used Short Time-Series Expression Miner (STEM) (Ernst and Bar-Joseph, 2006) with expression values log₂ normalized to the first time point. The algorithm generated a set of 50 distinct temporal expression profiles (maximum unit change = two) that were permuted to test for significantly higher numbers of genes assigned than predicted. Profiles were grouped when the pairwise correlation exceeded 0.7. We also tested for over-representation of Gene Ontology (GO) categories within the STEM results, applying a Bonferroni correction to the P-values.

Changes in gene expression largely tracked the metabolic rate landmarks of diapause termination with the most dramatic changes occurring at or after the first increase in metabolic rate (Table 1). Many more genes were significantly differentially expressed (FDR <0.01) over the 24 h and 48 h periods following the initial metabolic increase compared with the 14-day period between early and late diapause. Changes in expression for the 24–48 h period of initial metabolic increase were also more drastic than over the 5–7 days following the initial metabolic increase (48H vs PLAT), suggesting that many genes are transiently expressed as animals exit the metabolically depressed dormant phase. Expression differences were most pronounced late in the trajectory for rapidly developing pharate adults that had entered the final, exponential phase of metabolic increase (Table 1; supplementary material Table S1). This pattern most likely reflects the massive developmental reorganization associated with adult morphogenesis, which is not visible until the animal reaches the respiration rate plateau (Ragland et al., 2009).

Genes that change in expression between early and late diapause, before the clear increase in metabolism that is the first outward sign of diapause termination, are of particular interest because we expect that they represent some of the early steps in the termination process. Thirty-two genes were significantly differentially expressed (at FDR <0.01) between early and late diapause. Although relatively few compared with the hundreds of genes differentially expressed between time points later in the trajectory, many were highly differentially expressed (19 genes >2-fold), and 97 genes were differentially expressed at a less conservative significance threshold (FDR <0.05).

Gene set analysis revealed changes in a series of different KEGG biochemical pathways across the metabolic trajectory of diapause termination and adult morphogenesis. The majority of enriched KEGG pathways were differentially expressed between late diapause and 48 h after the initial metabolic increase (Table 2). The collective sets of co-ordinately regulated genes included pathways involved in: (1) amino acid metabolism (consistent with large-scale protein remodeling); (2) unsaturated fatty acid synthesis (important for cell membrane construction); and (3) target of rapamycin (TOR) signaling, i.e. a pathway that regulates cell growth and proliferation. Many metabolic pathways were also significantly differentially expressed later in the developmental trajectory as would be expected in order to support the large catabolic energetic and anabolic synthesis demands of adult morphogenesis.

Between early and late diapause three KEGG categories including (1) tyrosine metabolism, (2) biosynthesis of steroids, and (3) Wnt signaling were significantly differentially expressed, again supporting differential gene expression prior to the initial metabolic increase. Rhagoletis genes on the arrays in the tyrosine metabolism and biosynthesis of steroids pathways unfortunately are poorly functionally annotated in D. melanogaster. However, both of these pathways do contain genes that are putatively important in terminating pupal diapause and initiating adult morphogenesis. For example, the steroid hormone ecdysone is a crucial signal for terminating pupal diapause and initiating pharate adult development (Denlinger et al., 2005). Expression of steroid synthetic genes suggests priming of the ecdysteroid synthesis pathway prior to the production of the pulse of ecdysteroids that will terminate pupal diapause and initiate adult morphogenesis. Tyrosine and phenylalanine are crucial aromatic amino acids that play important roles in producing the new adult cuticle that will form at the onset of adult morphogenesis. Tyrosine and phenylalanine can be converted by enzymes of the tyrosine metabolic pathway into dopaquinone, and subsequently radicalized quinones that control cuticular hardening in Rhagoletis flies and many other insects (Andersen, 2005). Genes in Wnt signaling have experimentally confirmed roles in insect development, particularly for controlling the cell cycle and transitions from cell cycle arrest to active proliferation (Reya and Clevers, 2005). The upregulation of Wnt pathway genes, therefore, provides further evidence that preparation for the re-initiation of development occurs before the initial release from metabolic depression during diapause.

Time-series analysis, a method that groups genes according to their similarity in expression patterns across several time points, identified nine profiles with significant over-representation of gene members. We graphically depict only six significant profiles with the lowest P-values in Fig. 2, but discuss and provide data for all nine below and in supplementary material Table S2 and supplementary material Fig. S1. Profiles 3 and 4 have a similar shape and form a cluster (correlation >0.7), while profiles 1–2 and 5–6 were distinct (pairwise correlations <0.7) (supplementary material Table S3). All nine profiles were characterized by marked changes in expression coinciding with the first metabolic rate increase (within 24 h or 48 h). The relative magnitude of these changes over 1–2 days rivals the level of changes occurring over longer time periods earlier and later in the trajectory. All of the profiles except for profile 2 exhibited a significant over-representation of multiple GO categories, reflecting the co-ordination of expression of clusters of functionally related genes during diapause termination (Fig. 2; supplementary material Table S2). Although not detectably enriched in any GO categories, the pattern of expression in profile 2 is consistent with our expectations for genes involved in the maintenance phase of diapause; relatively high expression during diapause that decreases substantially as diapause is terminated and adult morphogenesis begins.

We used a physiological marker of developmental status (respiratory gas exchange) to characterize patterns of gene expression across a developmental trajectory from early pupal diapause through to diapause termination and subsequent adult development in the apple maggot fly R. pomonella. In contrast to many manipulative studies, we used animals that were allowed to break diapause in the laboratory without drastic pharmacological, endocrine or environmental manipulation. Furthermore, our examination of global transcriptional patterns of diapause termination is the most detailed time series of an insect diapause transition to date. Previous studies have mainly considered comparisons of diapausing and non-diapausing individuals at a few discontinuous time points (Baker and Russell, 2009; Emerson et al., 2010; Ragland et al., 2010). Thus, our expression data provide a roadmap for understanding the physiological mechanisms and developmental pathways involved in the process of insect diapause termination.

When considering the order of events in diapause termination, several authors have hypothesized that the initial increase in metabolic rates signaling the release from diapause-induced metabolic depression would precede the resumption of development and cell cycling (Ragland et al., 2009; Schneiderman and Williams, 1953). We expected that animals would first ramp up their metabolism for a period of time to prepare themselves for the anabolic and catabolic demands of morphogenesis before re-initiating cellular growth and proliferation. We show here that expression patterns of both metabolic and cell cycling/developmental genes change dramatically during the 24–48 h window in which animals exit metabolic depression characteristic of diapause. Expression continues to change throughout adult development in a pattern consistent with the regulatory events downstream of pupal diapause termination and adult morphogenesis. These data suggest that the release from metabolic suppression coincides with the initiation of post-diapause morphogenesis. Perhaps most interesting are the substantial transcriptional differences observed between diapausing pupae freshly removed from the overwintering treatment and diapausing pupae sampled several weeks later, nearer to completely terminating diapause. These early expressed sets of genes are likely involved in the preparatory programme for diapause termination and represent the earliest identified molecular signatures of diapause termination prior to the release from metabolic suppression. We discuss our results with an emphasis on endocrine and cellular mechanisms that we believe may regulate the transition from metabolic depression and cell cycle arrest in diapause to the resumption of development.

Identifying candidate genes and physiological pathways for diapause termination in the ancestral hawthorn host race provides a foundation for understanding how selection shapes naturally segregating variation in diapause termination timing of R. pomonella. Selection on existing variation likely drove the evolution of earlier termination timing in the derived apple host race. Once we have identified genes that regulate diapause termination, we can begin to explore variants among the host races and their association with diapause timing. More generally, exploring the machinery of diapause regulation will facilitate comparative studies of regulatory elements that are conserved across insects. Insects have diverse diapause strategies that include diapause responses in every conceivable life stage. Currently, large-scale transcriptomic data comparing diapausing animals with non-diapausing animals exists for just a small set of species: the adult diapausing fly D. melanogaster (Baker and Russell, 2009), the pupal diapausing fly S. crassipalpis (Ragland et al., 2010), and the larval diapausing mosquito Wyeomyia smithii (Emerson et al., 2010). An existing comparative analysis suggested a few key physiological elements of diapause that are conserved across broader phylogenetic divides, although overall diapause appears to highly evolutionarily labile on the transcriptional level (Ragland et al., 2010). Throughout the discussion we highlight some of the elements of diapause termination in R. pomonella that show clear similarities with other insect diapause programmes.

The timing of differential expression of thousands of genes supports a 24–48 h period of rapid physiological change that co-occurs with increased metabolism. Results of the time-series analysis show that all common trajectories of gene expression involve a relatively rapid change 24–48 h after metabolic rate begins to increase and that there is also an enrichment of important pathways for aerobic metabolism like mitochondrial electron transport chain activity (supplementary material Table S2). These patterns closely mirror patterns that we have previously described for pupal diapause in the flesh fly S. crassipalpis (Ragland et al., 2010). Our observations are consistent with those in numerous dormant animals wherein gluconeogenesis and other anaerobic pathways are upregulated during diapause (Burnell et al., 2005; Hahn and Denlinger, 2011; MacRae, 2010) and rapidly downregulated following diapause termination (Ragland et al., 2010). For example, in R. pomonella, expression of phosphoenolpyruvate carboxykinase (PEPCK), an enzyme that catalyzes a crucial irreversible step in gluconeogenesis, decreases 2.8-fold (FDR <0.001) during the transition from metabolic suppression in late diapause to 24 h after the initial metabolic rate increase. This result is consistent with previous observations of PEPCK decreasing at diapause termination in S. crassipalpis (Ragland et al., 2010) and also at termination of larval diapause in the mosquito W. smithii (Emerson et al., 2010).

Genes associated with environmental stress responses are also commonly upregulated during diapause to buffer against stressors such as extreme temperatures, aridity and oxidation (MacRae, 2010). Heat shock proteins (Hsps) have been particularly well studied during diapause, and upregulation of select Hsps is universally associated with diapause maintenance (Rinehart et al., 2007). Of the Hsp transcripts represented on our array that demonstrate differential regulation (Hsp70 is represented but not differentially expressed; supplementary material Table S1), all of the small Hsps (<30 kDa) are downregulated by 48 h following the first metabolic rate increase (Fig. 3A; Table 3), consistent with upregulation during diapause maintenance. This pattern is typical of small Hsp expression during pupal diapause in flesh flies (Hayward et al., 2005). However, two large Hsps (Hsp67 and Hsp83) show the opposite pattern and are upregulated by 48 h after the first metabolic rate increase. Hsps can be upregulated during periods of rapid development (Michaud et al., 1997), and Hsp83 (also commonly referred to as Hsp90) in particular shows expression patterns that track pulses of ecdysone during development in D. melanogaster (Thomas and Lengyel, 1986) and appears to be associated with molting in the moth Sesamia nonagrioides (Gkouvitsas et al., 2009). Flesh fly pupae also demonstrate upregulation of Hsp83 after diapause termination, suggesting that expression patterns of some Hsps may be driven primarily by post-diapause development and not by diapause regulation per se (Rinehart and Denlinger, 2000).

Diapause is a state of developmental suppression characterized by cell cycle arrest in either the G0/G1 or G2 phases (Kostal et al., 2009), and the re-initiation of the cell cycle and subsequent cellular proliferation clearly delineates the transition from pupal diapause back into active development and adult morphogenesis (Denlinger, 2002). The expression patterns of genes involved in cell cycling closely track the respirometric trajectory with clear changes in expression that coincide with the initial increase in metabolism. Cyclins and cyclin-dependent kinases (Cdks) are perhaps the best characterized, highly conserved proteins that regulate cell cycling. When bound by their specific cyclins, Cdks promote progression from one phase to the next by a series of protein phosphorylation steps. Several previous studies suggest that expression of the G1/S cyclins CycD and CycE does not differ between diapausing and actively developing stages (Kostal et al., 2009; Tammariello and Denlinger, 1998). CycE exhibits no detectable change in expression across the R. pomonella diapause developmental trajectory, but relative expression of CycD and CycA (both G1/S cyclins) significantly decreases coincident with the first metabolic rate increase (Fig. 3B; Table 3). Cdk4 (bound by CycD) increases marginally but not significantly during the initial 24–48 h of metabolic rate increase (Fig. 3B; Table 3). However, Cdk4 does increase significantly later during the metabolic rate plateau (supplementary material Table S1). In contrast, CycB, a G2/M cyclin, substantially and significantly increases in expression during the first metabolic rate increase as does its associated Cdk, Cdk1 (a.k.a. cdc2 in D. melanogaster; Fig. 3C). Cks30A, a protein that interacts with Cdk1 in cell cycle regulation (Swan et al., 2005), demonstrates an expression profile nearly identical to Cdk1. Levels of each of these G2/M-associated proteins increase and peak during the initial metabolic rate increase, and Cdk4 peaks before the final, exponential increase in metabolic rate. As in most developmental contexts, in R. pomonella tissue proliferation occurs first (thus greater cell cycling activity at an early stage) followed by a period consisting largely of cell growth and differentiation. The group of genes represented in profile 5 from the time-series analysis (Fig. 2) is highly enriched in cell cycle-related GO categories (supplementary material Table S2) and demonstrates exactly the expected pattern with a peak in activity followed by declining activity as cells that have already proliferated begin to differentiate.

Several additional regulators of the cell cycle support the coincidence of developmental changes with the initial metabolic rate increase in R. pomonella. Proliferating cell nuclear antigen (PCNA; a.k.a. mus209 in D. melanogaster), a protein integral to DNA replication, is downregulated during larval diapause in the drosophilid fly Chymomyza costata (Kostal et al., 2009) and upregulated following pupal diapause termination in flesh flies (Tammariello and Denlinger, 1998). We see a similar pattern of PCNA expression in R. pomonella, with strong upregulation during the first metabolic rate increase (Fig. 3D; Table 3). The timing and direction of changes in PCNA expression occurs nearly in lockstep with CycB and Cdk1. There is no evidence to suggest that PCNA directly interacts with either of these proteins, but the similar patterns of expression do suggest that all three genes play a potentially coordinated role in cell cycle progression from diapause termination to post-diapause development. Cdc25 (a.k.a. string in D. melanogaster), a gene that initiates mitosis by activating Cdk1 (Edgar et al., 1994), is absent on our array. But, two genes that interact with Cdc25 are present. Tribbles, a negative regulator of cell cycling (degrades Cdc25) (Mata et al., 2000), is upregulated more than 2-fold in late diapause compared with 48 h after the initial metabolic increase (Fig. 3D). Polo, a positive regulator of a Cdc25 homolog in Xenopus (Kumagai and Dunphy, 1996), is upregulated nearly 2-fold over the same time interval. Like CycB and Cdk1, Tribbles and Polo regulate the G2/M transition. Expression profiles of these secondary regulators in combination with profiles for cyclins and Cdks present a clear pattern of either a rapid increase in mitosis or rapid preparation for mitosis concomitant with elevated metabolic rate. The action of cyclins can be context-dependent (Kolonin and Finley, 2000), and our data point to an important role for G2/M cyclins and associated proteins in the transition from diapause to active development.

Endocrine control of diapause has been studied extensively in insects and diapause typically ends with a hormonal event, often changes in the titres of ecdysteroids and juvenile hormone (Denlinger et al., 2005). Pupal diapause is typically terminated by a pulse of ecdysteroids that triggers adult morphogenesis, and ecdysteroid sensitivity has been shown to be greater closer to termination than early in diapause (Denlinger et al., 2005). Gene set analysis shows enrichment of genes in the steroid biosynthesis KEGG pathway between early and late diapausing pupae, suggesting that diapausing pupae may be generating the machinery to shunt sterols into ecdysteroidogenesis even before the initial metabolic rate increase that signals diapause termination. Several genes in the ecdysone biosynthesis pathway (i.e. ‘Halloween’ genes) (Huang et al., 2008) show marked increases between late diapausing individuals and pupae 24–48 h after the initial metabolic increase (Fig. 3E; Table 3). Expression of prothoraciotropic hormone (PTTH), a crucial brain peptide hormone that induces ecdysteroid production (Rybczynski, 2005), also increases rapidly during the first metabolic rate increase (Fig. 3E; Table 3). Similar diapause PTTH patterns have been reported for the moths Heliothis virescens and Helicoverpa armigera (Wei et al., 2005; Xu and Denlinger, 2003). The combined observations of hormone biosynthesis genes and PTTH are consistent with R. pomonella preparing for ecdysteroidogenesis and then producing an ecdysone pulse that terminates diapause during the initial metabolic increase.

In addition to the control of ecdysteroid production, receptivity to ecdysteroids may be altered in diapause. For example, early diapausing pupae of the flesh fly S. crassipalpis require greater doses of exogenous ecdysteroids to terminate diapause than do pupae early in diapause (Denlinger et al., 2005). Changes in ecdysteroid sensitivity during diapause in S. crassipalpis pupae are correlated with the expression patterns of the two functional components of the ecdysteroid receptor, EcR and ultraspiracle, Usp. EcR is expressed continuously throughout S. crassipalpis diapause, but Usp transcripts are undetectable early in diapause and become abundant later in diapause near termination (Rinehart et al., 2001). Transcript abundance for both EcR and Usp then increase at diapause termination. We observed no change in Usp transcript abundance from early to late diapause in R. pomonella. In contrast to the observations in S. crassipalpis, transcript abundance drops noticeably during the initial metabolic rate increase then remains stably at a lower level throughout adult development (Fig. 3E; EcR is absent on the array). Further validation at the protein level is clearly needed in R. pomonella and other diapausing insects to determine if the observed transcriptional differences among species reflect actual differences in the presence and abundance of functional ecdysteroid receptors at diapause termination.

Expression of primary response genes downstream of ecdysone also supports an ecdysone pulse at or near the initial metabolic rate increase. Once activated by 20-hydroxyecdysone, the ecdysone receptor initiates the transcription of a suite of early transcription factors promoting adult morphogenesis (Fletcher and Thummel, 1995). Both of the two early genes represented on our array (Eip75, Eip78) show a dramatic increase in transcript abundance between the initial metabolic increase and the plateau phase (Fig. 3E; supplementary material Table S1), although their levels do not change immediately during the initial metabolic rate increase (Table 3). Ultimately we will require quantification of ecdysteroid titres along the same developmental trajectory to establish the precise timeline of the hormone pulse and downstream responses. However, our data strongly suggest that an ecdysone pulse does not occur any earlier than the point of release from metabolic depression.

Distinguishing the subset of differentially expressed genes that play a direct regulatory role in diapause termination from the thousands of differentially regulated genes downstream presents a major challenge. But, differentially expressed genes in well-characterized developmental pathways are obvious candidates. For example, the fork-head transcription factor FOXO plays a pivotal role in developmental arrest via the insulin signaling cascade (Puig et al., 2003; Schmidt et al., 2002). Insulin signaling is strongly implicated as a major regulator of dormancy mediated through the pathway's effects on developmental and metabolic suppression (Gerisch et al., 2001; Hahn and Denlinger, 2007; Ragland et al., 2010; Sim and Denlinger, 2008; Williams et al., 2006), and naturally segregating variation in the insulin signaling member PI3K has previously been associated with clinal variation in adult reproductive diapause in D. melanogaster (Williams et al., 2006). FOXO is significantly downregulated by 24 h into the initial metabolic rate increase (1.8-fold, FDR <0.001), but we have minimal evidence to suggest differential regulation of the primary insulin signaling pathway. Three of the major upstream genes are missing on our array (Inr, Pten, Akt), and PI3K (Pi3K92E) is not differentially regulated (supplementary material Table S1). However, insulin signaling regulates FOXO activity via protein phosphorylation (Lee et al., 2001). Thus, upregulation of insulin signaling genes would not provide a direct explanation for downregulation of FOXO transcription under the current models of pathway control.

Elements of the TOR signaling pathway (overlaps and interacts with insulin signaling) do appear to play a role in diapause termination that is consistent with regulation at the transcript level. Several members of the pathway were individually differentially expressed (supplementary material Table S1). Most notably, both Tsc1 and Tsc2 (gigas) were significantly upregulated in late diapause. Both are negative regulators of Tor, which in turn positively regulates cell growth. Outside of the main TOR pathway, recent studies have identified FOXO-mediated Tor inhibition via sestrins. Sestrins inhibit Tor activity through the Amp-activated protein kinase (AMPK), which in turn enhances Tsc1/Tsc2 inhibition of Tor (Lee et al., 2010) (Fig. 4), offering a potential mechanism for integrating both metabolic and growth responses. Our data suggest that this pathway represents a mechanism of Tor inhibition during diapause. The single D. melanogaster sestrin ortholog Sesn (dSesn) is upregulated in diapause, as are FOXO and the γ-regulatory subunit of the multimeric AMPK. In D. melanogaster, the Sesn gene contains eight FOXO-binding regions and is upregulated with the upregulation of FOXO (Lee et al., 2010). Combined with the upregulation of AMPK and Tsc1/Tsc2 in diapause in R. pomonella, our data are consistent with a diapause state that is inhibitory to Tor activity. Sesn is also significantly upregulated in diapause compared with non-diapause pupae in S. crassipalpis (Ragland et al., 2010) (supplementary material Table S1), providing comparative support for a general role for sestrins in diapause regulation. This FOXO cascade inhibits Tor at the protein level. Thus, observing no difference in Tor transcript abundance between late diapause and 48 h after the initial metabolic increase is still consistent with FOXO/Sesn-mediated inhibition of Tor during diapause.

Recently, studies suggest that genes involved in the circadian rhythm pathway may also contribute to diapause regulation. Genetic variation in timeless associates with propensity to enter photoperiodically induced diapause in the drosophilids D. melanogaster (Tauber et al., 2007) and Chymomyza costata (Stehlik et al., 2008), and variation in period similarly associates with propensity to diapause in the flesh fly Sarcophaga bullata (Han and Denlinger, 2009). Several gene expression and knockout studies also implicate circadian rhythm genes as potential diapause regulators, although whether the entire pathway is involved remains unresolved (Schiesari et al., 2010). None of the three main circadian rhythm genes represented on our R. pomonella array [cryptochrome (cry), cycle (cyc) and timeless (tim)] are significantly differentially regulated during the first metabolic rate increase (Fig. 3F; Table 3). However, jetlag (jet), clockwork orange (cwo) and PAR-domain protein 1 (Pdp1) genes, which serve as secondary regulators of the major clock genes, all significantly change expression by 48 h after the first metabolic rate increase; jet is upregulated while cwo and Pdp1 are downregulated. It is difficult to ascribe any functional significance to these patterns because any effects of these genes outside of circadian rhythm control are largely unexplored. In addition, it is possible that changes in circadian rhythms may simply be associated with a switch to active development and play no direct role in diapause regulation in R. pomonella.

The time-series clustering and cell cycle, Tor signaling and endocrine-related profiles support diapause termination and the resumption of development concomitant with the initial metabolic rate increase, but the data also suggest that upstream signals for diapause termination may occur earlier in the trajectory when metabolic rate has not measurably increased. Diapause itself is a molecularly dynamic state reflected in the 32 genes differentially expressed between early and late diapause prior to any measurable increase in metabolic rates (Denlinger, 2002; Kostal, 2006). Moreover, 28 out of the 32 genes differentially expressed between early and late diapause (FDR <0.01) are also significantly differentially expressed 24 h and 48 h into the initial metabolic rate increase (supplementary material Table S1), suggesting that processes occurring at the metabolic transition may be set in motion days to weeks earlier. Several of the differentially expressed early–late diapause genes play important roles in signaling and developmental processes and thus are potential candidates for regulating diapause termination. Bruchpilot (2.8-fold upregulated; FDR <0.001) plays a role in synapse function and its transcription is associated with neuronal differentiation (Wagh et al., 2006). Orct2 (2.3-fold upregulated; FDR <0.001), an organic cation transporter homolog, is transcriptionally upregulated via the insulin signaling pathway during normal tissue growth and proliferation in D. melanogaster (Herranz et al., 2006) and is significantly downregulated during pupal diapause in flesh flies (Ragland et al., 2010). In addition, Orct2 transcription is directly downstream from the Tsc1/Tsc2 complex, and Tsc1 is also differentially expressed between early and late diapause in R. pomonella (1.9-fold upregulated; FDR=0.014). These observations suggest that insulin/TOR signaling may play an early role in diapause termination in addition to the differential regulation of component genes that we observe during the first metabolic increase and apparent release from cell cycle arrest. Finally, Mapk-ak2, a substrate of MAPK, is differentially expressed between early and late diapause. MAPK signaling plays a primary role in cellular growth and proliferation (Seger and Krebs, 1995), and the ERK/MAPK signaling cascade has been closely linked to control of pupal diapause termination in the flesh fly S. crassipalpis (Fujiwara and Denlinger, 2007) and termination of embryonic diapause of the silkworm Bombyx mori (Fujiwara et al., 2006) and the melon beetle Atrachya menetriesi (Kidokoro et al., 2006).

Several KEGG pathways differentially expressed between early and late diapause provide additional evidence of preparatory activity before the resumption of active development and metabolism. The KEGG tyrosine metabolism pathway, crucial for producing the reactive quinones important for cross-linking the adult cuticle at the onset of morphogenesis (Andersen, 2005), was differentially expressed. This is consistent with preparation for new cuticle formation after a molt (from pupa to pharate adult for R. pomonella) that follows diapause termination (Emerson et al., 2010). The KEGG biosynthesis of steroids pathway was also differentially expressed, consistent with early synthesis of precursors for ecdysteroid production. Finally, the developmentally crucial KEGG Wnt signaling pathway was differentially expressed between early and late diapause, providing an attractive candidate for early control of termination. Wnt signaling affects numerous aspects of embryonic and pupal development (Logan and Nusse, 2004). Our data show that genes within the calcium-dependent Wnt signaling sub-pathway display patterns of differential expression consistent with diapause regulatory function. The KEGG Wnt signaling gene set is composed of three related but distinct sub-pathways: the canonical pathway, the planar cell polarity pathway, and the Wnt/Ca²⁺ pathway. Canonical Wnt signaling includes a large network of genes best studied in the context of tissue proliferation and segment polarity during development (Logan and Nusse, 2004). Fold changes of represented canonical Wnt pathway genes were minimal and directionality was inconsistent. But array representation of genes in this pathway was sparse (16 of the 29 genes represented), limiting any inference about canonical Wnt signaling. Patterns of expression in the planar cell polarity pathway were similarly inconsistent and sparsely represented on our array. However, the Wnt/Ca²⁺ pathway has high representation on our array and nearly all elements were upregulated in late diapause compared with early diapause including a 2-fold increase in the final member, the NFAT transcription factor (Fig. 5), which has been linked to stem cell quiescence (Horsley et al., 2008). None of the individual genes in the pathway was significantly differentially expressed after FDR correction (supplementary material Table S1), but gene set analysis is intended to circumvent problems of multiple comparison corrections by identifying collective patterns of functionally related genes that may be only moderately differentially expressed individually (Subramanian et al., 2005). Thus, Wnt signaling, particularly calcium-dependent Wnt, requires further exploration as a candidate pathway for regulating diapause termination.

Our results suggest that release from cell cycle arrest and metabolic depression occur simultaneously in R. pomonella, providing a tight temporal definition of the transition from diapause to active development. This transition returns a diapausing pupa to a physiological state comparable to a direct (non-diapause) developing pupa near the transition to imaginal disc differentiation and pharate adult development. Because diapause represents a graded, dynamic developmental trajectory, a precise definition of diapause termination has been elusive. Kostal defines diapause termination as a phase over which a diapausing individual becomes potentiated for release from metabolic/developmental suppression, culminating in the actual resumption of direct development (the ‘end’ of diapause) (Kostal, 2006). Under this definition the initial metabolic rate increase and concomitant changes in gene expression appear to be the ‘end’ of diapause, while changes during late diapause are consistent with an earlier transition to ‘termination’ that may ultimately determine when the ‘end’ occurs. Higher density temporal sampling of gene expression, endocrine events and cell cycling will be necessary to determine precisely when the termination process begins; in other words, how far in advance of the ‘end’ the machinery for diapause termination is set in motion.

Although this study provides strong evidence for the developmental timing of diapause termination, dynamic gene regulation during diapause requires further exploration and the specific roles of the candidate regulatory pathways require functional confirmation. Here we used whole body homogenates to assess transcript abundance, but early tissue differentiation and endocrine cascades begin in the central nervous system (Denlinger, 2002). Clearly, follow up studies should include brain and imaginal disc tissue-specific analyses. Transcript profiling highlights likely candidates, but ultimately additional observations of protein abundance and activation followed by functional manipulation (e.g. pharmacological or RNAi manipulations) will be needed to confirm the roles of candidate genes and pathways for diapause maintenance and termination. From a bottom-up, genetic perspective, developing population genomic markers and whole genome sequences will facilitate association and quantitative trait loci (QTL) studies of the diapause timing phenotype. Combining functional tools with quantitative genetics will provide a powerful approach for finding the genes and mechanisms that underlie host race formation and rapid life history evolution in the apple maggot fly.

We wish to thank John Fuller, Serra Goudarzi and Genevieve Ochs for assistance with respirometry and Yanping Zhang (UF genomics core) for labeling preparation and microarray hybridization and scanning. Two anonymous reviewers provided helpful comments on the manuscript.