Rapid Recovery Gene Downregulation during Excess-Light Stress and Recovery in Arabidopsis^[OPEN]

The Excess-Light Stress Experimental System

To interrogate transcriptome dynamics during stress perception, recovery, and memory, excess light was applied to plants at ×ばつ growth irradiance using a combination of metal halide and high pressure sodium lamps (1000 μmol photons m⁻² s⁻¹). This light, with a spectral quality similar to that of sunlight, is distinct from cold light sources such as water-filtered light or LEDs. We have previously reported on the relative induction of excess-light-responsive genes under different spectral qualities (Jung et al., 2013). For example, this "warm" light causes H₂O₂ accumulation and an increase in leaf temperature and is required for robust induction of specific excess-light-responsive transcripts such as APX2 (Jung et al., 2013). Here, we further demonstrate that exposing plants grown under standard growth conditions (100 μmol photons m⁻² s⁻¹; see Methods) to 1000 μmol photons m⁻² s⁻¹ irradiance for 60 min caused expected declines in chlorophyll fluorescence parameters (Galvez-Valdivieso et al., 2009; Karpinski et al., 1997; Supplemental Figure 1 and Supplemental Methods). These results confirmed the physiological relevance of this level of light stress, which significantly impaired PSII photochemistry measured at the adaxial side of the leaf, triggered photoprotective mechanisms, and caused a degree of photodamage, impacting leaf photosynthesis (Lichtenthaler et al., 2005). Importantly, it appears to be well within the capacity of a wild type, Col-0, to dissipate much of the excess light energy and repair any damage. At the physiological level, recovery occurs rapidly during the first 60 min following stress, with photodamage continuing to be repaired over 24 h; therefore, for subsequent detailed transcriptomic investigation, we focused on the first 60 min of recovery.

High coverage, strand-specific mRNA-seq [poly(A) selected] was initially performed on plants exposed to 60 min of light stress to identify and characterize the warm excess-light response. In parallel, a subset of plants was subjected to sustained and relatively slow soil-grown drought stress of 9 d without watering, causing a drop in relative water content (RWC) to around 60%. This allowed a comparison between short-term and sustained oxidative stresses. More than one-third of the excess-light-responsive genes at a false discovery rate (FDR) adjusted P value (padj) < 0.05 were similarly responsive to drought (469 commonly upregulated and 703 commonly downregulated transcripts; Figures 1A and 1B; Supplemental Data Sets 1 and 2). Thus, as expected, a short-term warm light stress has a strong overlap with a slow, long-term drought stress. There was a degree of asymmetry in the M verses A plots indicative of a set of lowly expressed genes becoming among the most abundant mRNAs in the transcriptome (indicated by boxes in Figures 1C and 1D).

Given that excess-light stress is often transient, we hypothesized that significantly upregulated transcripts would be highly unstable, which would facilitate a rapid response to alleviation of the stress or changing conditions. To test this hypothesis, for all transcripts elevated by more than 3-fold after 60 min of excess-light stress (padj < 0.05), we extracted half-life values from previous measurements made in Arabidopsis cell culture in the absence of stress (Narsai et al., 2007). Surprisingly, this subset of transcripts had a higher than average half-life at 6.4 h (median 4.5) compared with the transcriptome average of 5.3 h (median 3.8; Figure 1E). Transcripts upregulated by >8-fold had higher average half-lives again at 7.8 h (median 5.8). Drought-induced transcripts displayed the same trend (Figure 1E). The high stability of these transcripts is surprising and suggests that it could take hours for the transcriptome to reset following light stress. For example, HEAT SHOCK PROTEIN101 (HSP101) has a reported half-life of 6.4 h (Narsai et al., 2007). This transcript is upregulated 750-fold during excess-light stress; thus, in the absence of transcription, it could take >61 h to decay back to prestress levels. In light of this result, we directly tested how rapidly the transcriptome recovers.

Design of a Stress and Recovery Time Course for in-Depth Transcriptome Profiling

An excess-light stress and recovery time course, consisting of four phases, was designed to investigate complementary aspects of stress tolerance, recovery, and memory (Figures 2A and 2B). At each time point, the transcriptome was profiled by performing mRNA-seq, sRNA-seq, and parallel analysis of RNA ends (PARE). Whole-genome bisulfite sequencing (WGBS) was also conducted to examine DNA methylation at the three indicated time points.

RNA-seq data were analyzed to determine relative gene expression profiles over the time course, including differentially expressed genes, using edgeR (Robinson and Oshlack, 2010; McCarthy et al., 2012; Supplemental Data Set 3). The expression profiles of replicate time points formed discrete treatment groups upon clustering, reflecting the quality and reproducibility of the RNA-seq data (Supplemental Figure 2A). qPCR was performed, corroborating the normalized expression profiles (Supplemental Figures 2B and 2C).

Individual stress-recovery expression profile plots were created for the 1000 most rapidly downregulated transcripts (discussed below; Supplemental Data Set 4), and a selection are presented in Supplemental Figure 3. These individual time course plots are a resource for unexpected insights into potential novel roles of stress response genes. To give one example, the expression of the transcription factor SALT TOLERANCE ZINC FINGER (STZ ZAT10) peaked during recovery. Previously, ZAT10 has been reported to function specifically during stress exposure (Rossel et al., 2007; Sakamoto et al., 2004); these results suggest that beyond stress tolerance, ZAT10 could play a role in repressing expression of target loci during stress recovery, consistent with earlier reports that it binds a cis-acting element, the A(G/C)T repeat, and represses gene expression (Sakamoto et al., 2004).

Rapid Transcript Downregulation during Recovery

A striking feature of the stress and recovery RNA-seq expression profiles was that rapid induction can be followed by equally rapid decline in transcript abundance. For instance, HSP70 is expressed at relatively low levels of 28 counts per million (CPM) under control conditions (where control conditions refer to measurements made at time point I on untreated plants grown under standard growth conditions, as described in Methods). Upon excess-light stress, HSP70 becomes the 23rd most abundant transcript in the entire transcriptome at 10,746 CPM and then rapidly decreases back to prestress levels during recovery (Supplemental Data Set 4, page 575). We termed this phenomenon rapid recovery gene downregulation (RRGD).

To determine the prevalence of rapid transcript silencing during stress and recovery, the profiles of all transcripts activated by excess-light stress were examined. Many transcripts peak in expression after 30 min of excess-light exposure; thus, focus was turned to the 444 transcripts significantly (padj < 0.05) upregulated by 3-fold at 30 min of excess light (Figure 3A; Supplemental Data Set 5). The 444 excess-light-responsive transcripts were divided into three groups: Category 1 transcripts that maintain high expression into the recovery phase (Figure 3D); Category 2 transcripts that drop in expression specifically during recovery ("recovery-triggered decay"; Figure 3E); and Category 3 transcripts that peak in expression at 30 min and are subsequently downregulated during the remaining light stress period (Figure 3F). For example, HSP101 increased 750-fold by 30 min, with an equally rapid decay during and after light stress (Figure 3C). Interestingly, of the 444 transcripts upregulated by 3-fold or more at 30 min, 44% declined to below basal levels such that they were downregulated relative to untreated by 60-min recovery.

Unexpectedly, during stress recovery, a prominent spike in the relative expression level of many transcripts was observed at the 7.5-min recovery time point (e.g., Figure 3A). Using a high-resolution qPCR time course (Figures 4A to 4C) the expression spike during recovery was replicated a second time but with greater resolution for the three transcripts profiled. We observed that the recovery expression spiked within 3 min and could still be observed at the 6- and 9-min recovery time points. It is possible that this spike could be linked to a trigger for stress recovery, but this was not investigated further. Consistent with the RNA-seq data, the speed of induction for the HSPs observed in the high-resolution qPCR analysis combined with the speed of downregulation during recovery suggest that HSP101 and HSP20-like2 mRNAs are highly unstable (Figures 4B and 4C).

To explore the activation of the stress-recovery response, the genome-wide trajectories of transcripts in plants maintained in excess light for 120 min (IV) were compared with those of plants that were sent to recovery after 60 min (VIII). The plots in Figures 4D and 4E demonstrate that on average for excess-light activated transcripts, if excess-light is maintained, then higher expression is also maintained (red lines); conversely, if the stress is alleviated then expression is repressed (blue lines).

Unstable Transcripts

The rate at which transcripts recover from their peak expression defies the predicted half-lives of many excess-light responsive transcripts, which in cell culture averaged 6.4 h (Figure 1E). For instance, HSP101 transcripts with a reported half-life of 6.4 h in cell culture take just 60 min to recover in planta (Figure 3B). To measure stability during light stress and recovery, we considered two approaches: transcriptional shutoff (Park et al., 2012; Chiba et al., 2013) and in vivo estimation based on abundance measurements (Thomsen et al., 2010). First, we tested experimental procedures to stop transcription during stress. Photosynthetically active protoplasts (Yoo et al., 2007) presented a promising system for application of transcriptional inhibitors; however, we could neither observe reliable induction or recovery decay of any typical excess-light responsive transcripts, ruling out protoplasts as a suitable system for excess-light stress in planta (Supplemental Figure 4). We also performed cordycepin infiltrations into leaves during excess-light exposure with the aim of comparing half-lives before, during, and after stress; however, we could not reliably stop transcription during excess-light stress in intact leaves on the time scale of minutes employed in our time course.

To examine RNA half-lives in our excess-light recovery system, two strategies were undertaken. The initial strategy was to adopt an approach analogous to a prior study in Drosophila melanogaster that estimated RNA decay during development in vivo in the absence of exogenous application of transcriptional inhibitors (Thomsen et al., 2010). We examined the stability and speed of mRNA decay during stress and recovery by estimating maximum possible half-lives (half-life_max). This approach does not factor in ongoing transcription and accordingly will overestimate half-life if transcription is still occurring. For transcripts that exhibit rapid downregulation characteristic of predominant decay, this approach can provide an estimate of the maximum half-life with statistical confidence, but for other transcripts that are not rapidly decayed, this model may be inaccurate. The calculated half-life_max is an upper bound estimate based on the rate of downregulation, with the true half-lives being necessarily equal or lower (more unstable). For each transcript, the net decay was quantified and the in vivo decay constant (k) and half-life_max were estimated by applying an exponential decay model. The net decay values represent log₂ fold-changes between respective time points; hence, a net decay of –1 represents a 50% drop in transcript abundance.

For HSP101 mRNAs, we estimated a half-life_max of ∼12.1 min during recovery using the RNA-seq data. In the high-resolution qPCR analysis, a half-life of 5.2 was estimated. Moreover, given a half-life of 5.2 min, in the absence of transcription, it would take around 50 min for HSP101 transcripts to decay to prestress levels. Supporting our estimation and approach, 60 min after stress HSP101 expression levels returned to within 4- to 5-fold of prestress expression (Figures 3C and 4B).

The half-life_max values for each individual time interval of the time course were collated, revealing that the recovery period was associated with the highest transcript instability (shortest half-life_max). The transcripts with the shortest half-life_max were found to be downregulated most rapidly between 0 to 7.5 min and 7.5 to 15 min recovery, and transcripts downregulated during these periods also had the lowest median half-life_max (Figure 5A). One limitation of estimating half-life_max in vivo is that during short time intervals more stable transcripts do not decay significantly, so their half-life_max cannot be estimated, which could bias the distribution in Figure 5A. Therefore, we took the 1000 most rapid decay events across the whole time course, representing the 901 most unstable transcripts (Supplemental Data Set 6). Of these 1000 fastest decay events, 69% occur during the first 15 min of recovery (Figure 5B) in agreement with the median values (Figure 5A). Nevertheless, directly comparing half-life_max estimated for stress and recovery using the same time interval of 30 min revealed that the initial stress is also associated with very unstable transcripts (Figure 5C). Analysis of the Gene Ontology (GO) categories enriched among the most unstable transcripts revealed greatest enrichment for the Molecular Function transcription factor, which in turn likely accounted for a large diversity of stress responsive Biological Processes Terms that are enriched for this group (Supplemental Data Set 7).

Half-lives estimated here were directly correlated with their respective half-lives in cell culture in the absence of stress (Narsai et al., 2007). There is no correlation (r² < 0.01) between temporal half-life measured during stress and recovery and the half-life determined in the absence of stress (Figure 5D). For example, the measured half-life of APX2 is >17 h in the absence of stress, which is in the 99th percentile of stability (<1% of transcripts are more stable), whereas during the stress treatment between 7.5 and 15 min APX2 has an estimated half-life_max of 12.1 min (only 21 transcripts are less stable, 1st percentile). This result is consistent with the hypothesis that RNA stability is dynamically regulated during stress and recovery.

This initial half-life_max estimation method incorporates measurements made from just two consecutive time points. Given that it is not known if the stability of the mRNAs for a particular gene changes during excess-light stress and/or recovery, this approach was employed to capture any rapid decay event during the time course. However, an additional more robust calculation was undertaken to estimate half-life_max incorporating the additional available time points. Here, a least-squares linear model was applied to the time course expression values to estimate k_max, 95% confidence intervals, se of the decay rate, and recovery half-life_max. This calculation was applied to the 444 genes upregulated by 3-fold during excess-light stress (Figure 3A) and is reported in Supplemental Data Set 5. Category 3 genes were modeled from 30 min of stress until 60 min recovery to capture the period of rapid downregulation following the peak in expression. Given the expression spike between 0 and 7.5 min recovery, Category 2 (and Category 1) genes were modeled from 7.5 min recovery to similarly capture the period of most rapid downregulation. Using this combination of time periods also fitted the highest number of genes with statistical significance (FDR < 0.05). In total, 87% (388/444) of these stress upregulated genes had a recovery half-life_max of <60 min (regression padj < 0.05). This group of 388 highly unstable, rapidly recovering transcripts was defined as "RRGD loci" (Supplemental Data Set 5).

An important point of difference between this investigation and that conducted by Thomsen et al. (2010) is that the previous study analyzed Drosophila during early development, a time period reportedly devoid of transcription. Transcription undoubtedly still occurs during excess-light stress and recovery at some loci. Therefore, we undertook a second modeling approach to better estimate half-lives by accounting for ongoing transcription and decay. This approach also enabled us to test the suitability of the half-life_max estimation in our system. We used nonlinear least squares regression to fit Equation 2 derived from Equation 1 (see Methods) according to Pérez-Ortín et al. (2007), yielding estimates of transcription rates during recovery, and half-lives (termed half-life_NL to denote the model used in the calculation). This model exploits the kinetics of expression: Synthesis of an mRNA follows zero-order kinetics, whereas their decay follows first-order. Therefore, the model incorporates a linear term for transcription (zero-order) and an exponential term for decay (first-order).

The nonlinear model was applied to the 444 excess-light responsive transcripts to calculate half-life_NL. Using this more complex approach reduced our statistical power; however, the model converged and we were able to make estimates for 90 and 95% of Category 2 and 3 genes, respectively (Table 1). Example model results are plotted in Figure 6 for AMT1;2 and HSP101 along with the half-life_max model results. Both models exhibit a similar fit to the RNA-seq data with comparable half-lives. Applying the model to both the RNA-seq data and qPCR data also led to similar half-life estimates for HSP101 (Figure 6B and 6C). The linear model gives a larger half-life estimate than the nonlinear, as expected. However, for these examples gene transcription levels approach zero for both models, and the goodness of fit of the half-life_NL and half-life_max models is not significantly different (e.g., for HSP101, Akaike’s information criterion [AIC] 604 and 602, respectively; ANOVA, P = 0.92). The power of the nonlinear model is better demonstrated using simulated data, where the final steady state does not equilibrate to zero (Supplemental Figure 5). In this case, the goodness of fit of the half-life_NL models is significantly better than half-life_max model (AIC 604 and 611, respectively; ANOVA, P = 0.003).

In total, for 222 transcripts we could confidently (FDR < 0.05) estimate half-lives (all results listed in Supplemental Data Set 5). Among the high confidence stability estimates, Category 2 transcripts had an average half-life_NL of 8.8 min and Category 3 20.2 min; for instance, AMT1;2 had an estimated half-life 6.9 min (Figure 6; Supplemental Data Set 5). As expected, the half-life estimates were lower than the prior half-life_max calculations in all cases. Overall, following FDR adjustment, the linear transcription term was significantly greater than zero [FDR(α) < 0.05] for 30 Category 2 transcripts and 39 Category 3 transcripts, indicating that transcription is ongoing to some level at these loci (Supplemental Data Set 5). We also modeled the half-life of HSP101 during the period of stress (20–120 min) using the qPCR data, which gave a value of 25 min for HSP101 during continuous stress; however, the result was not statistically significant (FDR > 0.05).

Half-lives and k_decay provide a comparable measure of the stability of an mRNA; however, because RNA decay is a first order process, these parameters do not directly reflect actual rates of decay. The decay rate of an mRNA (molecules h⁻¹) is a product of the decay constant (k_decay) multiplied by concentration (C). Accordingly, it is informative to put into perspective the workload of the cellular RNA decay machinery by considering decay rates. Lee and Lee (2003) estimated that Arabidopsis has 100,000 transcripts per cell for average sized transcripts (citing estimates of 100,000 and 500,000 for higher plants; Kamalay and Goldberg, 1980; Kiper et al., 1979); accordingly, one transcript copy would approximate 10 CPM on average for RNA-seq data (without correcting for transcript length bias). This is consistent with estimates in animals: Mortazavi et al. (2008) proposed that one transcript copy approximated 0.5 to 5 FPKM and Marinov et al. (2014) also estimated ∼10 FPKM on average, noting that RNA quantities can vary by an order of magnitude between cells. Working with the estimate that one transcript per cell approximates 10 CPM in RNA-seq data from Arabidopsis leaves, we can thus calculate and compare molecular decay rates. HSP101 had the fastest decay rate of 585 molecules min⁻¹ cell⁻¹ at the 30-min time point (II) during excess-light stress, while RUBISCO ACTIVASE (RCA) was the second fastest decaying transcript. APX2 exhibits a rapid decrease between 7.5 and 15 min recovery, going from 148 to 40 CPM, corresponding to a half-life of 4.0 min. This translates to an estimated average decay rate of 1.4 molecules of APX2 min⁻¹ cell⁻¹. By contrast, during the same period an approximate 25% drop in RCA abundance was detected, which was statistically significant (padj < 0.05) and corresponded to a decay from 14,283 to 10,879 CPM over 7.5 min. At a minimum, this would require turning over in the order of 45.4 molecules of RCA min⁻¹ cell⁻¹ on average at an initial rate of 585 molecules min⁻¹ cell⁻¹ at the beginning of the decay period (assuming exponential decay dC/dt = −k_decay * C). Thus, small relative changes for highly expressed genes may require considerably more enzymatic activity than lower abundance genes.

Profiling RNA Degradome Dynamics during Stress and Recovery

Considering the high instability and rates of decay triggered by stress and/or recovery, one attractive mechanism for RRGD is cis-acting siRNAs, which could be generated from aberrant transcripts produced during the prolific transcription at stress-responsive loci. Therefore, we undertook PARE to profile the RNA degradome, capturing uncapped precursor siRNA template molecules. This procedure captures all uncapped and polyadenylated mRNA molecules in the transcriptome using a custom variation on the Illumina sRNA and mRNA library preparation protocols, providing a snapshot of the degradome (Zhai et al., 2014). Using these data, we can calculate the degradome abundance associated with the mRNA for a particular gene. First, we compared the PARE profile of the untreated sample to previously reported degradome sequencing. A significant correlation between PARE abundance and mRNA abundance (r² =0.64; Figure 7A; additional replicates provided in Supplemental Figure 6) in the untreated samples (I) was observed, consistent with previous reports that also examined untreated samples (Jiao et al., 2008; Addo-Quaye et al., 2008; Willmann et al., 2014). Of particular interest were the 1950 loci that displayed a greater than 2-fold enrichment for PARE tags relative to mRNA abundance (Figure 7A, upper left). As expected, PARE-enriched loci include known microRNA (miRNA) targets and represent canonical miRNA cleavage signatures in the PARE data (Figure 7A, red dots). Some PARE-enriched loci correspond to structural RNAs such as transfer RNAs, small nuclear RNAs, and small nucleolar RNAs as previously observed (Hou et al., 2014), and many also correspond to protein-coding genes. These observations raise the possibility that the abundance and pattern of PARE reads mapping to protein-coding transcripts could be of functional relevance.

The relationship between PARE abundance and mRNA stability was first examined using existing mRNA half-life data (Narsai et al., 2007). We observed no correlation (r² < 0.01) between the stability of a transcript (as measured in a transcriptional inhibition experiment) and the degradome abundances measured in our experiments (Figure 7B). Similarly, there was no correlation (r² < 0.01) between the half-lives measured in our experiments and their PARE abundance at 60 min recovery, a time point when RRGD loci approach equilibrium (Figure 7C).

During the periods where expression levels change, the time-course PARE analysis herein revealed that rapid transcript upregulation is associated with a proliferation in uncapped mRNA molecules (Figures 7D to 7I). One possibility is that unstable transcripts undergoing higher rates of decay could exhibit an increased PARE:RNA-seq ratio. Therefore, PARE and mRNA-seq abundances were profiled for RRGD loci that display rapid induction followed by rapid decay during excess-light stress, and a selection is shown in Figures 7D to 7I. Overall, the degradome abundance largely mirrors mRNA abundance. There is some evidence in these plots that periods of decay are associated with an increased PARE ratio, with PARE abundance appearing to "lag" behind mRNA abundance. PARE abundance spikes to relative levels higher even than mRNA abundance at the 7.5-min recovery time point (V), which is a period of very rapid downregulation (Figures 4A to 4C).

Recently, Hou et al. (2016) and Yu et al., (2016) demonstrated that PARE data from plants can be used to identify loci undergoing cotranslational decay by virtue of ribosome footprints detectable in the degradome. To validate the quality of our degradome data and investigate whether RRGD loci undergo cotranslational decay, we developed a software tool called PAREphase (see Methods). At a genome-wide level, consistent with the previous findings of Hou et al. (2016) and Yu et al. (2016), we also observe 3-nucleotide phasing over coding sequences. Our software tool also replicates the findings of both previous investigations using the published data sets (Supplemental Figures 7C to 7F). For all three data sets, frame 1 is enriched over the other frames (where 0 is the translation frame), with remaining reads distributed more in frame 2 than frame 0 in the direction of the translating ribosome. With respect to ribosome stalling, in agreement with Yu et al. (2016), we find that the peak accumulation of reads occurs at 16 nucleotides and a secondary peak at 17 nucleotides upstream of the mRNA stop codon, corresponding precisely to the 5ʹ boundary of ribosomes stalled with their A site at the stop codon (Figure 8A; Supplemental Figures 7A and 7B).

RRGD loci exhibit clear 3-nucleotide phasing and the 16- to 17-nucleotide peak, demonstrating that these loci undergo cotranslational decay under control conditions (Figure 8B). The ratio of the abundance of the protected frame 1 to the average abundance in frames 0 and 2, is defined as the cotranslational decay index (CDI) and can be used as an indication of relative levels of cotranslation decay (Yu et al., 2016). Here, we further divide the CDI into a phasing decay index for the translating ribosome and a stalled ribosome decay index because the ribosome footprints of translating and stalled ribosome are enriched in different frames (Figure 8C; see Methods). During the recovery period, the RRGD loci exhibit little change in the phased CDI; however, a pronounced increase in their stalled CDI occurs (Figure 8D). This is suggestive of an increase in cotranslation decay during stress and recovery for these loci compared with the genome-wide trend. While these PARE data are only a single replicate per time point, the trend is consistent over the recovery time course.

Stress Recovery and XRNs

PARE specifically captures uncapped, polyadenylated RNA that is in the process of 5ʹ end decay. The increased abundance of PARE tags for excess-light-activated transcripts also suggested the involvement of the 5ʹ-3ʹ decay pathway. To test this, qPCR was undertaken for a selection of RRGD loci using nuclear (xrn2xrn3) and cytosolic (xrn4) XRN mutants. In addition, we generated a triple xrn2 xrn3 xrn4 mutant and included the SAL1 mutant alleles alx8 and fry1-6 that have concurrently impaired function of all the XRNs (Estavillo et al., 2011; Gy et al., 2007). Wild-type plants exhibited strong induction of gene expression during the 60 min of light stress followed by a rapid restoration of mRNA expression levels during recovery. Collectively, stress recovery for the sal1 and xrn mutants proceeded with dynamics similar to those of wild-type plants as opposed to the hypothesis that abundance would decrease only minimally or even continue to increase during recovery (Figure 9; Supplemental Figure 8).

Endogenous Transcripts Protected from siRNA Silencing

The rapid proliferation of uncapped RNAs during stress at excess-light-responsive loci raised the possibility of siRNA biogenesis. Small RNAs were profiled by sRNA-seq in triplicate for each time point of the time course (using the same RNA on which RNA-seq was conducted). De novo annotation and quantification of siRNAs was then performed using ShortStack (Axtell, 2013; Shahid and Axtell, 2014) to determine regulatory siRNA clusters and identify siRNA producing loci. In total, between 16,731 and 27,495 siRNA loci were identified per time point (Supplemental Data Set 8), consistent with previous reports of siRNA cluster number in Arabidopsis aerial tissue (Axtell 2013). On average 76.5% of the siRNA clusters overlap with transposable elements (TEs), accounting for 14,160, or 45.4% of all TEs. By contrast, only 8.6% of the siRNA clusters overlap with protein-coding loci. A total of 1979 genes have siRNA clusters, with 35.5% of these gene-associated siRNA clusters conserved across every sample. Two analyses were then undertaken: an overlap with stress-responsive transcripts and an unbiased global analysis. First, of the 444 significantly light-stress-activated transcripts (Figure 3A), 12 were associated with siRNAs (Supplemental Data Set 9). However, none of these siRNA clusters were found to be stress or recovery responsive: All 12 siRNA clusters were present in untreated, treated, and recovery samples with no clear treatment-dependent expression patterns. Thus, none of the light-stress-induced loci exhibited siRNA production coincident with mRNA or degradome upregulation.

Globally, we found no evidence of stress- or recovery-induced siRNA clusters. Many stress-responsive and rapidly decayed loci do show an increase in small RNA tags; however, these were generally of mixed size classes (18–25) and overwhelmingly aligned to the sense strand. Accordingly, these small RNA reads are most likely degradation intermediates (small fragments of degraded RNA), rather than functional siRNAs. These small RNA decay intermediates are contrasted in Figure 10 with examples of two of the major classes of siRNAs, 24-nucleotide heterochromatic small interfering RNAs and 21-nucleotide phased, secondary, small interfering RNAs (Fei et al., 2013). On occasion, siRNA loci were observed neighboring stress-responsive loci (Figure 10C); however, the siRNA loci were not responsive to excess-light stress and no obvious influence on the neighboring genes was observed.

To complement the global siRNA profiling, qPCR of multiple transcripts in various combinatorial rdr mutant backgrounds was performed to determine whether impairing the siRNA biogenesis pathway led to defects in mRNA decay following excess-light induction (Supplemental Figure 9). REDOX RESPONSIVE TRANSCRIPTION FACTOR1 (RRTF), HEAT SHOCK TRANSCRIPTION FACTOR A7A (HSFA7A), HSP101, and LIGHT-HARVESTING CHLOROPHYLL-PROTEIN COMPLEX II SUBUNIT B1 (LHCB1.4) mRNAs were analyzed, with LHCB1.4 included as a representative excess-light-repressed transcript. For APX2, rdr1 mutants displayed attenuated induction but no defect in APX2 decay during recovery (Supplemental Figure 9A). Backgrounds with combinations of the rdr2 and rdr6 mutations unexpectedly displayed severely delayed induction of APX2 (Supplemental Figure 9B). In fact, other loci (RRTF1, HSP101, and HSF7A) are also characterized by mildly attenuated induction. Nevertheless, with respect to recovery, expression levels decreased in the rdr mutants with high similarity to the profile of wild-type plants.

Excess-Light Stress DNA Methylome

In addition to the key role of RNA decay in the rapid downregulation of transcripts during recovery, transcriptional repression will also contribute to this process. This could directly affect future stress response if repression is maintained or inherited epigenetically. DNA methylation is associated with repressive chromatin states and a diverse range of environmental cues can alter the methylome (Eichten et al., 2014). Thus, whole-genome bisulfite sequencing was performed on four biological replicates for the control (I), excess-light (III), and 24-h recovery (X) samples.

Globally, DNA methylation patterns were highly conserved across all samples (r = 0.96–0.98). However, when correlations between all samples were clustered, two putative preexisting epitypes among our Col-0 plants were observed (Supplemental Figure 10A). That is to say that the samples had one of two distinct epigenome profiles. These epitypes could be distinguished displaying unique DNA methylation states, independent of the experimental treatments received (Supplemental Figure 10B). The seed stock for the plants used in this study was maintained by bulk harvests (rather than an individual line propagated by single seed descent); thus, individuals may be genetically, and epigenetically, separated by multiple generations. We reasoned that any robust and reproducible excess-light response should be evident in a bulk seed stock.

Differentially methylated regions (DMRs) were identified using a Bayesian hierarchical model, incorporating the biological replicates, implemented in the R-based software DSS (Feng et al., 2014). In total, there were 484 DMRs including 389 CG DMRs between the epitypes (Supplemental Data Sets 10 and 11), and hierarchical clustering of methylation levels at these DMR loci revealed two clades distinguished by distinct and consistent DNA methylation differences (Supplemental Figure 10B; representative DMRs are depicted in Supplemental Figures 10C and 10D). Of these, 79/484 (16%) of the epitype DMRs overlap with the 3018 previously reported spontaneous epiallele loci (Becker et al., 2011; Schmitz et al., 2011; Supplemental Data Set 10). Given that these two putative epitypes might be attributed to genetic variation, we identified single nucleotide polymorphisms (SNPs) from the bisulfite sequencing data using BS-SNPer (Gao et al., 2015). Only a few putative SNPs could be identified between samples (Supplemental Data Set 12), ruling out widespread cis-acting genetic variation as an underlying cause of the DMRs and confirming that the plants were all closely related. Similarly, none of the SNPs disrupted the methylation machinery or genes that are known to regulate patterns of DNA methylation; in fact, no SNPs were found in protein-coding genes. Thus, our epitype DMRs likely correspond to potentially spontaneous epiallele formation leading to the emergence of distinct epitypes within our seed stock. Accordingly, these preexisting DMR loci were excluded from the analysis of excess-light effects.

Next, we identified DMRs between treatments (filtering out preexisting epitype loci). Contrasting untreated (I) with 60-min excess-light (III) revealed only six DMRs (FDR adjusted P < 0.05; two example DMRs are depicted in Supplemental Figures 10E and 10F) and contrasting untreated (I) with 24-h recovery from excess-light (X) similarly revealed minimal changes in DNA methylation (Table 2; Supplemental Data Set 13). Thus, under our conditions, a single exposure to excess-light stress has negligible impact on patterns of DNA methylation.

Excess-Light Stress Transcriptional Memory and Recovery-Specific Expression Patterns

Given that RRGD can cause extremely rapid repression, coupled with the observation that many transcripts fall below prestress levels after a 60-min recovery, we examined whether RRGD impaired subsequent stress responses. A generalized linear model (GLM) was implemented in edgeR to identify transcripts that display a significantly different (padj < 0.05) response to excess light upon a second exposure. In total, 547 transcripts responded differently to excess-light applied 24 h after the first exposure, compared with their initial response at the first exposure (based on the GLM test of significance). Of the 444 3-fold upregulated mRNAs (Figure 3A), only two transcripts exhibited a diminished response, indicating that in the vast majority of cases RRGD does not impair subsequent induction. In fact, 72 RRGD transcripts displayed an enhanced response, typical of transcriptional memory. Inspection of individual stress-recovery expression profiles also revealed that transcripts respond differently to subsequent stress exposure, for instance, APX2 and DEHYDRATION-RESPONSIVE ELEMENT BINDING PROTEIN2 (Supplemental Figure 3). Genome wide, the first 60-min excess-light stress (III) led to the differential expression of 1841 transcripts (log₂ > 1, padj < 0.05; 953 up, 888 down). Both repeat stress treatments (IX and XI) elicited a greater global transcriptome perturbation compared with the initial stress, as evident in the heat map (Figure 11A). In particular, reexposure to excess light after just 60-min recovery (IX) caused a significant increase in the number and magnitude of differential expression.

Four dehydration stress memory types were proposed (Ding et al., 2012b, 2013; Avramova, 2015). Building upon these criteria we observe six profiles (Supplemental Figure 11), including hyper, rapid, persistent (Figure 11), gain (e.g., AT1G21940), loss (e.g., AT3g03260), and inversion (e.g., AT1G66700; see Supplemental Data Set 4). Consistent with reports for dehydration stress in Arabidopsis (Ding et al., 2012b, 2013), significant occurrences of hyperresponsiveness were identified; nevertheless, the response to the second stress was predominately the same as the response to the first stress (III versus I compared with XI versus IX). Of the 8800 transcripts responsive to the second stress exposure (XI versus X, padj < 0.05), 394 exhibited hyperinduction and 192 hyperrepression (Figures 11B and 11C; all stress memory transcripts are listed in Supplemental Data Set 14). In total, 189 transcripts displayed persistent-induction memory and 284 persistent repression (Figures 11E and 11F). The expression of many hyperinduction transcripts attained comparable levels in both the first and second exposure; however, expression peaked sooner in the first exposure and subsequently declined to lower levels at the 60-min comparison point (Figure 11D). Moreover, many transcripts in the "persistent" memory class peaked in expression during the recovery period, pointing to a role in stress recovery (Figure 11G). In summary, there is prevalent genome-wide excess-light transcriptional memory in Arabidopsis, independent of DNA methylation.

Beyond rapid resetting, hierarchical clustering of the expression patterns and directionality over the time course revealed a gradient of treatment effects in the heat map (Figure 12A) and also in the multidimensional scaling plot (Supplemental Figure 2A). The expression profiles of samples before stress, during stress, and during recovery were clearly distinguishable. This may point to the activation of recovery-specific expression networks as opposed to merely a relaxation back to the prestressed state. Four general patterns were observed, two are anticipated, namely, induction then recovery and repression then recovery (Figures 12B and 12D). By contrast, the two gene sets in Figures 12C and 12E depict transcripts with stable expression prior to and during excess light (II and III); however, during recovery (V–VIII) expression is specifically upregulated or repressed by the reversion to standard light intensity. Among the group of recovery-induced transcripts, GO term analysis revealed enrichment for transcripts associated with synthesis of ethylene (GO:0009693) and jasmonic acid (GO:0009695; all enriched GO terms are listed in Supplemental Data Set 15).

Rapid Excess-Light Recovery Is Underpinned by Highly Unstable Transcripts

Given that slow relaxation of NPQ and constitutive expression of factors that confer stress tolerance deleteriously impact yield, we propose the response of the transcriptome both during and after stress may be critical in determining the rate of return to optimal growth, while also governing the extent of priming and memory. A key question then is how responsive are mRNAs to "deactivation" or downregulation? An established idiom of transcript kinetics is that RNA stability determines the rapidity of the transition between expression states (Pérez-Ortín et al., 2007). Surprisingly, many excess-light and drought-induced transcripts were reported to be relatively stable in cell culture, with exemplar light stress-responsive transcripts such as APX2 reported to have half-lives of >17 h (Narsai et al., 2007; Figure 1 E). However, as these stability measurements were made in cell culture in the absence of stress (Narsai et al., 2007), it was necessary to measure the transcriptome response during recovery in planta.

In planta, transcripts exhibited exceptional responsiveness. For example, the excess-light treatment upregulated HSP101 and AMT1;2 after 30 min; then, by an hour after stress they recovered to 4.4- and 0.6-fold prestress levels, corresponding to very rapid downregulation (Figures 3B and 3C). Moreover, a high-resolution qPCR, with 3-min intervals, corroborated the speed of this recovery (Figures 4B and 4C). Initially, we estimated mRNA half-lives in planta using the equations that describe mRNA decay in the absence of transcription, as has been applied previously in vivo (Thomsen et al., 2010), providing an estimate of the maximum possible half-life for transcripts for which the model gave a significant fit. Of the 444 transcripts upregulated by excess light, the 388 RRGD transcripts were estimated to have maximum half-lives of <60 min during recovery (Supplemental Data Set 5). This revealed an unexpected lack of correlation (r² < 0.01) between stability in cell culture in the absence of stress versus in planta during recovery, albeit calculated by different methods and experimental conditions (Figure 5D). We also kept in mind that if transcription is still taking place, which is possible for many transcripts, then our estimated stress half-lives represent an underestimation. A contributing factor for this stark difference may be that it is difficult to accurately measure the half-life of lowly expressed transcripts in the absence of stress in a transcriptional shutoff experiment (Narsai et al., 2007), and transcripts could have altered stability in cell culture (Papadakis and Roubelakis-Angelakis, 2002; Xu et al., 2013). Nevertheless, while the most unstable mRNA in cell culture has a half-life of 13.0 min, we identified 71 transcripts with an estimated half-life_NL of lower than 13.0 min and 222 transcripts with a half-life_NL of under 60.0 min. It is clear that the RRGD subset of excess-light-induced transcripts can be highly unstable.

Given the uncertainty surrounding possible ongoing transcription specifically at RRGD loci during excess-light recovery, we performed a more sophisticated least squares nonlinear modeling approach to account for any ongoing transcription. This model is based upon transcription being zero order and decay first order. It also assumes (as do all half-life equations) that a given transcription rate does not gradually change over the modeled period. This means that this approach will not be appropriate in all cases and is recommended only where the rate of change in transcription is very rapid, as described by Pérez-Ortín et al. (2007). For this reason, and due to noise in the data and the limited number of time points, we could only derive half-life estimates that were significant for 222 of the 388 RRGD loci. Statistical power to estimate half-lives could be increased with more sampling points and replicates. Our analyses indeed revealed that a degree of transcription is ongoing for some genes, and once accounted for, estimated half-lives were calculated to be the same or less than the maximum estimate. Indeed, there are instances when the two models fitted the data equally well (Figure 6). However, the power of the nonlinear model was revealed using the simulated data for which it had a significantly better fit (Supplemental Figure 5).

Mechanistic Insights into RNA Decay during Stress Recovery

An intriguing further question concerns the timing of destabilization of RRGD transcripts. For instance, are the transcripts that are generated under conditions of stress unstable from the outset, or can stability change during the course of stress or during recovery? In animals, transcription and decay can also be coupled to enhance the responsiveness of transcripts. Through an mRNA imprinting mechanism, RNA Polymerase II subunits Rpb4 and Rpb7 chaperone stress-induced mRNAs into the cytoplasm and promote degradation (Shalem et al., 2008; Molin et al., 2009; Amorim et al., 2010; Elkon et al., 2010). It is not yet known if similar mechanisms operate in plants, although several studies have previously demonstrated that stress can alter mRNA stability (Petracek et al., 1998; Park et al., 2012; Chiba et al., 2013). Analysis of transcription rates could provide a complementary approach toward determining the relative contributions of transcription and decay (Park et al., 2012). For instance, Pol II chromatin immunoprecipitation and global run-on sequencing methods could be implemented to measure transcription (Hetzel et al., 2016; Erhard et al., 2015). However, these approaches may be problematic on the time scale of minutes due to the time delay required to prepare, fix, or cross-link samples (Madsen et al., 2015). Regardless, our data suggest that posttranscriptional RNA regulation plays an important role in mediating the stress response and rapid recovery, and the answers to these questions will likely require a more mechanistic understanding of the RNA decay pathways involved.

Possibilities for regulating the rapid decrease in mRNA abundance of RRGD transcripts include targeted action by the RNA decay pathways, production of (de)stabilization factors, such as RNA binding proteins (Reichel et al., 2016), and the action of small RNAs. PARE enables a snapshot of RNA decay through profiling of uncapped RNA molecules (Gregory et al., 2008) that are potential templates for siRNA production (Christie et al., 2011; Gazzani et al., 2004; Gregory et al., 2008; Gy et al., 2007). We observed that upregulation of mRNA levels during excess-light stress was coincident with a significant upregulation in PARE abundance, which could suggest activation of a 5ʹ-3ʹ RNA decay pathway with a concomitant increase in the pool of potential RDR substrates for siRNA biogenesis against these transcripts.

Interestingly, multiple lines of evidence indicate that impairing RNA decay (Zhang et al., 2015) or decapping (Martínez de Alba et al., 2015) in a range of mutant backgrounds can lead to siRNA production at endogenous loci and potential silencing of endogenous genes. Several examples of environmentally responsive siRNAs have been reported (Yao et al., 2010; Moldovan et al., 2010). However, the vast majority of studied siRNAs in plants do not respond to environmental perturbations, and environmental conditions have yet to be identified where endogenous genes become subject to cis-PTGS in this manner. Here, we find despite prolific and rapid transcription during excess light, no endogenous loci attracted cis-PTGS. Moreover, we observed upregulation of uncapped RNAs, which are siRNA biogenesis triggers (Christie et al., 2011; Gazzani et al., 2004; Gregory et al., 2008; Gy et al., 2007); yet, it appears that this alone was not sufficient to trigger PTGS. Consistent with these observations, mRNA decay during recovery was not impaired in the rdr mutants analyzed (Supplemental Figure 9). In summary, we found no evidence for a role for RDRs and siRNAs in the recovery following excess-light stress.

Previously, PARE/degradome data sets have primarily been used to investigate the activity of miRNAs (Addo-Quaye et al., 2008; Gregory et al., 2008; German et al., 2008), with only a few reports considering the broader use of PARE data (Jiao et al., 2008; Hou et al., 2014; Zhang et al., 2013). One possibility is that loci undergoing higher rates of decay could exhibit an increased PARE:RNA-seq ratio. However, we could not find convincing evidence for such a relationship. PARE:RNA-seq ratios had low correlation to mRNA half-lives (Figure 7B), and we did not see evidence of a shift in the PARE:RNA-seq ratio during time intervals when transcript abundance declined rapidly (Figures 7D to 7I). Recently, PARE analysis was also used to study ribosome stalling and cotranslational decay in plants (Yu et al., 2016; Hou et al., 2016). Extending on these discoveries, we applied this method to a time course during environmental perturbations for the first time. We discovered that RRGD loci exhibit patterns of PARE read phasing characteristic of targets of cotranslation decay (Figure 8B). Thus, a significant mechanism regulating the abundance and decay of RRGD loci is likely to involve the cotranslation decay pathway. Genome wide, there appears to be a decline in the ribosome stalling index average for all loci during stress; the significance of this is not clear. Interestingly, we observed a significant increase in the relative cotranslational decay index at RRGD loci during the recovery period, compared with all loci (Figure 8D). At this stage, we can only speculate as to the cause, but these data are consistent with the specific targeting of these transcripts by this decay mechanism.

The increased abundance of PARE tags for excess-light-activated transcripts (Figures 7D to 7I), coupled with the targeting by cotranslation decay suggests the involvement of the 5ʹ-3ʹ decay pathway in the regulation of the excess-light response, especially as XRN4 is known to have substrate specificity related to conserved sequence motifs in target mRNAs (Rymarquis et al., 2011). XRN4 and LA RELATED PROTEIN 1A associate during heat stress (15 min at 38°C) to facilitate a massive heat-induced mRNA decay process targeting more than 4500 mRNAs (Merret et al., 2013). Under these conditions, xrn4 mutants are highly susceptible to prolonged moderate heat stress. By contrast, another investigation using a different heat-stress regime of 44°C for 3 to 5 h found that xrn4 mutants are more tolerant to this transitory heat shock (Nguyen et al., 2015). Of note, HSFA2 expression recovers to prestress levels more slowly in xrn4 mutants.

Given these findings, it is surprising that xrn4 mutations do not impede recovery for the transcripts analyzed herein. Similar results were observed for xrn2 xrn3, the xrn2 xrn3 xrn4 triple mutant, and sal1 mutants. There is now convincing evidence that cotranslational decay occurs in Arabidopsis (Merret et al., 2015; Hou et al., 2016; Yu et al., 2016). Together, these prior reports have established that cotranslational decay is mediated by the 5ʹ-3ʹ exonuclease activity of XRN4. Loss of XRN4 has been shown to lead to a loss of cotranslational decay, demonstrated through degradome analysis (Yu et al., 2016). Accordingly, loss of XRN4 also leads to the accumulation of transcripts that normally undergo cotranslation decay (Merret et al., 2015). Considering these prior findings, we hypothesized that the decay of RRGD transcripts during stress and recovery would be impaired in xrn4 mutants. In fact, we hypothesized that RRGD transcripts would accumulate to high levels in xrn4 mutants. By contrast, we observed that lesions in XRN4 had little effect on the rate of rapid downregulation of RRGD loci, which was consistent across multiple xrn4 alleles. On the one hand, this result is in stark contrast to prior observations by Merret et al. (2015) where xrn4 lesions had a significant effect on gene expression. Our hypothesis to reconcile these seemingly contradictory observations is that, in the case of excess-light stress recovery, the cotranslational decay pathway is redundant with other RNA decay pathways, which operate in parallel and/or independently. One possibility is 3ʹ-5ʹ decay by the exosome. To establish the existence of an independent pathway, PARE analysis could be performed in xrn4 mutants to confirm that cotranslation decay has been abolished during recovery. While we have inferred this genetically using xrn4 mutants, it would be valuable to confirm this, for instance, by conducting PARE sequencing and subsequent ribosome footprinting analysis using PAREphase upon xrn4 mutants subject to stress and recovery. If cotranslational decay is abolished and rapid downregulation still occurs, then this would point to an additional mechanism. Elucidation of the parallel mechanisms for RRGD decay could then be investigated through further analysis of stress recovery in mutants impaired in other aspects of RNA decay, for instance, 3ʹ-5ʹ exosome-mediated decay.

An intriguing observation in our stress recovery data is the prominent spike in relative expression levels observed for hundreds of transcripts between 3 to 9 min of recovery (Figures 3A and 4A to 4C). It is attractive to speculate that this could correspond to the release of a repressor or a change in chromatin state, for instance, histone exclusion or variants (To and Kim, 2014; Coleman-Derr and Zilberman, 2012). While we do not yet have an explanation for this expression spike, it does point to the possibility of a recovery-activated mechanism preceding and triggering the rapid decay during recovery.

DNA Methylation

Sixty minutes of excess light and a 24-h recovery did not induce any substantive changes in DNA cytosine methylation. It remains to be determined whether any of the six DMRs associated with excess-light stress are of functional significance, but none are associated with known abiotic stress loci. This result was surprising given that DNA methylation changes can be induced by drought, flooding, nutrient limitation, temperature shock, pathogen infection, high salinity, heavy metal exposure, UV radiation, and herbivory (reviewed in Herman and Sultan, 2011). However, consistent with our results, phosphate starvation-induced changes in DNA methylation observed in rice (Oryza sativa) were not found in Arabidopsis, likely owing to the low prevalence of TEs compared with rice (Secco et al., 2015). Additionally, transcriptional changes preceded changes in DNA methylation, which fully manifest only after 24 d (Secco et al., 2015). Thus, the low frequency of TEs in Arabidopsis compared with other species and the short-term and transient nature of the excess-light stress applied herein might explain the lack of DMRs, and it is possible that a longer term or more repetitive excess-light regime would alter the DNA methylome.

Toward Understanding and Enhancing Stress Recovery

We have demonstrated that an aspect of stress recovery involves the rapid resetting of transcript levels that may be underpinned by high RNA instability and decay. In this regard, the recovery process could be altered by manipulating RNA stability or RNA decay processes. Slowing down transcriptome recovery could facilitate a stress memory, increasing future tolerance, although likely with a yield trade-off. More rapid or timely recovery, if possible, might be advantageous. An intriguing possibility is that "memories" could be encoded through retained changes in RNA instability/responsiveness.

We hypothesize that stress recovery is an active process, rather than explained solely by relaxation of the mechanisms engaged during stress. The rapid recovery response is not simply triggered by the length of the stress or upon reaching a certain expression level, suggesting that alleviation of the light stress could be the trigger for the recovery response (Figures 4D and 4E). Whether recovery cues could trigger destabilization of mRNAs or if this response is wholly dependent on transcriptional repression is a fascinating question. In addition, of particular interest are the atypical expression profiles, for instance, hundreds of transcripts maintain steady expression levels before and during the stress, then become up- or downregulated during the recovery phase (Figures 12C and 12E). The new expression states could be a consequence of the stress, reflecting a new homeostasis; alternatively, this could be evidence of a specific function during recovery. As such, these transcripts and expression networks are a key target for future investigations that, like the RRGD loci, might be manipulated to define the mechanism and determine the potential for optimizing stress recovery.

Plant Growth and Stress Conditions

For all experiments Arabidopsis thaliana Col-0 plants were cultivated in soil under a 12-h photoperiod, 100 (±25) μmol photons m⁻² s⁻¹ PAR, and 23/22°C (±2°C) day/night temperatures (standard growth conditions). For excess-light treatments, 10 ×ばつ growth irradiance PAR was applied, ∼1000 μmol photons m⁻² s⁻¹ using a temperature controlled Conviron chamber (21°C, 55% humidity) fitted with 4 ×ばつ 250 W metal halide lamps and 4 ×ばつ 250 W high-pressure sodium lamps. For recovery phases, plants were returned to their standard growth cabinets at 100 (±25) μmol photons m⁻² s⁻¹. Light stress treatments were always performed between 14:00 and 16:00, which is midday in the day/night cycle (08:00–20:00), using mature 3-week-old Arabidopsis plants. Whole rosettes were excised from their roots and immediately frozen in liquid nitrogen for RNA or DNA extraction. Drought treatment was applied by withholding water from 3-week-old plants for 9 d, leading to wilting and a RWC of around 60% in representative plants. RWC was calculated as Inline graphic. A subgroup of drought-stressed plants were rewatered to confirm viability and to assess the capacity for recovery. Rewatered plants readily recover in 3 to 4 d. Arabidopsis mutants alx8, fry1-6 (SALK_020882), xrn2-1 (SALK_041148), xrn3-3 (SAIL_1172_C07), xrn2-1 xrn3-3 (xrn2 xrn3), xrn4-6 (SALK_ 014209), and ein5-6 were as previously described (Estavillo et al., 2011). Other mutant lines, xrn4-5 (SAIL_681_E01), rdr1-1 (salk_112300), rdr2-1 (SAIL_1277H08), rdr2-2 (salk_059661), and rdr6-15 (SAIL 617_H07), were obtained through TAIR. An xrn2 xrn3 xrn4 mutant line was generated by crossing xrn2 xrn3^+/− plants with xrn4-6. The xrn2 xrn3 and xrn2 xrn3 xrn4 seed stocks were maintained with xrn3 in a heterozygous state and homozygous plants identified each planting by PCR genotyping and morphological phenotyping.

Chlorophyll Fluorescence Measurements

Chlorophyll fluorescence measurements were performed using an Imaging-PAM chlorophyll fluorometer K-4 and ImagingWin software application (Walz). For 60-min excess-light treatments, 1000 PAR actinic light was provided either by a Conviron chamber or by the PAM LEDs. Ground fluorescence (F₀), variable fluorescence (F_v), maximum fluorescence (F_m), as well as light-adapted F, F_m′, and F₀′ were measured then the following parameters were calculated using R according to the following formulae (Baker, 2008): Inline graphic (maximum PSII quantum efficiency); Inline graphic (maximum PSII quantum efficiency, light adapted); Inline graphic (PSII operating efficiency, ∆F/F_m′ or Inline graphic_PSII); Inline graphic (PSII efficiency factor, mathematically identical to qP); Inline graphic (nonphotochemical quenching); Inline graphic(coefficient of nonphotochemical quenching). F_v/F_m measurements were taken following a 20-min dark adaptation. The constituents of qN were estimated as described by Lichtenthaler et al. (2005) based on the dark relaxation of F_m measured at 1, 5, and 20 min following the excess-light stress. These values only reflect the approximate constituents of qN due to variations and overlap in relaxation kinetics of each component, which also vary with environmental factors (Baker, 2008; Maxwell and Johnson, 2000); however, they provide an indication of the level of photodamage and are useful for comparing different excess-light stress regimes.

RNA Isolation and Quality Assessment

For transcriptome experiments, RNA was isolated using Trizol reagent (Life Technologies) using a procedure adapted from Allen et al. (2010). Briefly, up to 100 mg of snap-frozen tissue was ground then lysed in 1 mL of Trizol with gentle agitation. Following 5 min incubation at room temperature, the organic phase was extracted twice with 200 μL of chloroform. The RNA was precipitated by addition of an equal volume of 100% isopropanol and incubated overnight at –20°C. RNA was recovered by centrifugation (20,000 relative centrifugal force) at 4°C for 20 min and washed with 70% ethanol, air-dried at room temperature, resuspended in water or Tris buffer, and stored at –80°C. All RNA-seq experiments, i.e., RNA-seq, sRNA, and PARE, were performed in biological triplicate on paired samples from the same RNA extraction. RNA quality was assessed using a Bioanalyzer (Agilent) or a LabChip GXII (Perkin-Elmer) for RIN > 7.0. The size and concentration of all sequencing libraries were assessed using a MultiNA (Shimatzu) or LabChip GXII capillary electrophoresis instrument and a Qubit Fluorometer (Life Technologies).

RT-qPCR

For RT-qPCR analysis, RNA was either extracted using the Trizol procedure or using the Spectrum Total RNA kit (Sigma-Aldrich). RNA was reverse transcribed into cDNA using SuperScript III (Invitrogen, Life Technologies) and oligo(dT18VN) primers. Gene expression was analyzed by quantitative real-time PCR (qPCR) on a Roche LightCycler480 using Sybr Green I (Roche Diagnostics). Raw fluorescence data were exported and analyzed using LinRegPCR (Ramakers et al., 2003; Ruijter et al., 2009) to perform background subtraction, determine PCR efficiency, and calculate starting concentration (N0; in arbitrary fluorescence units). Samples were then normalized against CYCLOPHILIN5 (CYP5, AT2G29960), GLYCERALDEHYDE-3-PHOSPHATE DEHYDROGENASE C2 (GAPC2, AT1G13440), and PROTEIN PHOSPHATASE 2A SUBUNIT A3 (PP2AA3, AT1G13320) and expressed as fold changes against untreated wild-type plants. At least three biological replicates (individual plants) per genotype per experiment were sampled, and each qPCR reaction was run in duplicate or triplicate. Primer sequences and descriptions are provided in Supplemental Data Set 16.

Sequencing and Analysis of mRNAs

A summary of all sequencing experiments is provided in Supplemental Data Set 17. For the RNA-seq experiment comparing drought and excess-light stress, 3-week-old control and excess-light treated plants were directly compared with the drought stressed plants (4 weeks old following drought stress). The Illumina TruSeq V2 Kit was used for RNA-seq library preparation as per the manufacturer’s instructions (TruSeq User Guide Revision D) with a custom dUTP incorporation step during second strand synthesis to introduce strand specificity (adapted from Parkhomchuk et al., 2009). Four micrograms of total RNA was used for input and fragmented for 8 min. Following first-strand synthesis, cDNA was precipitated to remove dNTPs; per 25 μL first-strand reaction, 5 μL of 5 M NH₄OAc (ammonium discriminates against precipitation of dNTPs), 2 μL of the coprecipitant Glycolblue (Ambion), and 110 μL 100% ethanol were added and samples incubated overnight at –80°C. The following day samples were recovered by centrifugation (20,000 relative centrifugal force) for 20 min at 4°C, and pellets were washed with 70% ethanol, air-dried for 5 min, and then resuspended in custom second strand buffer incorporating a dNTP/dUPT mix. The second strand mix contained 4 μL of ×ばつ RT buffer from the SuperScript III kit (Life Technologies), 2 μL 50 mM MgCl₂, 2 μL 0.1 M DTT, 1 μL 50 ng/μL random hexamers, 49 μL water, 8 μL Second Strand Synthesis buffer (NEBNext; B6117S), 10 μL dNTP/dUTP mix (Fermentas; R0251), and 4 μL Second Strand enzyme mix (NEBNext; E6111S). Samples were incubated for 2.5 h at 16°C. Prior to the PCR enrichment step (after adapter ligation), the USER enzyme was used to degrade the dUTP containing coding strand to give a reverse stranded final library. Per sample, 1 μL of USER enzyme mix was added and samples incubated at 37°C for 15 min. Using pilot PCRs, an optimal PCR cycle number of 15 cycles was determined empirically for amplification of all samples and final amplification was performed with a half volume PCR. Following PCR, clean-up samples were pooled in equal molar ratios and sequenced (100 bp paired end) on one lane of the HiSeq 2000 at the Biomolecular Resource Facility (BRF) at the Australian National University.

All bioinformatic analysis pipelines are freely available online on GitHub (https://github.com/pedrocrisp/NGS-pipelines). Quality control was performed with FastQC v.0.11.2. Adapters were removed using scythe v.0.991 with flags -p 0.01 for the prior and reads were quality trimmed with sickle v.1.33 with flags -q 20 (quality threshold) -l 20 (minimum read length after trimming). The trimmed and quality filtered reads were aligned to the Arabidopsis genome (TAIR10) using the subjunc v.1.4.6 aligner with -u and -H flags to report only reads with a single unambiguous best mapping location (Liao et al., 2013). Reads were then sorted, indexed, and compressed using samtools v1.1-26-g29b0367 (Li et al., 2009a) and strand-specific bigwig files were generated using bedtools genomecov v2.16.1 (Quinlan and Hall, 2010) and the UCSC utility bedGraphToBigWig for viewing in the Integrated Genomics Viewer (IGV; Robinson et al., 2011). A summary of all Illumina sequencing experiments performed, describing sample names, indexing, lane pooling, yield of quality filtered reads, and number of genome mapped reads is provided in Supplemental Data Set 17.

For standard differential gene expression testing, the number of reads mapping per TAIR10 gene loci was summarized using featureCounts v.1.4.6 with flags -p and -c to discard read pairs mapping to different chromosomes and the -s 2 flag for (reverse) strand specificity, multi-mapping reads and multi-overlapping reads (i.e., reads mapping to overlapping regions of more than one gene loci) were not counted (Liao et al., 2014). Reads were summarized to parent TAIR10 gene loci rather than individual splice variants by summarizing to the genomic coordinates defined by the feature "gene" in the TAIR10 GFF3 reference. There are 28,775 gene loci in the current TAIR10 release and 19,444 genes were detected as expressed in this RNA-seq data set. Gene loci are annotated with the description of the representative gene model listed in TAIR. Statistical testing for relative gene expression was performed in R using edgeR v.3.4.2 (McCarthy et al., 2012; Robinson and Oshlack, 2010; Robinson et al., 2010; Robinson and Smyth, 2007, 2008). Reads mapping to rRNA were removed (contamination rate <0.1% for all samples); organelle transcripts were removed and only loci with an abundance of at least 1 CPM in at least three samples (∼10–20 reads for each replicate in one sample group) were retained, yielding data for 17,402 gene loci. GO of the concordantly responsive genes was examined for enriched ontology categories using AmiGO v1.8 (Boyle et al., 2004).

For the stress and recovery time course, RNA-seq libraries were prepared using the Illumina TruSeq V2 Kit (RS-122-2001) as per the manufacturer’s instructions (TruSeq User Guide Revision D). Four micrograms of total RNA was used for input and fragmented for 7 min. Using pilot PCRs, an optimal PCR cycle number of 10 cycles was determined empirically for amplification of all samples and final amplification was performed with a half-volume PCR. Following PCR cleanup, samples were pooled in equal molar ratios and sequenced over two lanes, 18 samples per lane (including three samples from an unrelated project) on the HiSeq 2000 in 100-bp paired-end mode at the BRF. Sample pooling was such that at least one replicate from each time point was represented in each lane; no effect of sequencing lane was observed in sample clustering (Supplemental Figure 2A). RNA-seq analysis was performed as described above for the excess light and drought comparison, except at the featureCounts step, reads were summarized in a non-strand-specific manner in accordance with the library preparation. In total, 17,991 loci passed filtering and were tested for differential expression. Multidimensional scaling plots (PCA) were created using edgeR, using a variation on PCA appropriate for digital gene expression. A set of top genes were chosen that have largest biological variation between the libraries (those with largest tag-wise dispersion treating all libraries as one group). Then the distance between each pair of libraries (columns) is the biological coefficient of variation (square root of the common dispersion) between those two libraries alone, using the top genes. Among the transcripts upregulated by >3-fold (padj < 0.05) after 30 min of excess-light stress (II), Categories 1 to 3 were defined as follows: Category 3 transcripts reduced in abundance by >2-fold between 30 and 60 min stress (II versus III); Category 2 reduced in abundance by <2-fold between 30 and 60 min stress (II versus III) and >1.5-fold between 60 min stress and 30 min recovery (III versus VII); and Category 1 reduced in abundance by <2-fold between 30 and 60 min stress (II versus III) and again <1.5-fold between 60 min stress and 30 min recovery (III versus VII).

Transcriptional Stress Memory Types

To identify transcripts that display a significantly different response to excess light upon a second exposure, first a GLM was implemented in edgeR (padj < 0.05). The prevalence of hyperinduction memory (Supplemental Figure 11) was further analyzed by filtering for transcripts that are differentially expressed after the first exposure (III versus I, padj < 0.05), relaxed to prestress levels during recovery (X versus I, padj > 0.05) then were differentially expressed during the second stress when contrasted with the first stress (III versus XI, padj > 0.05). Transcripts displaying a persistent memory type were defined as those that (1) were differentially expressed at all three time points compared with time 0 and (2) fluctuated by less than ±2-fold between each successive time point.

RNA Half-Life Calculation

Maximum half-lives (half-life_max) were calculated according to the exponential decay model using the equation Inline graphic; that is, half-life_max is equal to the natural log of a half divided by the negative of the decay constant k. The decay constant is proportional to the natural log of the change in mRNA abundance over time, Inline graphic (Ross, 1995; Lam et al., 2001; Gutierrez et al., 2002; Narsai et al., 2007; Thomsen et al., 2010). Using this equation, half-lives_max were calculated based on changes in steady state abundance determined using RNA-seq (CPM values) or qPCR. This procedure assumes that decay during recovery is constant; however, we cannot rule out the possibility that stepwise changes in half-life occur during the recovery phase. For the initial survey of the entire transcriptome, all decay events were quantified using the change in abundance between any two consecutive time points. To identify RRGD loci, the decay constant was determined using a linear regression on the time course III, V-VIII for Category 2 transcripts, and II, III,V-VIII for Category 3 transcripts. Least squares linear regressions were performed in R on log-transformed CPM values to give the coefficient (K_decay), 95% confidence interval, and se of K_decay. An R² value was calculated for the fitted model of each transcript with an accompanying F-statistic and P value, testing the null hypothesis that the fitted regression is equivalent to a y-intercept-only model (a significant P value indicating that the fit of a y-intercept only model is reduced compared with the fitted regression). RRGD loci were defined as those that were upregulated by >3-fold (padj < 0.05) after 30 min of excess-light stress (time point II) and had a recovery half-lives_max of < 60 min (padj <0.05).

To more accurately estimate half-lives in vivo and to account for ongoing transcription, we undertook a further nonlinear modeling approach. The decay constant of a transcript may be estimated from data of transcript abundance over time. Pérez-Ortín et al. (2007) showed that the concentration of a transcript m as it changes from an initial steady state Inline graphic toward equilibrium at a final steady state Inline graphic, can be modeled as:

where k is the decay constant and t is time. Using nonlinear least squares regression, we can estimate the three terms: Inline graphic, Inline graphic, and k according to:

This assumes an instantaneous change in the transcription rate TR at t = 0 and that both the new TR and k remain constant over the modeled period. Since we can solve for k and m, and given that steady state concentration is proportional to transcription rate divided by decay rate Inline graphic (Pérez-Ortín et al., 2007), then:

Similarly, where transcript abundance is at equilibrium (steady state) at the initial time point, one could also solve Inline graphic for Inline graphic (because Inline graphic at steady state). This model (Equation 1) is consistent with the half-life_max model used to estimate decay constants in transcriptional halt experiments, such as in Narsai et al. (2007), because it is mathematically equivalent when the transcription rate is zero such that the final steady state concentration Inline graphic also reaches zero (see Supplemental Methods for derivation). Importantly, this model can estimate the true decay constant when transcription is ongoing, as may occur in in vivo systems in the absence of transcriptional inhibitors.

Experimental Procedure for Small RNA-Seq Analysis

sRNAs were size fractionated from 10 μg of total RNA using PAGE on 15%, 8 M TBU gels, by excising the region from ∼18 to 35 bp. Sequencing libraries were prepared using the NEB Next sRNA kit (E7300) according to the manufacturer’s instructions from 5 μg of isolated sRNA. Ligations were performed for 1 h and pilot PCRs (20 μL) were performed on the cDNA to determine an optimal PCR cycle number of 13, and final PCR was performed using a 40 μL reaction in duplicate. Amplified libraries were size selected by PAGE on 6% ×ばつ TBE gels. Libraries were prepared in biological triplicate, indexed, pooled, and sequenced over two lanes of the Illumina HiSeq 2000 in 50-bp single-end mode at the BRF.

Read quality control was performed with FastQC v0.11.2 before and after trimming. Given that reads were sequenced with 50-bp read lengths, all sRNAs were expected to contain adapters on the 3ʹ. Adapters were trimmed using cutadapt v1.8 with the flags: -e 0.1, first nine bases of adapter must match perfectly; -m 18 -M 25, only keep reads lengths from 18 to 25 nucleotides after trimming; -O 10, adapter must overlap by at least 10 bases; -a AGATCGGAAGAGC, adapter sequence (Python v2.7.8; Martin, 2011; note that the cutadapt algorithm has changed significantly since publication). Histograms of read-length distribution were calculated using textHistogram confirming enrichment for 21- and 24-nucleotide sRNAs as expected. Reads were then aligned to the Arabidopsis genome (TAIR10) using Bowtie2 v2.2.5 (Langmead and Salzberg, 2012) using the flags: -a to report all matches for multi-mapped reads; -D 20 and -R 3, increases the likelihood that Bowtie2 will report the correct alignment for a read that aligns many places and -i S,1,0.50, reduces the substring interval, further increasing sensitivity; -end-to-end, preventing trimming of reads to enable alignment; -L 10 reduces substring length to 10 (default 22) as these are short reads and -N 0 requires exact match in the seed; -score-min L,0,0 reports only exact matches in -end-to-end mode (alignment score of 0 required which is max possible in end mode). Reads were then sorted, indexed, and compressed using samtools v1.1-26-g29b0367 (Li et al., 2009a) and strand-specific bigwig files were generated using bedtools genomecov v2.16.1 (Quinlan and Hall, 2010) and the UCSC utility bedGraphToBigWig for viewing in IGV (Robinson et al., 2011). To identify regulatory siRNAs clusters, de novo annotation and quantification of siRNAs were also performed using ShortStack v3.3.3 (Axtell, 2013; Shahid and Axtell, 2014). Here, trimmed sRNAs of 7 to 34 nucleotides were mapped to TAIR10 using the ShortStack wrappers with default settings except to make cluster sensitivity more consistent across samples the–mincov flag was used, set to 0.5 RPM, equivalent to around seven reads per library. siRNA clusters were overlapped with TAIR10 gene and TE annotations, downloaded from TAIR, using bedtools intersect (v2.24.0) with default parameters.

PARE

PARE libraries were prepared as per Zhai et al. (2014) with minor modifications included below. PARE libraries capture a 20- to 21-nucleotide tag of the 5ʹ end of all uncapped (5ʹ monophosphate) mRNA molecules by ligation of a 5ʹ adapter incorporating the recognition site of the Type II restriction endonuclease MmeI. Libraries were prepared from poly(A)⁺ RNA captured from 75 μg of total RNA. Following ligation of the double-stranded DNA 3ʹ adapter, the ligation product was PAGE purified to remove abundant adapter dimers, yielding higher quality final libraries. Following tag generation, libraries were amplified with 3ʹ index primers from the Illumina Truseq sRNA Kit (RS-200-0012) and a custom PARE-TruSeq primer, facilitating incorporation into the Illumina workflow for sequencing. Pilot PCRs were employed to determine an optimal PCR cycle number of eight cycles for the final amplification, which was performed as a half volume reaction (25 μL). Following the final PCR, the indexed libraries were pooled equal-molar then purified using AMPure beads using a bead ratio of 1:1.7 determined empirically using a standard curve to maximize recovery of the 130-bp library and minimize carryover of 80-bp contaminants. Libraries were sequenced on the Illumina HiSeq 2500 in 50-bp single-end sequencing over one lane at the BRF using a custom read 1 primer (5ʹ-CCACCGACAGGTTCAGAGTTCTACAGTCCGAC-3ʹ, PAGE purified).

FastQC v0.11.2 was used to examine fastq files for quality statistics; typically, read quality is exceptional for bases 1 to 21 of all reads (Q > 32, 99.9% accuracy). The PARE protocol produces very consistent fragments of 20 to 21 nucleotides; therefore, using longer read technology (50 or 100 bp) leads to low quality sequencing of the 3ʹ adapter owing to color imbalance. To trim this often low-quality adapter sequence, reads were first hard trimmed to 20 nucleotides using cutadapt, then the remaining 20-nucleotide fragments were scanned for adapter contamination using cutadapt v1.8 with the flags: –cut -31, trim 31 bases from the 3ʹ end; -e 0.1, first nine bases of adapter must match perfectly; -m 14 -M 22 only keep reads between 14 and 22 nucleotides after trimming; -O 3, adapter must overlap by at least three bases; -a TGGAATTCTC, Illumina small RNA TruSeq 5ʹ adapter sequence. This approach gave a higher yield of reads passing filtering and ultimately mapping to the genome compared with quality trimming or adapter removal without a fixed length trim first (>90% alignment rate). Reads were then aligned to the Arabidopsis TAIR10 genome using Bowtie2 v2.2.5 (Langmead and Salzberg, 2012) using the same settings as described above for the sRNA analysis. The PARE data from Hou et al. (2016) were obtained from SRA (SRR3143654) and processed in the same manner above. GMUCT data from Yu et al. (2016) (SRR863540 and SRR863541) were obtained from SRA and mapped using subjunc v.1.4.6 owing to the longer read length of 50 bp. To analyze ribosome footprints in PARE data, we implemented a novel tool PAREphase, which records the position of the 5ʹ end of each aligned read relative to the stop codon. PAREphase is implemented in Python3 and produces gene-wise summaries in a tabular form from reads in BAM format and an annotated reference genome in FASTA/GFF3 format. PAREphase version 1.0 was used for all analyses; PAREphase is available at https://github.com/kdmurray91/PAREphase. The phased CDI was determined by considering counts between –95 to –19 and taking the sum of the counts in frame 1 (protected frame) and dividing by the average of the sum of counts in frames 0 and 2. The stalled CDI was determined using the sum of counts between –18 and –16 (stalled codon) divided by the average counts of all codons between 0 and +47 of the untranslated region, which represent background PARE coverage.

WGBS Analysis

WGBS analysis was performed on four replicates for the time points: control (I), 60-min excess light (III), and 24 h recovery (X); however, these samples were not tissue paired with the RNA sequencing analysis. DNA was extracted using the Qiagen DNeasy Plant Mini Kit as per the manufacturer’s instructions. In total, 1000 ng of genomic DNA was fragmented using a Covaris column purified using the MinElute PCR Cleanup (Qiagen), end repaired using the End-It Kit (Epicentre), A-tailed then bead cleaned using AMPure XP beads (Beckman Coulter). NEXTflex methylated adapters (Bioo Scientific) were ligated and following bead cleanup the ligation product was quantified on a Qubit (Life Technologies). Then, 450 ng of clean ligation was bisulfite converted using the MethylCode kit (Life Technologies). Bisulfite -converted DNA was amplified with six cycles of PCR using a Kapa HiFi HotStart U+ ×ばつ readymix (Kapa Biosystems) and an Illumina compatible NEXTflex Primer mix (Bioo Scientific). Samples were bead-cleaned, quantified by qPCR using a Library Quantification Kit for Illumina/LightCycler 480 (Kapa Biosystems), pooled equal-molar, and six-sample pools were sequenced per lane on a HiSeq 2500 (100 bp single end) at the BRF; we saw no evidence of lane effects (Supplemental Figure 10).

Raw sequencing reads were quality filtered using Trim Galore! v0.3.7 (using Cutadapt v1.9). Reads were the aligned to TAIR10 using Bismark v0.14.5 (Krueger and Andrews, 2011) with the flags –bowtie1 -n 2, -l 20. Methylated cytosines were extracted from aligned reads using the Bismark methylation extractor with default parameters. Bisulfite conversion rate was calculated from the proportion of unconverted cytosines in the CHH context from the chloroplast genome and ranged from 99.4 to 99.7% per sample. The proportion of CG, CHG, and CHH methylation was determined as mean methylation across reads at single cytosine resolution and across 100-bp windows for genome-wide comparisons. Correlations (R²) between methylation levels were performed on the 100-bp bin mean methylation levels. DMRs were identified using Bayesian modeling, across individual cytosines, in R (v3.2.5) using the package DSS (v2.10.0; Feng et al., 2014) following the recommended default settings except for a reduced smoothing window size (smoothing.span = 100). The threshold methylation difference for DMRs in each sequence context (delta) was 0.4 for CG, 0.2 for CHG, and 0.2 for CHH. DSS calculates an FDR-adjusted P value and the posterior probability that the difference specified (delta) is significant; DMRs were considered significant at FDR < 0.05. Methylation levels were assigned to annotated gene and transposable element features of the TAIR10 assembly using Bedtools (v2.21.0; Quinlan and Hall, 2010).

SNPs were identified using the WGBS data and BS-SNPer (v1.0; Gao et al., 2015). Samples were compared with the TAIR10 Col-0 reference using default parameters with the following modifications, minhetfreq = 0.15, minhomfreq = 0.85, Minquali = 30, and Mincover = 15. Low-confidence SNPs were filtered out by removing thousands of C-T and G-A transition mutations (mostly likely caused by bisulfite conversion); SNPs not given a PASS for FILTER status and SNPs with a Phred-scaled quality score < 30, leaving on average 20 putative SNPs per sample. Distance correlations and hierarchical clustering of samples (using "complete" distance) was performed in R.

Accession Numbers

RNA-seq, sRNA-seq, PARE, and WGBS reads are deposited with the NCBI Short Read Archive under BioProject accession number PRJNA391262. Bioinformatics analysis pipelines are available online on GitHub (https://github.com/pedrocrisp/NGS-pipelines). Accession numbers are as follows: HSP101 (AT1G74310), AMT1;2 (AT1G64780), APX2 (AT3G09640), HSP20like1 (AT1G59860), HSP20like2 (AT1G54050), RRTF1 (AT4G34410), HSF.A7A (AT3G51910), ZAT10/STZ (AT1G27730), HSP70 (AT3G12580), RCA (AT2G39730), and LHCB1.4 (AT2G34430).2

Supplemental Data

• Supplemental Figure 1. Characterization of warm excess-light stress.

• Supplemental Figure 2. RNA-seq stress and recovery time-course sample clustering and qPCR corroboration.

• Supplemental Figure 3. A selection of example excess-light response and recovery profiles for individual transcripts.

• Supplemental Figure 4. Excess-light stress in Arabidopsis protoplasts.

• Supplemental Figure 5. Comparison of RNA decay models using simulated data.

• Supplemental Figure 6. Degradome abundance correlates with mRNA abundance but not half-life.

• Supplemental Figure 7. PARE QC and cotranslational decay profiling in published data sets.

• Supplemental Figure 8. Additional gene expression profiling during recovery in 5ʹ RNA decay mutants.

• Supplemental Figure 9. Stress recovery in rdr mutants.

• Supplemental Figure 10. Excess-light stress DNA methylome and epitypes.

• Supplemental Figure 11. Transcriptional priming and memory types.

• Supplemental Methods.

• Supplemental Data Set 1. List of genes differentially expressed during drought determined by RNA-seq.

• Supplemental Data Set 2. List of genes differentially expressed during excess-light stress determined by RNA-seq.

• Supplemental Data Set 3. RNA-seq relative gene expression for the excess-light stress and recovery time course.

• Supplemental Data Set 4. Individual stress and recovery plots for 1000 mRNAs.

• Supplemental Data Set 5. Identification and half-lives of RRGD loci among excess-light upregulated genes.

• Supplemental Data Set 6. List of the 1000 most rapid decay events and corresponding half-lives.

• Supplemental Data Set 7. Enriched GO terms among the 1000 most unstable transcripts.

• Supplemental Data Set 8. Total numbers and genomic locations of siRNA clusters.

• Supplemental Data Set 9. siRNA clusters overlapping excess-light upregulated transcripts.

• Supplemental Data Set 10. Summary epitype DMRs and their overlap with spontaneous epiallele loci.

• Supplemental Data Set 11. Genomic coordinates of epitype DMRs.

• Supplemental Data Set 12. SNPs between Col-0 plants used for WGBS.

• Supplemental Data Set 13. Excess-light stress associated DMRs.

• Supplemental Data Set 14. List of transcriptional memory loci.

• Supplemental Data Set 15. Enriched GO terms among recovery activated transcripts.

• Supplemental Data Set 16. Primers used in this study.

• Supplemental Data Set 17. Summary of sequencing experiments.

Glossary

NPQ

nonphotochemical quenching

siRNA

small interfering RNA

RDR

RNA-dependent RNA polymerases

PTGS

posttranscriptional gene silencing

RWC

relative water content

FDR

false discovery rate

padj

adjusted P value

PARE

parallel analysis of RNA ends

WGBS

whole-genome bisulfite sequencing

CPM

counts per million

RRGD

rapid recovery gene downregulation

GO

Gene Ontology

CDI

cotranslational decay index

TE

transposable element

SNPs

single nucleotide polymorphisms

DMR

differentially methylated region

GLM

generalized linear model

miRNA

microRNA

1. Addo-Quaye C., Eshoo T.W., Bartel D.P., Axtell M.J. (2008). Endogenous siRNA and miRNA targets identified by sequencing of the Arabidopsis degradome. Curr. Biol. 18: 758–762. [DOI] [PMC free article] [PubMed] [Google Scholar]
2. Allen R.S., Li J., Alonso-Peral M.M., White R.G., Gubler F., Millar A.A. (2010). MicroR159 regulation of most conserved targets in Arabidopsis has negligible phenotypic effects. Silence 1: 18. [DOI] [PMC free article] [PubMed] [Google Scholar]
3. Amorim M.J., Cotobal C., Duncan C., Mata J. (2010). Global coordination of transcriptional control and mRNA decay during cellular differentiation. Mol. Syst. Biol. 6: 380. [DOI] [PMC free article] [PubMed] [Google Scholar]
4. Apel K., Hirt H. (2004). Reactive oxygen species: metabolism, oxidative stress, and signal transduction. Annu. Rev. Plant Biol. 55: 373–399. [DOI] [PubMed] [Google Scholar]
5. Armbruster U., Carrillo L.R., Venema K., Pavlovic L., Schmidtmann E., Kornfeld A., Jahns P., Berry J.A., Kramer D.M., Jonikas M.C. (2014). Ion antiport accelerates photosynthetic acclimation in fluctuating light environments. Nat. Commun. 5: 5439. [DOI] [PMC free article] [PubMed] [Google Scholar]
6. Avramova Z. (2015). Transcriptional ‘memory’ of a stress: transient chromatin and memory (epigenetic) marks at stress-response genes. Plant J. 83: 149–159. [DOI] [PubMed] [Google Scholar]
7. Axtell M.J. (2013). ShortStack: comprehensive annotation and quantification of small RNA genes. RNA 19: 740–751. [DOI] [PMC free article] [PubMed] [Google Scholar]
8. Baker N.R. (2008). Chlorophyll fluorescence: a probe of photosynthesis in vivo. Annu. Rev. Plant Biol. 59: 89–113. [DOI] [PubMed] [Google Scholar]
9. Becker C., Hagmann J., Müller J., Koenig D., Stegle O., Borgwardt K., Weigel D. (2011). Spontaneous epigenetic variation in the Arabidopsis thaliana methylome. Nature 480: 245–249. [DOI] [PubMed] [Google Scholar]
10. Belostotsky D.A., Sieburth L.E. (2009). Kill the messenger: mRNA decay and plant development. Curr. Opin. Plant Biol. 12: 96–102. [DOI] [PubMed] [Google Scholar]
11. Boyko A., Kovalchuk I. (2011). Genome instability and epigenetic modification--heritable responses to environmental stress? Curr. Opin. Plant Biol. 14: 260–266. [DOI] [PubMed] [Google Scholar]
12. Boyle E.I., Weng S., Gollub J., Jin H., Botstein D., Cherry J.M., Sherlock G. (2004). GO:TermFinder--open source software for accessing Gene Ontology information and finding significantly enriched Gene Ontology terms associated with a list of genes. Bioinformatics 20: 3710–3715. [DOI] [PMC free article] [PubMed] [Google Scholar]
13. Castells-Roca L., García-Martínez J., Moreno J., Herrero E., Bellí G., Pérez-Ortín J.E. (2011). Heat shock response in yeast involves changes in both transcription rates and mRNA stabilities. PLoS One 6: e17272. [DOI] [PMC free article] [PubMed] [Google Scholar]
14. Chan K.X., Phua S.Y., Crisp P., McQuinn R., Pogson B.J. (2016). Learning the languages of the chloroplast: Retrograde signaling and beyond. Annu. Rev. Plant Biol. 67: 25–53. [DOI] [PubMed] [Google Scholar]
15. Chekanova J.A., et al. (2007). Genome-wide high-resolution mapping of exosome substrates reveals hidden features in the Arabidopsis transcriptome. Cell 131: 1340–1353. [DOI] [PubMed] [Google Scholar]
16. Chiba Y., Mineta K., Hirai M.Y., Suzuki Y., Kanaya S., Takahashi H., Onouchi H., Yamaguchi J., Naito S. (2013). Changes in mRNA stability associated with cold stress in Arabidopsis cells. Plant Cell Physiol. 54: 180–194. [DOI] [PubMed] [Google Scholar]
17. Christie M., Brosnan C.A., Rothnagel J.A., Carroll B.J. (2011). RNA decay and RNA silencing in plants: competition or collaboration? Front. Plant Sci. 2: 99. [DOI] [PMC free article] [PubMed] [Google Scholar]
18. Coleman-Derr D., Zilberman D. (2012). Deposition of histone variant H2A.Z within gene bodies regulates responsive genes. PLoS Genet. 8: e1002988. [DOI] [PMC free article] [PubMed] [Google Scholar]
19. Crisp P.A., Ganguly D., Eichten S.R., Borevitz J.O., Pogson B.J. (2016). Reconsidering plant memory: Intersections between stress recovery, RNA turnover, and epigenetics. Sci. Adv. 2: e1501340. [DOI] [PMC free article] [PubMed] [Google Scholar]
20. Demmig-Adams B., Adams W. (1992). Photoprotection and other responses of plants to high light stress. Annu. Rev. Plant Physiol. Plant Mol. Biol. 43: 599–626. [Google Scholar]
21. Dietz K.-J. (2015). Efficient high light acclimation involves rapid processes at multiple mechanistic levels. J. Exp. Bot. 66: 2401–2414. [DOI] [PubMed] [Google Scholar]
22. Ding J., Li D., Ohler U., Guan J., Zhou S. (2012a). Genome-wide search for miRNA-target interactions in Arabidopsis thaliana with an integrated approach. BMC Genomics 13 (suppl 3): S3. [DOI] [PMC free article] [PubMed] [Google Scholar]
23. Ding Y., Fromm M., Avramova Z. (2012b). Multiple exposures to drought ‘train’ transcriptional responses in Arabidopsis. Nat. Commun. 3: 740. [DOI] [PubMed] [Google Scholar]
24. Ding Y., Liu N., Virlouvet L., Riethoven J.-J., Fromm M., Avramova Z. (2013). Four distinct types of dehydration stress memory genes in Arabidopsis thaliana. BMC Plant Biol. 13: 229. [DOI] [PMC free article] [PubMed] [Google Scholar]
25. Eichten S.R., Schmitz R.J., Springer N.M. (2014). Epigenetics: beyond chromatin modifications and complex genetic regulation. Plant Physiol. 165: 933–947. [DOI] [PMC free article] [PubMed] [Google Scholar]
26. Elkon R., Zlotorynski E., Zeller K.I., Agami R. (2010). Major role for mRNA stability in shaping the kinetics of gene induction. BMC Genomics 11: 259. [DOI] [PMC free article] [PubMed] [Google Scholar]
27. Erhard K.F. Jr., Talbot J.-E.R.B., Deans N.C., McClish A.E., Hollick J.B. (2015). Nascent transcription affected by RNA polymerase IV in Zea mays. Genetics 199: 1107–1125. [DOI] [PMC free article] [PubMed] [Google Scholar]
28. Estavillo G.M., et al. (2011). Evidence for a SAL1-PAP chloroplast retrograde pathway that functions in drought and high light signaling in Arabidopsis. Plant Cell 23: 3992–4012. [DOI] [PMC free article] [PubMed] [Google Scholar]
29. Fan J., Yang X., Wang W., Wood W.H. III, Becker K.G., Gorospe M. (2002). Global analysis of stress-regulated mRNA turnover by using cDNA arrays. Proc. Natl. Acad. Sci. USA 99: 10611–10616. [DOI] [PMC free article] [PubMed] [Google Scholar]
30. Fei Q., Xia R., Meyers B.C. (2013). Phased, secondary, small interfering RNAs in posttranscriptional regulatory networks. Plant Cell 25: 2400–2415. [DOI] [PMC free article] [PubMed] [Google Scholar]
31. Feng H., Conneely K.N., Wu H. (2014). A Bayesian hierarchical model to detect differentially methylated loci from single nucleotide resolution sequencing data. Nucleic Acids Res. 42: e69. [DOI] [PMC free article] [PubMed] [Google Scholar]
32. Frei dit Frey N., Muller P., Jammes F., Kizis D., Leung J., Perrot-Rechenmann C., Bianchi M.W. (2010). The RNA binding protein Tudor-SN is essential for stress tolerance and stabilizes levels of stress-responsive mRNAs encoding secreted proteins in Arabidopsis. Plant Cell 22: 1575–1591. [DOI] [PMC free article] [PubMed] [Google Scholar]
33. Galvez-Valdivieso G., Fryer M.J., Lawson T., Slattery K., Truman W., Smirnoff N., Asami T., Davies W.J., Jones A.M., Baker N.R., Mullineaux P.M. (2009). The high light response in Arabidopsis involves ABA signaling between vascular and bundle sheath cells. Plant Cell 21: 2143–2162. [DOI] [PMC free article] [PubMed] [Google Scholar]
34. Gao S., et al. (2015). BS-SNPer: SNP calling in bisulfite-seq data. Bioinformatics 31: 4006–4008. [DOI] [PMC free article] [PubMed] [Google Scholar]
35. Gazzani S., Lawrenson T., Woodward C., Headon D., Sablowski R. (2004). A link between mRNA turnover and RNA interference in Arabidopsis. Science 306: 1046–1048. [DOI] [PubMed] [Google Scholar]
36. German M.A., et al. (2008). Global identification of microRNA-target RNA pairs by parallel analysis of RNA ends. Nat. Biotechnol. 26: 941–946. [DOI] [PubMed] [Google Scholar]
37. Gregory B.D., O’Malley R.C., Lister R., Urich M.A., Tonti-Filippini J., Chen H., Millar A.H., Ecker J.R. (2008). A link between RNA metabolism and silencing affecting Arabidopsis development. Dev. Cell 14: 854–866. [DOI] [PubMed] [Google Scholar]
38. Gutierrez R.A., Ewing R.M., Cherry J.M., Green P.J. (2002). Identification of unstable transcripts in Arabidopsis by cDNA microarray analysis: rapid decay is associated with a group of touch- and specific clock-controlled genes. Proc. Natl. Acad. Sci. USA 99: 11513–11518. [DOI] [PMC free article] [PubMed] [Google Scholar]
39. Gutzat R., Mittelsten Scheid O. (2012). Epigenetic responses to stress: triple defense? Curr. Opin. Plant Biol. 15: 568–573. [DOI] [PMC free article] [PubMed] [Google Scholar]
40. Gy I., Gasciolli V., Lauressergues D., Morel J.B., Gombert J., Proux F., Proux C., Vaucheret H., Mallory A.C. (2007). Arabidopsis FIERY1, XRN2, and XRN3 are endogenous RNA silencing suppressors. Plant Cell 19: 3451–3461. [DOI] [PMC free article] [PubMed] [Google Scholar]
41. Herman J.J., Sultan S.E. (2011). Adaptive transgenerational plasticity in plants: case studies, mechanisms, and implications for natural populations. Front. Plant Sci. 2: 102. [DOI] [PMC free article] [PubMed] [Google Scholar]
42. Hetzel J., Duttke S.H., Benner C., Chory J. (2016). Nascent RNA sequencing reveals distinct features in plant transcription. Proc. Natl. Acad. Sci. USA 113: 12316–12321. [DOI] [PMC free article] [PubMed] [Google Scholar]
43. Hirsch J., et al. (2011). A novel fry1 allele reveals the existence of a mutant phenotype unrelated to 5′->3′ exoribonuclease (XRN) activities in Arabidopsis thaliana roots. PLoS One 6: e16724. [DOI] [PMC free article] [PubMed] [Google Scholar]
44. Horton P., Johnson M.P., Perez-Bueno M.L., Kiss A.Z., Ruban A.V. (2008). Photosynthetic acclimation: does the dynamic structure and macro-organisation of photosystem II in higher plant grana membranes regulate light harvesting states? FEBS J. 275: 1069–1079. [DOI] [PubMed] [Google Scholar]
45. Hou C.-Y., Lee W.-C., Chou H.-C., Chen A.-P., Chou S.-J., Chen H.-M. (2016). Global analysis of truncated RNA ends reveals new insights into ribosome stalling in plants. Plant Cell 28: 2398–2416. [DOI] [PMC free article] [PubMed] [Google Scholar]
46. Hou C.-Y., Wu M.-T., Lu S.-H., Hsing Y.-I., Chen H.-M. (2014). Beyond cleaved small RNA targets: unraveling the complexity of plant RNA degradome data. BMC Genomics 15: 15. [DOI] [PMC free article] [PubMed] [Google Scholar]
47. Iwasaki M., Paszkowski J. (2014). Identification of genes preventing transgenerational transmission of stress-induced epigenetic states. Proc. Natl. Acad. Sci. USA 111: 8547–8552. [DOI] [PMC free article] [PubMed] [Google Scholar]
48. Jiao Y., Riechmann J.L., Meyerowitz E.M. (2008). Transcriptome-wide analysis of uncapped mRNAs in Arabidopsis reveals regulation of mRNA degradation. Plant Cell 20: 2571–2585. [DOI] [PMC free article] [PubMed] [Google Scholar]
49. Jung H.-S., Crisp P.A., Estavillo G.M., Cole B., Hong F., Mockler T.C., Pogson B.J., Chory J. (2013). Subset of heat-shock transcription factors required for the early response of Arabidopsis to excess light. Proc. Natl. Acad. Sci. USA 110: 14474–14479. [DOI] [PMC free article] [PubMed] [Google Scholar]
50. Kamalay J.C., Goldberg R.B. (1980). Regulation of structural gene expression in tobacco. Cell 19: 935–946. [DOI] [PubMed] [Google Scholar]
51. Karpinski S., Escobar C., Karpinska B., Creissen G., Mullineaux P.M. (1997). Photosynthetic electron transport regulates the expression of cytosolic ascorbate peroxidase genes in Arabidopsis during excess light stress. Plant Cell 9: 627–640. [DOI] [PMC free article] [PubMed] [Google Scholar]
52. Kasahara M., Kagawa T., Oikawa K., Suetsugu N., Miyao M., Wada M. (2002). Chloroplast avoidance movement reduces photodamage in plants. Nature 420: 829–832. [DOI] [PubMed] [Google Scholar]
53. Kiper M., Bartels D., Herzfeld F., Richter G. (1979). The expression of a plant genome in hnRNA and mRNA. Nucleic Acids Res. 6: 1961–1978. [DOI] [PMC free article] [PubMed] [Google Scholar]
54. Krueger F., Andrews S.R. (2011). Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics 27: 1571–1572. [DOI] [PMC free article] [PubMed] [Google Scholar]
55. Lam L.T., et al. (2001). Genomic-scale measurement of mRNA turnover and the mechanisms of action of the anti-cancer drug flavopiridol. Genome Biol. 2: research0041.1–research0041.11. [DOI] [PMC free article] [PubMed]
56. Langmead B., Salzberg S.L. (2012). Fast gapped-read alignment with Bowtie 2. Nat. Methods 9: 357–359. [DOI] [PMC free article] [PubMed] [Google Scholar]
57. Lawlor D.W. (2013). Genetic engineering to improve plant performance under drought: physiological evaluation of achievements, limitations, and possibilities. J. Exp. Bot. 64: 83–108. [DOI] [PubMed] [Google Scholar]
58. Lee J.-Y., Lee D.-H. (2003). Use of serial analysis of gene expression technology to reveal changes in gene expression in Arabidopsis pollen undergoing cold stress. Plant Physiol. 132: 517–529. [DOI] [PMC free article] [PubMed] [Google Scholar]
59. Li H., Handsaker B., Wysoker A., Fennell T., Ruan J., Homer N., Marth G., Abecasis G., Durbin R.; 1000 Genome Project Data Processing Subgroup (2009a). The Sequence Alignment/Map format and SAMtools. Bioinformatics 25: 2078–2079. [DOI] [PMC free article] [PubMed] [Google Scholar]
60. Li Z., Wakao S., Fischer B.B., Niyogi K.K. (2009b). Sensing and responding to excess light. Annu. Rev. Plant Biol. 60: 239–260. [DOI] [PubMed] [Google Scholar]
61. Liao Y., Smyth G.K., Shi W. (2014). FeatureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30: 923–930. [DOI] [PubMed] [Google Scholar]
62. Liao Y., Smyth G.K., Shi W. (2013). The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote. Nucleic Acids Res. 41: e108. [DOI] [PMC free article] [PubMed] [Google Scholar]
63. Lichtenthaler H.K., Buschmann C., Knapp M. (2005). How to correctly determine the different chlorophyll fluorescence parameters and the chlorophyll fluorescence decrease ratio RFd of leaves with the PAM fluorometer. Photosynthetica 43: 379–393. [Google Scholar]
64. Luo Z., Chen Z. (2007). Improperly terminated, unpolyadenylated mRNA of sense transgenes is targeted by RDR6-mediated RNA silencing in Arabidopsis. Plant Cell 19: 943–958. [DOI] [PMC free article] [PubMed] [Google Scholar]
65. Madsen J.G.S., Schmidt S.F., Larsen B.D., Loft A., Nielsen R., Mandrup S. (2015). iRNA-seq: computational method for genome-wide assessment of acute transcriptional regulation from total RNA-seq data. Nucleic Acids Res. 43: e40. [DOI] [PMC free article] [PubMed] [Google Scholar]
66. Marinov G.K., Williams B.A., McCue K., Schroth G.P., Gertz J., Myers R.M., Wold B.J. (2014). From single-cell to cell-pool transcriptomes: stochasticity in gene expression and RNA splicing. Genome Res. 24: 496–510. [DOI] [PMC free article] [PubMed] [Google Scholar]
67. Martin M. (2011). Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal 17: 10–12. [Google Scholar]
68. Martínez de Alba A.E., Moreno A.B., Gabriel M., Mallory A.C., Christ A., Bounon R., Balzergue S., Aubourg S., Gautheret D., Crespi M.D., Vaucheret H., Maizel A. (2015). In plants, decapping prevents RDR6-dependent production of small interfering RNAs from endogenous mRNAs. Nucleic Acids Res. 43: 2902–2913. [DOI] [PMC free article] [PubMed] [Google Scholar]
69. Maxwell K., Johnson G.N. (2000). Chlorophyll fluorescence--a practical guide. J. Exp. Bot. 51: 659–668. [DOI] [PubMed] [Google Scholar]
70. McCarthy D.J., Chen Y., Smyth G.K. (2012). Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 40: 4288–4297. [DOI] [PMC free article] [PubMed] [Google Scholar]
71. Merret R., Descombin J., Juan Y.T., Favory J.-J., Carpentier M.-C., Chaparro C., Charng Y.Y., Deragon J.-M., Bousquet-Antonelli C. (2013). XRN4 and LARP1 are required for a heat-triggered mRNA decay pathway involved in plant acclimation and survival during thermal stress. Cell Reports 5: 1279–1293. [DOI] [PubMed] [Google Scholar]
72. Merret R., Nagarajan V.K., Carpentier M.-C., Park S., Favory J.-J., Descombin J., Picart C., Charng Y.Y., Green P.J., Deragon J.-M., Bousquet-Antonelli C. (2015). Heat-induced ribosome pausing triggers mRNA co-translational decay in Arabidopsis thaliana. Nucleic Acids Res. 43: 4121–4132. [DOI] [PMC free article] [PubMed] [Google Scholar]
73. Mickelbart M.V., Hasegawa P.M., Bailey-Serres J. (2015). Genetic mechanisms of abiotic stress tolerance that translate to crop yield stability. Nat. Rev. Genet. 16: 237–251. [DOI] [PubMed] [Google Scholar]
74. Mirouze M., Paszkowski J. (2011). Epigenetic contribution to stress adaptation in plants. Curr. Opin. Plant Biol. 14: 267–274. [DOI] [PubMed] [Google Scholar]
75. Moldovan D., Spriggs A., Yang J., Pogson B.J., Dennis E.S., Wilson I.W. (2010). Hypoxia-responsive microRNAs and trans-acting small interfering RNAs in Arabidopsis. J. Exp. Bot. 61: 165–177. [DOI] [PMC free article] [PubMed] [Google Scholar]
76. Molin C., Jauhiainen A., Warringer J., Nerman O., Sunnerhagen P. (2009). mRNA stability changes precede changes in steady-state mRNA amounts during hyperosmotic stress. RNA 15: 600–614. [DOI] [PMC free article] [PubMed] [Google Scholar]
77. Molina-Navarro M.M., Castells-Roca L., Bellí G., García-Martínez J., Marín-Navarro J., Moreno J., Pérez-Ortín J.E., Herrero E. (2008). Comprehensive transcriptional analysis of the oxidative response in yeast. J. Biol. Chem. 283: 17908–17918. [DOI] [PubMed] [Google Scholar]
78. Mortazavi A., Williams B.A., McCue K., Schaeffer L., Wold B. (2008). Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 5: 621–628. [DOI] [PubMed] [Google Scholar]
79. Müller P., Li X.-P., Niyogi K.K. (2001). Non-photochemical quenching. A response to excess light energy. Plant Physiol. 125: 1558–1566. [DOI] [PMC free article] [PubMed] [Google Scholar]
80. Murchie E.H., Niyogi K.K. (2011). Manipulation of photoprotection to improve plant photosynthesis. Plant Physiol. 155: 86–92. [DOI] [PMC free article] [PubMed] [Google Scholar]
81. Murchie E.H., Pinto M., Horton P. (2009). Agriculture and the new challenges for photosynthesis research. New Phytol. 181: 532–552. [DOI] [PubMed] [Google Scholar]
82. Narsai R., Howell K.A., Millar A.H., O’Toole N., Small I., Whelan J. (2007). Genome-wide analysis of mRNA decay rates and their determinants in Arabidopsis thaliana. Plant Cell 19: 3418–3436. [DOI] [PMC free article] [PubMed] [Google Scholar]
83. Nguyen A.H., Matsui A., Tanaka M., Mizunashi K., Nakaminami K., Hayashi M., Iida K., Toyoda T., Nguyen D.V., Seki M. (2015). Loss of Arabidopsis 5′-3′ exoribonuclease AtXRN4 function enhances heat stress tolerance of plants subjected to severe heat stress. Plant Cell Physiol. 56: 1762–1772. [DOI] [PubMed] [Google Scholar]
84. Nilkens M., Kress E., Lambrev P., Miloslavina Y., Müller M., Holzwarth A.R., Jahns P. (2010). Identification of a slowly inducible zeaxanthin-dependent component of non-photochemical quenching of chlorophyll fluorescence generated under steady-state conditions in Arabidopsis. Biochim. Biophys. Acta 1797: 466–475. [DOI] [PubMed] [Google Scholar]
85. Ort D.R. (2001). When there is too much light. Plant Physiol. 125: 29–32. [DOI] [PMC free article] [PubMed] [Google Scholar]
86. Papadakis A.K., Roubelakis-Angelakis K.A. (2002). Oxidative stress could be responsible for the recalcitrance of plant protoplasts. Plant Physiol. Biochem. 40: 549–559. [Google Scholar]
87. Park S.-H., Chung P.J., Juntawong P., Bailey-Serres J., Kim Y.S., Jung H., Bang S.W., Kim Y.-K., Do Choi Y., Kim J.-K. (2012). Posttranscriptional control of photosynthetic mRNA decay under stress conditions requires 3′ and 5′ untranslated regions and correlates with differential polysome association in rice. Plant Physiol. 159: 1111–1124. [DOI] [PMC free article] [PubMed] [Google Scholar]
88. Parkhomchuk D., Borodina T., Amstislavskiy V., Banaru M., Hallen L., Krobitsch S., Lehrach H., Soldatov A. (2009). Transcriptome analysis by strand-specific sequencing of complementary DNA. Nucleic Acids Res. 37: e123. [DOI] [PMC free article] [PubMed] [Google Scholar]
89. Pastenes C., Porter V., Baginsky C., Horton P., González J. (2004). Paraheliotropism can protect water-stressed bean (Phaseolus vulgaris L.) plants against photoinhibition. J. Plant Physiol. 161: 1315–1323. [DOI] [PubMed] [Google Scholar]
90. Pérez-Ortín J.E., Alepuz P.M., Moreno J. (2007). Genomics and gene transcription kinetics in yeast. Trends Genet. 23: 250–257. [DOI] [PubMed] [Google Scholar]
91. Petracek M.E., Dickey L.F., Nguyen T.T., Gatz C., Sowinski D.A., Allen G.C., Thompson W.F. (1998). Ferredoxin-1 mRNA is destabilized by changes in photosynthetic electron transport. Proc. Natl. Acad. Sci. USA 95: 9009–9013. [DOI] [PMC free article] [PubMed] [Google Scholar]
92. Quinlan A.R., Hall I.M. (2010). BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26: 841–842. [DOI] [PMC free article] [PubMed] [Google Scholar]
93. Ramakers C., Ruijter J.M., Deprez R.H.L., Moorman A.F.M. (2003). Assumption-free analysis of quantitative real-time polymerase chain reaction (PCR) data. Neurosci. Lett. 339: 62–66. [DOI] [PubMed] [Google Scholar]
94. Raven J. (1994). The cost of photoinhibition to plant communities. In Photoinhibition of Photosynthesis – Molecular Mechanisms to the Field, Baker N., Bowyer J., eds (Oxford, UK: Bios Scientific; ), pp. 449–464. [Google Scholar]
95. Reichel M., Liao Y., Rettel M., Ragan C., Evers M., Alleaume A.-M., Horos R., Hentze M.W., Preiss T., Millar A.A. (2016). In planta determination of the mRNA-binding proteome of Arabidopsis etiolated seedlings. Plant Cell 28: 2435–2452. [DOI] [PMC free article] [PubMed] [Google Scholar]
96. Rivers J., Warthmann N., Pogson B.J., Borevitz J.O. (2015). Genomic breeding for food, environment and livelihoods. Food Secur. 7: 375–382. [Google Scholar]
97. Robinson J.T., Thorvaldsdóttir H., Winckler W., Guttman M., Lander E.S., Getz G., Mesirov J.P. (2011). Integrative genomics viewer. Nat. Biotechnol. 29: 24–26. [DOI] [PMC free article] [PubMed] [Google Scholar]
98. Robinson M.D., McCarthy D.J., Smyth G.K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26: 139–140. [DOI] [PMC free article] [PubMed] [Google Scholar]
99. Robinson M.D., Oshlack A. (2010). A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 11: R25. [DOI] [PMC free article] [PubMed] [Google Scholar]
100. Robinson M.D., Smyth G.K. (2007). Moderated statistical tests for assessing differences in tag abundance. Bioinformatics 23: 2881–2887. [DOI] [PubMed] [Google Scholar]
101. Robinson M.D., Smyth G.K. (2008). Small-sample estimation of negative binomial dispersion, with applications to SAGE data. Biostatistics 9: 321–332. [DOI] [PubMed] [Google Scholar]
102. Romero-Santacreu L., Moreno J., Pérez-Ortín J.E., Alepuz P. (2009). Specific and global regulation of mRNA stability during osmotic stress in Saccharomyces cerevisiae. RNA 15: 1110–1120. [DOI] [PMC free article] [PubMed] [Google Scholar]
103. Ross J. (1995). mRNA stability in mammalian cells. Microbiol. Rev. 59: 423–450. [DOI] [PMC free article] [PubMed] [Google Scholar]
104. Rossel J.B., Wilson P.B., Hussain D., Woo N.S., Gordon M.J., Mewett O.P., Howell K.A., Whelan J., Kazan K., Pogson B.J. (2007). Systemic and intracellular responses to photooxidative stress in Arabidopsis. Plant Cell 19: 4091–4110. [DOI] [PMC free article] [PubMed] [Google Scholar]
105. Ruckle M.E., Burgoon L.D., Lawrence L.A., Sinkler C.A., Larkin R.M. (2012). Plastids are major regulators of light signaling in Arabidopsis. Plant Physiol. 159: 366–390. [DOI] [PMC free article] [PubMed] [Google Scholar]
106. Ruijter J.M., Ramakers C., Hoogaars W.M.H., Karlen Y., Bakker O., van den Hoff M.J.B., Moorman A.F.M. (2009). Amplification efficiency: linking baseline and bias in the analysis of quantitative PCR data. Nucleic Acids Res. 37: e45. [DOI] [PMC free article] [PubMed] [Google Scholar]
107. Rymarquis L.A., Souret F.F., Green P.J. (2011). Evidence that XRN4, an Arabidopsis homolog of exoribonuclease XRN1, preferentially impacts transcripts with certain sequences or in particular functional categories. RNA 17: 501–511. [DOI] [PMC free article] [PubMed] [Google Scholar]
108. Sakamoto H., Maruyama K., Sakuma Y., Meshi T., Iwabuchi M., Shinozaki K., Yamaguchi-Shinozaki K. (2004). Arabidopsis Cys2/His2-type zinc-finger proteins function as transcription repressors under drought, cold, and high-salinity stress conditions. Plant Physiol. 136: 2734–2746. [DOI] [PMC free article] [PubMed] [Google Scholar]
109. Schmitz R.J., Schultz M.D., Lewsey M.G., O’Malley R.C., Urich M.A., Libiger O., Schork N.J., Ecker J.R. (2011). Transgenerational epigenetic instability is a source of novel methylation variants. Science 334: 369–373. [DOI] [PMC free article] [PubMed] [Google Scholar]
110. Secco D., Wang C., Shou H., Schultz M.D., Chiarenza S., Nussaume L., Ecker J.R., Whelan J., Lister R. (2015). Stress induced gene expression drives transient DNA methylation changes at adjacent repetitive elements. eLife 4: 09343. [DOI] [PMC free article] [PubMed] [Google Scholar]
111. Seeley K.A., Byrne D.H., Colbert J.T. (1992). Red light-independent instability of oat phytochrome mRNA in vivo. Plant Cell 4: 29–38. [DOI] [PMC free article] [PubMed] [Google Scholar]
112. Shahid S., Axtell M.J. (2014). Identification and annotation of small RNA genes using ShortStack. Methods 67: 20–27. [DOI] [PMC free article] [PubMed] [Google Scholar]
113. Shalem O., Dahan O., Levo M., Martinez M.R., Furman I., Segal E., Pilpel Y. (2008). Transient transcriptional responses to stress are generated by opposing effects of mRNA production and degradation. Mol. Syst. Biol. 4: 223. [DOI] [PMC free article] [PubMed] [Google Scholar]
114. Skirycz A., Inzé D. (2010). More from less: plant growth under limited water. Curr. Opin. Biotechnol. 21: 197–203. [DOI] [PubMed] [Google Scholar]
115. Suzuki N., Devireddy A.R., Inupakutika M.A., Baxter A., Miller G., Song L., Shulaev E., Azad R.K., Shulaev V., Mittler R. (2015). Ultra-fast alterations in mRNA levels uncover multiple players in light stress acclimation in plants. Plant J. 84: 760–772. [DOI] [PMC free article] [PubMed] [Google Scholar]
116. Thomsen S., Anders S., Janga S.C., Huber W., Alonso C.R. (2010). Genome-wide analysis of mRNA decay patterns during early Drosophila development. Genome Biol. 11: R93. [DOI] [PMC free article] [PubMed] [Google Scholar]
117. To T.K., Kim J.-M. (2014). Epigenetic regulation of gene responsiveness in Arabidopsis. Front. Plant Sci. 4: 548. [DOI] [PMC free article] [PubMed] [Google Scholar]
118. Vogel M.O., Moore M., König K., Pecher P., Alsharafa K., Lee J., Dietz K.-J. (2014). Fast retrograde signaling in response to high light involves metabolite export, MITOGEN-ACTIVATED PROTEIN KINASE6, and AP2/ERF transcription factors in Arabidopsis. Plant Cell 26: 1151–1165. [DOI] [PMC free article] [PubMed] [Google Scholar]
119. Willmann M.R., Berkowitz N.D., Gregory B.D. (2014). Improved genome-wide mapping of uncapped and cleaved transcripts in eukaryotes--GMUCT 2.0. Methods 67: 64–73. [DOI] [PubMed] [Google Scholar]
120. Wilson P.B., Estavillo G.M., Field K.J., Pornsiriwong W., Carroll A.J., Howell K.A., Woo N.S., Lake J.A., Smith S.M., Harvey Millar A., von Caemmerer S., Pogson B.J. (2009). The nucleotidase/phosphatase SAL1 is a negative regulator of drought tolerance in Arabidopsis. Plant J. 58: 299–317. [DOI] [PubMed] [Google Scholar]
121. Xu X., Xie G., He L., Zhang J., Xu X., Qian R., Liang G., Liu J.-H. (2013). Differences in oxidative stress, antioxidant systems, and microscopic analysis between regenerating callus-derived protoplasts and recalcitrant leaf mesophyll-derived protoplasts of Citrus reticulata Blanco. Plant Cell Tissue Organ Cult. 114: 161–169. [Google Scholar]
122. Yao Y., Ni Z., Peng H., Sun F., Xin M., Sunkar R., Zhu J.-K., Sun Q. (2010). Non-coding small RNAs responsive to abiotic stress in wheat (Triticum aestivum L.). Funct. Integr. Genomics 10: 187–190. [DOI] [PMC free article] [PubMed] [Google Scholar]
123. Yoo S.-D., Cho Y.-H., Sheen J. (2007). Arabidopsis mesophyll protoplasts: a versatile cell system for transient gene expression analysis. Nat. Protoc. 2: 1565–1572. [DOI] [PubMed] [Google Scholar]
124. Yu X., Willmann M.R., Anderson S.J., Gregory B.D. (2016). Genome-wide mapping of uncapped and cleaved transcripts reveals a role for the nuclear mRNA cap-binding complex in cotranslational RNA decay in Arabidopsis. Plant Cell 28: 2385–2397. [DOI] [PMC free article] [PubMed] [Google Scholar]
125. Zhai J., Arikit S., Simon S.A., Kingham B.F., Meyers B.C. (2014). Rapid construction of parallel analysis of RNA end (PARE) libraries for Illumina sequencing. Methods 67: 84–90. [DOI] [PubMed] [Google Scholar]
126. Zhang J., Mao Z., Chong K. (2013). A global profiling of uncapped mRNAs under cold stress reveals specific decay patterns and endonucleolytic cleavages in Brachypodium distachyon. Genome Biol. 14: R92. [DOI] [PMC free article] [PubMed] [Google Scholar]
127. Zhang X., et al. (2015). Plant biology. Suppression of endogenous gene silencing by bidirectional cytoplasmic RNA decay in Arabidopsis. Science 348: 120–123. [DOI] [PubMed] [Google Scholar]
128. Zhu X.-G., Ort D.R., Whitmarsh J., Long S.P. (2004). The slow reversibility of photosystem II thermal energy dissipation on transfer from high to low light may cause large losses in carbon gain by crop canopies: a theoretical analysis. J. Exp. Bot. 55: 1167–1175. [DOI] [PubMed] [Google Scholar]