Growth-limiting drought stress induces time-of-day-dependent transcriptome and physiological responses in hybrid poplar

Abstract Drought stress negatively impacts the health of long-lived trees. Understanding the genetic mechanisms that underpin response to drought stress is requisite for selecting or enhancing climate change resilience. We aimed to determine how hybrid poplars respond to prolonged and uniform exposure to drought; how responses to moderate and more severe growth-limiting drought stresses differed; and how drought responses change throughout the day. We established hybrid poplar trees (Populus × ‘Okanese’) from unrooted stem cutting with abundant soil moisture for 6 weeks. We then withheld water to establish well-watered, moderate and severe growth-limiting drought conditions. These conditions were maintained for 3 weeks during which growth was monitored. We then measured photosynthetic rates and transcriptomes of leaves that had developed during the drought treatments at two times of day. The moderate and severe drought treatments elicited distinct changes in growth and development, photosynthetic rates and global transcriptome profiles. Notably, the time of day of sampling produced the strongest effect in the transcriptome data. The moderate drought treatment elicited global transcriptome changes that were intermediate to the severe and well-watered treatments in the early evening but did not elicit a strong drought response in the morning. Stable drought conditions that are sufficient to limit plant growth elicit distinct transcriptional profiles depending on the degree of water limitation and on the time of day at which they are measured. There appears to be a limited number of genes and functional gene categories that are responsive to all of the tested drought conditions in this study emphasizing the complex nature of drought regulation in long-lived trees.


Background
Episodic drought can be catastrophic for forest trees and for the ecosystems they support (Assal et al. 2016;Brodribb et al. 2020). In response to drought, trees may slow their photosynthetic and stomatal conductance rates (Xiao et al. 2009;Hamanishi et al. 2012;Attia et al. 2015), reduce their overall growth rates (Han et al. 2019;Yi et al. 2020) and experience increased rates of crown thinning and mortality (Worrall et al. 2013). Because of their long lives, trees experience drought as a spatially and temporally heterogeneous stressor (Assal et al. 2016). For example, drought responses are interwoven with underlying daily and seasonal changes in water demands (McLaughlin et al. 2003), light availability (Aranda et al. 2005) and changes in gene expression (Filichkin et al. 2011). During periods of drought, trees develop morphologically distinct leaves with altered stomatal density (Hamanishi et al. 2012;McKown et al. 2014) and hormone profiles (Paul et al. 2018) that reflect the water status of the plant during their development and which directly influence the ability of the tree to withstand a period of drought (Himes et al. 2021). Additionally, because drought occurs over time, it is often coincident with other environmental stressors that may modulate drought response synergistically, additively or antagonistically (Mittler 2006;Plessis et al. 2015).
Creating growth-limiting soil water conditions in controlled environments like greenhouses and growth chambers that approximate drought conditions in nature is essential for studying physiologically relevant drought-stress responses (Verelst et al. 2013;Lovell et al. 2016). Experiments that have directly compared drought-stress responses in controlled and field environments have reported significant differences in physiological responses and gene expression profiles depending on the type of environments where plants were grown (Lovell et al. 2016;Wilkins et al. 2016). Some of these differences can be attributed to the degree or duration of the drought (Monroe et al. 2021), or to the underlying complexity of the environment (Mittler 2006;Plessis et al. 2015). In nature, the onset of drought is often gradual, and plants adjust their physiology and development to minimize damage and optimize growth (Sohn et al. 2013). In contrast, in controlled environments, trees are grown in relatively small volumes of soil and dry out relatively quickly (Wilkins et al. 2009). As such, a water-withholding treatment in controlled environments may be experienced as a shock rather than a drought and may elicit responses that diverge strongly from drought responses in the field (Lovell et al. 2016). Moreover, long-lived trees in nature can use their large root volumes to exploit water from deeper soil horizons thereby extending the period of water sufficiency which may allow for the development of new tissues and organs during soil drying (Mackay et al. 2020). Differences between responses in controlled and field environments are also critical considerations for genetic enhancement of stress tolerance emphasizing the need to study stress responses in conditions in which plants are grown (Zaidem et al. 2019;Groen et al. 2020).
Fluctuations in diel signals, including light, temperature and atmospheric humidity, have a substantial influence on plant water demand and drought responses. These environmental oscillations are linked to daily patterns of physiological and developmental processes, such as photosynthesis and respiration ) and leaf expansion (Pantin et al. 2011). In addition to these physiological diel oscillations, the transcriptome also varies throughout the day, even while soil water conditions remain constant (Filichkin et al. 2011). These predictable daily transcription patterns, together with oscillations in water demand, result to distinct drought transcriptomes at different times of day (Wilkins et al. 2009;Hamanishi et al. 2010;Dubois et al. 2017). Thus, to capture the full spectrum of stress responses elicited in poplars, drought-responsive alterations in physiology and gene expression must be examined at multiple times of the day.
Trees of the genus Populus are well suited for studying molecular drought responses as they are propagated clonally, have a high-quality reference genome (Tuskan et al. 2006) and vary with respect to their tolerance to drought (Barchet et al. 2014;Attia et al. 2015). Moreover, trees of this genus are foundational to many temperate ecosystems and are important as a short rotation fibre crop and for carbon sequestration throughout the Northern hemisphere (Fang et al. 2007). Understanding the genetic mechanisms underpinning these responses and which create the adaptive variation in drought tolerance in forest trees may illuminate potential avenues for enhancing drought tolerance (Estravis- Barcala et al. 2020). To gain broader insight into how poplar trees respond to drought stress, we designed an experiment to overcome some of the limitations of controlled environment drought experiments. We grew hybrid poplar trees (Populus × 'Okanese') (Schroeder et al. 2013) in greenhouse conditions for 42 days under well-watered conditions. We then transferred the trees to one of three gravimetric soil water conditions for 3 weeks. After the trees had been acclimated to the new conditions, we measured a variety of plant growth and photosynthesis parameters, in addition to gene expression profiles in fully expanded leaves to characterize elements of the drought response. The aims of this work were to determine how hybrid poplar responds to prolonged and uniform exposure to water-limiting conditions; to determine if the responses to moderate and more severe drought stresses were qualitatively or quantitatively different; and to determine how the response to drought changes throughout the day.

Plant growth conditions
Unrooted hybrid poplar stems (Populus × 'Okanese') were provided by Dr Raju Soolanakayahally from the hybrid poplar collection at the Agriculture Indian Head AgriFood Canada Research Farm (Saskatchewan, Canada). The cuttings were soaked in water for 24 h at 4 °C in the dark, trimmed to remove dead ends and planted in 4-L pots filled with Sungro Sunshine Professional Growing Mix (Sun Gro Horticulture, Agawam, MA, USA). Trees were grown in a climate-controlled greenhouse with 22-24 °C days and 15-17 °C nights. Artificial supplemental lighting was turned on at 0600 h to coincide with sunrise and a photoperiod of 15 h was maintained throughout the experiment. All trees were fertilized with Miracle-Gro All-Purpose Water-Soluble Plant Food 24-8-16 (MiracleGro, Marysville, OH, USA) once before the water-withholding treatments were initiated.

Sample collections
To ensure that all measurements were made on fully expanded leaves that developed entirely during the treatment period, the uppermost leaf on each tree was tagged 42 days after planting. On this day, all trees were assigned to one of three treatment groups with gravimetric soil water contents (SWCg) of 30 %, 50 % or 80 % of field capacity. For the subsequent 21 days, pots were weighed daily and watered as needed to maintain them at their target soil water contents [see Supporting Information- Fig. S1]. Photosynthetic measurements were made 21 days after water treatments were initiated, once between 1000 and 1200 h, and again between 1400 and 1600 h, using a Li-Cor LI-6400XT Infrared Gas Analyser (LI-COR, Lincoln, NE, USA). Leaves for transcriptome analysis were harvested at 1000 h (late morning) and 1600 h (early evening), flash-frozen in liquid nitrogen and stored at −80 °C until further processing.

Image analysis
The area of leaves that had developed and expanded during the treatments was measured using PlantCV (Gehan et al. 2017). Leaf images were captured in a lightbox using an iPhone-X camera. The images were processed using the multiplant workflow to estimate leaf area. First, the leaf images were masked from the background pixels to isolate the leaf pixels and then, the area was quantified in pixels using the shape_img function of PlantCV. Analysis of variance was carried out in R to determine statistical significance between the leaf area of the treatment groups.

Transcriptome library preparation and sequencing
Leaf tissue was ground to a fine powder under liquid nitrogen and total RNA was extracted using RNeasy Plant MiniKits (Qiagen, Toronto, ON, Canada). RNA was sent to Genome Québec (Montréal, QC, Canada) where NEBNext-stranded sequencing libraries were prepared according to manufacturer's instructions (New England Biolabs, Whitby, ON, Canada) and single-end 100 base pair reads were sequenced on an Illumina HiSeq4000. A minimum of 23 million reads were generated for each library [see Supporting Information- Table S1]. Raw and processed data are available for download as accession number GSE191155 on NCBI's Gene Expression Omnibus (Edgar et al. 2002).

RNA-seq data analysis
Processing of the raw read files was accomplished using highperformance computing clusters provided by WestGrid (www. westgrid.ca) and Compute Canada (www.computecanada. ca). Read quality was assessed using FastQC (Andrews 2010). Low-quality reads and adaptor sequences were removed with Trimmomatic (Bolger et al. 2014). The reads were aligned to the Populus trichocarpa v4.1 reference genome (Tuskan et al. 2006) (downloaded from the Phytozome database; Goodstein et al. 2012) using HISAT2 (Kim et al. 2019) and converted to BAM format with SAMtools ). Transcript abundance was quantified using featureCounts (Liao et al. 2014). All downstream analyses were performed in R (R Core Team 2021) using edgeR (Robinson et al. 2010) and plotted using ggplot2 (Wickham 2011). In total, 31 033 genes were detected from the RNA-seq assay. Genes were filtered at an average log counts per million (logCPM) > −2 and < 9 to remove low gene counts and outliers, resulting in 23 700 genes (removed 226 with avgLogCPM > 9, and 7107 with avgLogCPM < −2; see Supporting Information-Tables S2 and S3). Raw counts were normalized using the trimmed mean of M-values method (Robinson et al. 2010), and dispersion estimates were calculated using the estimateGLMRobustDisp function. Principal component analysis (PCA) was performed using the plotPCA function in DESeq2 (Love et al. 2014) on the filtered, normalized counts (n = 23 700 genes). Normalized counts were converted to logCPM, and all downstream analyses were performed on the logCPM values. Differentially expressed genes (DEGs) were identified as those with a false discovery rate (FDR)-adjusted P-value < 0.05 using the Benjamini-Hochberg method (Benjamini and Hochberg 1995). The code for the stacked bar plot (Fig. 5A) figure was modified from Bjornson et al. (2021). All code for this project is available at https://github.com/Wilkins-Lab/Poplar2021.

Clustering and Gene Ontology term enrichment analysis
Transcript profiles were clustered using cutTreeDynamic from the package dynamicTreeCut (Langfelder et al. 2008) using a minimum cluster size of 20 and a deepsplit value (sensitivity) of 1. Gene set enrichment analysis was performed using topGO (Alexa and Rahnenfuhrer 2021), using Fisher's exact test (P < 0.05). Gene Ontology terms (P. trichocarpa V3.0 GO annotation) were downloaded from agriGO v2.0 (Tian et al. 2017). Only 21 638 genes (<50 % of all genes) have annotated GO terms; as such, some functional categories may not be reported as enriched because of limitations on the GO annotation for poplar.

Results
Growth and photosynthesis rates are reduced with declining soil water content After 21 days growing in soil media with SWCg of 20-30, 40-50, 70-80 % of field capacity, trees in the 80 % SWCg group were taller ( Fig. 1A) and had developed more new leaves (Fig. 1C) than trees in the other two treatment groups. Notably, this was a graduated response, where trees in the 50 % SWCg group were taller than trees in the 30 % SWCg group in both measured aspects. The difference in total leaf area for trees between 50 % and 80 % SWCg groups ( Fig.  1E) is principally the product of the number of leaves ( Fig.   1C) as the size of newly developed leaves in these treatment groups was indistinguishable ( Fig. 1B and D). This is in contrast to trees in the 30 % SWCg group, which developed fewer (Fig. 1C) and smaller ( Fig. 1B and D) leaves than either of the other treatment groups. Additionally, the growth rate of the trees decreased upon water withholding, with the 50 % and 80 % SWCg groups stabilizing after ~14 days, while the 30 % SWCg group growth rate continued to decline throughout the experiment [see Supporting Information- Fig. S2].
To explore the physiological consequences of water limitation, we measured photosynthetic, transpiration and stomatal conductance rates in the morning (1000-1200 h) and afternoon (1400-1600 h) of the 21st day of the water-withholding treatments (Fig. 2). All three of these parameters decreased gradually with decreasing soil water availability at both times of day ( Fig. 2A-C); however, for a given treatment group, none of the parameters differed significantly between the morning and afternoon measurements (P < 0.05). Median water-use efficiency (WUE) was slightly higher in trees in the 50 % SWCg group than in trees in either of the other treatments, but this difference was not significant (Fig. 2D). Wateruse efficiency was extremely variable for trees in the 30 % SWCg treatment group, with both the highest and lowest calculated WUEs for any treatment group in the morning and afternoon. To explore the source of this variation, we plotted the relationship between the actual soil water content of each pot, which was below the upper limit of the target range at the time of the physiological measurement, and WUE [see Supporting Information- Fig. S3]. This analysis determined that WUE increased in the trees in the 30 % SWCg group, as long as the overall weight of the wet soil media was greater than ~1000 g; WUE rapidly declined if the water content was below this threshold. Because the photosynthetic parameters were measured before plants were re-watered to the upper limit of their soil water content, the photosynthesis measurements were made when plants may have been in even more water-limiting conditions.

Time of day is the strongest source of variation in the transcriptome
To understand molecular mechanisms underlying the long-term response to water deficit, we measured the transcriptomes of leaves harvested from hybrid poplar trees grown in three soil water content treatments at two times of day as described above. In total, more than 600 million 100-bp single-end reads were generated across 18 samples [see Supporting Information- Table S1]. Approximately 90 % of reads passing the quality threshold ('surviving reads') mapped to unique loci in the P. trichocarpa v4.1 reference genome (Tuskan et al. 2006) [see Supporting Information- Table S1]. To identify the main sources of variation in the transcriptome data, we conducted a PCA on normalized read counts. The first principal component (PC1) explained 94 % of the variance in the transcriptome data (Fig. 3A). The samples were cleanly divided along this axis by the time of day at which they were collected, indicating that time of day was the greatest source of variation in the transcriptome data. The second (PC2, 2 % of variation) and third component (PC3, 1 % of variation) explained the variance in response to the drought treatments ( Fig. 3A; see Supporting Information- Fig. S4). The samples are imperfectly stratified by the water treatment group along this axis. The transcriptomes of leaves from 30 % SWCg were strongly separated from the transcriptomes of the other treatment groups in the late morning; this separation was less apparent in the transcriptome data from the early evening.
Given the strong effect of time of day in the data set, we mapped the transcript abundance profiles of genes with wellcharacterized daily expression patterns to explore their expression in the present experimental conditions. We plotted the transcript abundance of 18 circadian clock genes (Fig.  3B). As expected, the morning-phase clock genes CCA1 and PRR9, for example, have higher expression values in the late morning, and other afternoon/evening genes have higher expression in the early evening (Fig. 3B). Paralogues of PRR7 and PRR9 were found to have different expression patterns and had different expression profiles in response to the drought treatments. This was consistent with divergent expression of clock gene paralogues in Brassica rapa (Greenham et al. 2020). For each water treatment, transcripts that were differentially expressed between the two sampling times ( Fig. 3C; FDR-adjusted P < 0.05) were identified. In total, 16 148 DEGs were identified of which approximately half were higher in the morning, and half were higher in the evening ( Fig. 3C; see Supporting Information-Tables S4 and S5). More than 9500 genes were differentially expressed between the sampling times in all treatment groups, while only 3809 genes were differentially expressed at a single level of drought stress (Fig. 3C). While most genes that were differentially expressed between the morning and evening samples were unaffected by the soil water conditions, several small groups of transcripts were affected by both the time of sampling and by water availability (Fig. 3D; see Supporting Information- Fig. S5).
Genes that responded to both drought and time of day were enriched for several GO terms including RNA processing and

DEGs at each time of day have distinct functions in the drought-stress response
For each time of day, we identified DEGs between the droughttreated (30 % and 50 % SWCg) and the control samples (80 % SWCg, FDR < 0.05). In total, 971 DEGs were identified between the late morning samples, and 1448 DEGs were found between the early evening [see Supporting Information- Fig.  S7]. Transcript profiles were then hierarchically clustered using the complete linkage method, resulting in 960 of the morning DEGs and 1423 of the evening DEGs assigned to clusters. Clusters that included more than 20 genes and had highly correlated (pairwise correlation > 0.99) transcript profiles were further characterized (Fig. 4A-D; see Supporting Information- Fig. S7 and Supporting Information- Table  S6). Amongst these clusters, four principal expression patterns were observed ( Fig. 4C and D): transcripts that had increased abundance only at the lowest level (30 % SWCg) of water availability (patterns I: 226 genes, and VI: 27 genes), transcripts whose abundance increased gradually with decreasing soil water availability (patterns II: 227 genes, and VII: 560 genes), transcripts that had decreased abundance only at the lowest level (30 % SWCg) of water availability (patterns III: 389 genes, and VIII: 40 genes) and transcripts whose abundance decreased gradually with decreasing water availability (patterns IV: 118 genes, and IX: 631 genes). In the early evening DEGs, however, two additional distinct patterns are seen: upregulated and downregulated in low and moderate water availability (pattern V: 66 genes, and X: 99 genes). These two patterns account for 12 % (165/1423) of the clustered early evening genes. Additionally, in the early evening, 84 % (1191/1423) of clustered DEGs had a graded Measurements were taken between 1000 and 1200 h and again between 1400 and 1600 h. Water-use efficiency was calculated as the ratio of P n to g s . Significance was determined by Tukey's Honest Significant Difference test (P < 0.05). response to water availability, wherein transcripts in the 50 % SWCg accumulated to a moderate level between the 30 % and 80 % SWCg values (pattern VI: 560 genes; pattern IX: 631 genes). In contrast, in the late morning, 64 % (615/960) of the clustered genes were increased (pattern I) or decreased (pattern III) only in the 30 % SWCg treatments.
We next performed a gene set enrichment analysis of each of the clustered gene patterns. The enriched GO terms revealed triangles for EE samples); colour corresponds to soil water treatment (blue for SWCg 80 %; green for SWCg 50 %; red for SWCg 30 %). (B) Heatmap of transcript profiles of selected clock genes. The heatmap colours correspond to the z-score for each gene in each condition. The z-score is a scaling method used to display relative expression levels. Blue tones indicate that expression is below the mean expression value and yellow tones indicate higher than the mean expression value. (C) 16 148 genes were differentially expressed between the morning and evening sampling points. Venn diagram showing that many genes were differentially expressed (|log 2 fc| > 0, FDR < 0.05) in all soil water conditions. (D) Heatmap of the 16 148 genes that were differentially expressed between the morning and afternoon time points (same genes as in B). Heatmap colours correspond to z-scores. The arrows identify clusters of genes whose expression profiles are influenced by both time of day of sampling and water status. distinct functional responses within the gene expression patterns of DEGs ( Fig. 4E and F; see Supporting Information- Fig. S8). Notably, we found that functionally different gene sets were enriched for the same expression patterns observed in the morning and evening. For example, genes that were downregulated with decreasing water availability were enriched for translational initiation and posttranslational modification in the morning, but cell signalling and DNA binding in the evening. Moreover, we show that common GO terms were enriched in different expression profiles in the morning and evening samples. For instance, functions related to cell wall modification were expressed only in low water availability in the morning, but in both low and moderate water availability in the evening.

Global analysis of time-of-day transcriptome identifies poplar response to water stress
To determine if any genes were universally responsive to the water stress treatment in our experiments, we selected the top 1000 DEGs (FDR rank-ordered) for each water stress treatment (30 % and 50 % SWCg) compared to the control (80 % SWCg) at each time of day [see Supporting Information- Table  S7]. This analysis found just over half of the possible 4000 genes that could be generated by this analysis (2362 genes) were differentially expressed in only one contrast (Fig. 5A). The remaining genes were differentially expressed in two (565 genes), three (152 genes) or all four (13 genes) contrasts ( Fig.  5A; see Supporting Information- Table S8). The list of DEGs identified in the late morning 50 % vs. 80 % SWCg contrast was dissimilar to the other three DEG lists; more than 80 % of DEGs (812/1000) on this list were unique to this treatment. Of the 565 DEGs that were common to two patterns, 435 (77 %) were common to different treatments at the same time of day.
To explore the quantitative relationships between the milder and more severe water-limiting treatments, we calculated Spearman's rank correlation coefficient for the top 1000 DEGs for each treatment as described above. We found a weaker correlation between the expression of the genes in response to milder (50 % SWCg) and more severe water deficits (30 % SWCg) in samples collected in the morning (rho = 0.65, Fig. 5B) compared to samples collected in the evening (rho = 0.81, Fig. 5C), although the slope of these relationships was similar. We also determined a moderate correlation between the expression of DEGs in the more severe water-limiting conditions (30 % SWCg) between the morning and evening (rho = 0.64, Fig. 5D), while no relationship existed between the expression of DEGs in the milder waterlimiting conditions (50 % SWCg) across sampling time (rho = 0.01, Fig. 5E). We examined the distribution of the FDRadjusted P-values for each list of top 1000 DEGs and determined that contrasts including the 30 % SWCg treatment had much lower P-values than contrasts which included the 50 % SWCg drought treatment [see Supporting Information- Fig.  S9]. This was consistent with the very small number of DEGs identified in response to the 50 % SWCg treatment.

Discussion
Plants can acclimate to water-limiting conditions during prolonged exposure In the current study, we show that reduced water availability slows stem growth in hybrid poplar and that the rate of growth stabilizes after 2-3 weeks in the tested conditions. This finding is consistent with a poplar drought study that found growth rates stabilized following a 15-day water-deficit stress. Previous studies on drought stress in poplars have also shown that radial stem growth continued to decline over 11 days of water withholding (Bloemen et al. 2016) and that the overall growth rate did not stabilize until 8-15 days after the water deficit was initiated (Bogeat-Triboulot et al. 2019). In our study, growth rate stabilization was concomitant with robust decreases in photosynthetic and stomatal conductance rates in plants grown with limited water availability. Studies in poplar (Attia et al. 2015), Kentucky bluegrass (Poa pratensis) (Saud et al. 2016) and Arabidopsis thaliana (Bechtold et al. 2016) that have measured photosynthesis and gas exchange parameters throughout longer-term soil drying experiments have shown that significant changes in these parameters occur only after plants have acclimated to their new conditions. In this study, we showed that leaves that developed and matured entirely during stress conditions were smaller and fewer than leaves that developed without water limitation. Detailed studies of the effects of drought on poplar leaf development have demonstrated that drought leads to variations in morphological and anatomical features including thicker leaf tissues, increased stomatal density (Regier et al. 2009) and reduced leaf size (Hamanishi et al. 2012;Viger et al. 2016), in addition to biochemical changes, such as total non-structural carbohydrate concentrations (Regier et al. 2009). Such differences in leaves give rise to differences in photosynthetic parameters (Cao et al. 2012) and responses to water availability (McKown et al. 2014). Temporal and spatial variation in the physiological and morphological features of leaves during drought suggests that molecular measurements made during the initial period of stress may reflect transitional states rather than adaptive responses to stressful environmental conditions.
Despite the clear effects that water limitation had on the growth and gas exchange parameters in the current study, global RNA sequencing revealed the number of genes that were differentially expressed and the magnitude of gene activity were lower than had been reported in previous transcriptome drought studies in poplar (Wilkins et al. 2009;Hamanishi et al. 2010;Raj et al. 2011). We hypothesized that this difference may exist because, unlike in the previous studies, the current study measured the transcriptomes of leaves that developed in water-limiting conditions, where the trees had time to acclimate to the soil water treatments before sampling. However, the results of our transcriptomic analysis are consistent with a study in switchgrass (Lovell et al. 2016) that identified fewer DEGs after long-term drought than in short-term dry-down experiments. These disparities may be a result of differential gene expression between short-term and long-term drought responses. This explanation is consistent with the concurrent stabilization of growth and changes in photosynthetic parameters that we observed in the current study. Although trees may experience abrupt changes in water availability in the field, longer-term water limitations due to erratic precipitation patterns are more ecologically relevant to tree growth and productivity (Brodribb et al. 2020).

Time of day and the extent of drought stress affect photosynthesis and global changes in gene expression
In the current study, as in other studies of drought response in poplar, distinct patterns of gene expression were invoked depending on the time of day of observation (Wilkins et al.  Hamanishi et al. 2010), in addition to the extent of water deficit experienced by the trees (Yan et al. 2012;Cao et al. 2014). By far, the largest source of variation in gene expression in the present study was the time of day of sampling. This was anticipated since an estimated 30-50 % of Arabidopsis genes (Bläsing et al. 2005;Covington et al. 2008), 74 % of B. rapa genes (Greenham et al. 2020) and 60 % of poplar genes (Wilkins et al. 2009;Filichkin et al. 2011) have predictable and reproducible daily changes in transcript abundance in non-stress conditions. Responses to any additional perturbations, such as environmental stressors, must be invoked on this underlying, dynamic transcriptional landscape (Greenham and McClung 2015;Bonnot and Nagel 2021). Studies in poplar (Wilkins et al. 2009;Hamanishi et al. 2010;Raj et al. 2011), Arabidopsis Dubois et al. 2017;Blair et al. 2019) and rice (Oryza sativa) (Wilkins et al. 2016;Grinevich et al. 2019;Desai et al. 2021) have shown that transcriptional responses to stress are different at different times of day. This was consistent with our findings here that the majority of drought-responsive genes had distinct expression profiles in the morning and early evening. Our observations may be attributable to underlying diurnal variation in expression patterns. For example, global changes in gene expression in response to water limitation or a given level of water availability are perceived as a stressor only at some times of day (McLaughlin et al. 2003).
Here we show that different functional classes of droughtresponsive genes were enriched at different times of day and with different degrees of water limitation. For example, genes involved in the regulation of iron transport were induced by water limitation in the morning only. The transport of iron from the vasculature to the chloroplasts is regulated by the circadian clock (Salomé et al. 2013), and is affected by water availability (Pan et al. 2015). We also showed that genes involved in protein folding and translation initiation were reduced in water-limited samples in the morning but were unchanged in response to drought later in the day. This finding may be consistent with the growing literature showing that stress responses may be regulated in part by translation (Reynoso et al. 2019;Bonnot and Nagel 2021) and posttranslational modifications (Rodriguez et al. 2021). Together, these findings highlight the importance of including multiple time points in experimental design for capturing the full suite of stress response mechanisms.
In contrast to the strong time-of-day effect on gene expression profiles, we determined that gas exchange parameters, including carbon assimilation and transpiration rates, were most strongly affected by the extent of water limitation. We showed that in more severe water-limiting conditions, WUE was highly variable and included high and low WUE values. As previously described in maize (Blankenagel et al. 2018), in the present study there appeared to be a water content threshold below which WUE declined precipitously. This observation may be related to exponential decreases in soil water potential at low SWCg (Dowd et al. 2019), whereby a small change in water volume leads to a large change in water availability. Moreover, physiological and growth traits, including specific leaf area and leaf carbon content, are affected non-linearly by soil water availability (Monroe et al. 2021) highlighting the importance of studying these responses under a limited number of discrete soil moisture conditions. To fully characterize plant responses to drought stress, it will be necessary to incorporate intermediate water-limiting conditions to understand how the degree of water stress impacts plant status, rather than an 'all-or-nothing' approach.
Response to drought stress relies on qualitative and quantitative changes in the transcriptome Response to drought stress in poplar is not limited to a subset of high-amplitude drought-responsive genes; rather, it is a transcriptome-wide adjustment that is affected by the time of day and by the extent of the stress applied. In this study, we show that drought-induced transcriptome profiles across different water levels and time of day lead to qualitative and quantitative changes, as seen in previous studies (Lovell et al. 2016). Qualitative changes refer to the specific genes invoked, as the poplars expressed different genes depending on their water status and time of day, while quantitative changes are represented by the degree of change observed, as the water status induced variable proportions of genes. We find that genes differentially expressed in more severe water-limiting conditions were also differentially expressed, though not significantly so, in the moderate water-limiting conditions at both times of day, and that this relationship was stronger in the early evening than in the late morning. Notably, low water availability in the morning and evening elicited similar degrees of differential expression across a wide panel of genes, but there was no correlation between DEGs in moderate water availability at the two times of day.
These findings demonstrate that water availability and diel signals interact to determine what conditions will be perceived as stress and which biological mechanisms will be invoked in response to environmental perturbations. Previous studies in Arabidopsis (Legnaioli et al. 2009) and soybean (Marcolino-Gomes et al. 2014) have shown the interconnectivity of the circadian clock and the drought response, likely due to the reciprocal regulation of abscisic acid with clock components (Legnaioli et al. 2009;Pokhilko et al. 2013). Widespread changes in the timing, magnitude and diversity of gene regulation are necessary for plants to adequately respond to stress (Greenham and McClung 2015) and so experimental studies must take this complexity into account.

Supporting Information
The following additional information is available in the online version of this article- Figure S1. Daily weights of pots after initiation of waterdeficit treatments. Once pots had reached their target weights, pots were re-watered to the upper limit of the target weights each morning. Figure S2. The growth rate of poplars stabilizes at 14 days after initiation of water deficit, except in 30 % gravimetric soil water content (SWCg) conditions. Figure S3. Water-use efficiency (WUE) of three gravimetric soil water content (SWCg) conditions on the last day of water deficit shows variability in the 30 % SWCg group. Figure S4. Principal component analysis of filtered normalized counts for all detected transcripts (23 700 genes). Shape of the data point corresponds to time of sampling (diamonds for late morning (LM) samples; triangles for early evening (EE) samples); colour corresponds to soil water treatment (blue for gravimetric soil water content (SWCg) 80 %; green for SWCg 50 %; red for SWCg 30 %). PC1 vs. PC2 is shown in Fig. 3A. Figure S5. Hierarchical clustering and expression patterns of genes that were differentially expressed between the late morning (LM) and early evening (EE) (gene lists in Supporting Information-Tables S4 and S5). The coloured bar below the dendrogram represents the clusters identified by WGCNA. The scaled expression profile of the identified clusters is presented as a line graph with mean values shown as a red line. Expression profiles are ordered according to the clusters (coloured bars) from left to right. Figure S6. Selected enriched Gene Ontology terms for gene clusters that responded to water deficit and time of day. Genes in each cluster are listed in Supporting Information- Table S5. Figure S7. Hierarchical clustering and scaled gene expression profiles of differentially expressed genes (DEGs) in the (A) late morning (LM) and (B) early evening (EE). Expression profiles are ordered according to the coloured bars from left to right. Mean expression of the cluster is shown with a red line. Figure S8. Full list of Gene Ontology (GO) terms for the (A) late morning and (B) early evening gene expression profiles. See Fig. 4 for list of patterns (I-X). Figure S9. False discovery rate (FDR)-adjusted (Benjamini-Hochberg (BH)) P-values of the top 1000 differentially expressed genes (DEGs) in each drought treatment (30 % and 50 % gravimetric soil water content (SWCg)) compared to the well-watered treatment (80 % SWCg) from Fig. 5A. Table S1. Processing of raw RNA-sequencing reads. Surviving and dropped reads are from Trimmomatic, the overall alignment rate is from HISAT2 and overall gene assignment rate is from featureCounts. Table S2. Raw expression counts of all the genes passing quality filtering step (average logCPM > −2 and < 9). Table S3. Scaled expression counts (logCPM) of all the genes passing quality filtering step (average logCPM > −2 and < 9). Table S4. Raw expression count of all differentially expressed genes between the late morning (LM) and early evening (EE). Table S5. Scaled expression (logCPM) of all differentially expressed genes between the late morning (LM) and early evening (EE) including cluster assignments. Table S6. Clusters of late morning (LM) differentially expressed genes (DEGs) (LM30vsLM80 & LM50vsLM80; n = 960) and early evening (EE) DEGs (EE30vsEE80 & EE50vsEE80; n = 1423). Table S7. Top 1000 differentially expressed genes (false discovery rate (FDR)-ranked) from each water stress treatment compared to control. Table S8. Log 2 fold-changes of shared top 1000 differentially expressed genes (false discovery rate (FDR)-ranked) between comparisons of time of day and extent of water stress.