Time-Series Transcriptomics Reveals That AGAMOUS-LIKE22 Affects Primary Metabolism and Developmental Processes in Drought-Stressed Arabidopsis

In Arabidopsis thaliana , changes in metabolism and gene expression drive increased drought tolerance and initiate diverse drought avoidance and escape responses. To address regulatory processes that link these responses, we set out to identify genes that govern early responses to drought. To do this, a high-resolution time series transcriptomics data set was produced, coupled with detailed physiological and metabolic analyses of plants subjected to a slow transition from well-watered to drought conditions. A total of 1815 drought-responsive differentially expressed genes were identi ﬁ ed. The early changes in gene expression coincided with a drop in carbon assimilation, and only in the late stages with an increase in foliar abscisic acid content. To identify gene regulatory networks (GRNs) mediating the transition between the early and late stages of drought, we used Bayesian network modeling of differentially expressed transcription factor (TF) genes. This approach identi ﬁ ed AGAMOUS - LIKE22 ( AGL22 ), as key hub gene in a TF GRN. It has previously been shown that AGL22 is involved in the transition from vegetative state to ﬂ owering but here we show that AGL22 expression in ﬂ uences steady state photosynthetic rates and lifetime water use. This suggests that AGL22 uniquely regulates a transcriptional network during drought stress, linking changes in primary metabolism and the initiation of stress responses.


INTRODUCTION
Water limitation in agriculture is poised to intensify in the coming decades due to urbanization, industrialization, depletion of aquifers, and increasingly erratic rainfall patterns exacerbated by climate change (Easterling et al., 2000;Christensen et al., 2007;Seager et al., 2007;Famiglietti and Rodell, 2013). Reduced water availability leads to drought stress, which is a major constraint on the physiology, growth, development, and productivity of plants (Boyer, 1970(Boyer, , 1982Lobell and Field, 2007;Roberts and Schlenker, 2009;Skirycz et al., 2010;Lobell et al., 2011;Verelst et al., 2013). Therefore, understanding the mechanisms of drought response in plants is essential for the improvement of plant performance under water-limiting conditions and has been the subject of many investigations over the years Yamaguchi-Shinozaki, 1997, 2007;Chaves et al., 2009;Nakashima et al., 2009;Pinheiro and Chaves, 2011). Water deficit responses are complex and require stress sensing and signaling to adjust plant growth, maintain water status through osmoregulation, prevent water loss through decreases in stomatal conductance, and activate detoxification processes (Passioura, 1996;Chaves et al., 2003;Pinheiro and Chaves, 2011). An important consideration is that even a slight reduction in water availability can elicit stomatal closure and a reduction in CO 2 assimilation and in combination with the diversion of resources toward drought defense mechanisms will affect plant productivity (Chaves et al., 2003).
Plants have adopted different strategies to respond to water limitation, such as drought escape through early flowering and reducing the size of plants to increase water use efficiency or drought avoidance through enhanced soil moisture capture or reduced transpiration (Ludlow, 1989;Blum, 2005;Aguirrezabal et al., 2006;Franks, 2011). In this context, the influence of drought on plant development and growth through its effects on developmental processes such as germination, seedling growth, and leaf development has been studied extensively in the past decade (van der Weele et al., 2000;Finkelstein et al., 2002;Xiong et al., 2006;Yaish et al., 2011). Optimal timing of flowering and inflorescence development are important traits essential in determining plant yield, and these can vary greatly in response to water limitation (Eckhart et al., 2004;Franke et al., 2006;Su et al., 2013;Ma et al., 2014).
At the cellular level, plants respond to drought with changes in gene expression and protein and metabolite abundances (Charlton et al., 2008;Harb et al., 2010;Wilkins et al., 2010;Baerenfaller et al., 2012), which are part of defense mechanisms and detoxification processes (Shinozaki and Yamaguchi-Shinozaki, 2007;Begcy et al., 2011;Ozfidan et al., 2012). Recent progress in genomics, transcriptomics, and bioinformatics has paved the way for dissecting drought-response mechanisms and has enabled the targeted manipulation of drought-responsive genes in plants. For example, the overexpression of a number of genes that code for transcription factors (TFs) leads to drought resistance (Sakuma et al., 2006;Nelson et al., 2007;Chen et al., 2008;Quan et al., 2010;Tang et al., 2012).
In many studies to identify genes important in the regulation of drought responses, the effects of water limitation at the transcriptional level have been analyzed by exposing plants to severe dehydration. This involves treatments such as cutting and air drying leaves and/or roots or induction of osmotic shock through the application of highly concentrated osmotica such as polyethylene glycol or mannitol (Kreps et al., 2002;Seki et al., 2002;Kawaguchi et al., 2004;Kilian et al., 2007;Weston et al., 2008;Fujita et al., 2009;Abdeen et al., 2010;Deyholos, 2010;Mizoguchi et al., 2010). These experiments have substantially increased our knowledge of molecular responses under severe drought stress, but they do not always reflect physiological conditions experienced by drought-stressed soil-grown plants (Bechtold et al., 2010(Bechtold et al., , 2013Harb et al., 2010;Wilkins et al., 2010;Lawlor, 2013;Zhang et al., 2014). Physiological responses such as stomatal conductance, photosynthetic performance, and metabolic changes are usually not measured during the progression of the drought stress, and the varied nature of the stress induction treatments makes comparative analysis between experiments problematic. Slow developing soil water deficits have different physiological consequences than those induced by rapid tissue dehydration and therefore possibly utilize different gene networks (Chaves et al., 2003(Chaves et al., , 2009Pinheiro and Chaves, 2011).
This inconsistency among experiments was first noted in a metaanalysis of microarray experiments comparing air drying, soil drying, and mannitol treatments (Bray, 2004). This analysis found very few differentially expressed genes (DEGs) common to all treatments (Bray, 2004). Consequently, recent experiments have focused on soil-grown plants (Harb et al., 2010;Wilkins et al., 2010;Zhang et al., 2014). From these studies, an overall integrative picture of the temporal responses to drought is emerging slowly, and it is clear that use of a single or a small number of time points and different types of experimental conditions lead to very different outcomes. One consequence of this is that very little is known about the early events in the perception of drought stress signals (Ueguchi et al., 2001;Wohlbach et al., 2008;Pinheiro and Chaves, 2011).
To address the above issues, we set out to gain detailed information on the processes that occur during the transition from wellwatered to drought conditions, in which the intensity of the stress becomes gradually greater. We monitored the physiological and metabolic status of plants through a progressive drought experiment and mapped onto these data the temporal responses of the transcriptome. Our intention was to use the highly resolved transcriptional profiling data to construct gene regulatory networks (GRNs) using dynamic Bayesian network modeling (Beal et al., 2005;Breeze et al., 2011;Penfold and Wild, 2011) with the aim of identifying regulatory genes functional during drought perception and signaling. The goal was to link early physiological and metabolic drought avoidance responses with later drought escape and/or tolerance responses (Claeys and Inzé, 2013). This initially required testing of the network modeling to evaluate the capability of these approaches to identify genes important in the regulation of drought responses. This was achieved by selecting a highly connected candidate gene, AGAMOUS-LIKE22 (AGL22; also known as SHORT VEGETATIVE PROTEIN), from the GRNs. AGL22 has an established function in plant development (Gregis et al., 2013;Méndez-Vigo et al., 2013), but in this study, it was shown to play a thus far undiscovered role in the critical early stages of the plant's response to drought. These results demonstrated the potential value of experimental strategies that combine time-series transcriptomics data with dynamic modeling as a means of identifying stress-responsive genes.

RESULTS
Time-series experiments were performed analyzing physiological, metabolic, and transcriptional changes in Arabidopsis thaliana to reveal the chronology of plant responses to drought stress. A progressive slow-drying experiment starting at 95% relative gravimetric soil water content (rSWC) and drying down to 17% rSWC was performed on 5-week-old Arabidopsis plants ( Figure 1). To determine the severity of the stress, daily measurements of relative leaf water content (RWC; Figure 1A) and leaf water potential were also performed ( Figure  1B). During the experiment, the average rSWC loss was ;10% per day, but RWC was maintained throughout the progressive drying period until the point of wilting at 17% rSWC ( Figure 1A).

Maximum Photosynthetic Capacity Responds Similarly in Well-Watered and Drought-Stressed Plants during a Progressive Drought Experiment
Stomatal conductance (g s ) and photosynthetic carbon assimilation (A) were measured daily on well-watered and drought-stressed plants through the progressive drought treatment ( Figures 1C and  1D). Stomatal conductance declined at ;60% rSWC (day 5; Figure  1C), which was followed by a decline in carbon assimilation at ;45% rSWC (day 7), indicating that stomatal diffusional limitations affected carbon assimilation ( Figure 1D). Plant growth evaluated as rosette fresh weight and rosette area ceased at ;40% rSWC (Supplemental Figures 1A to 1C).
The light and CO 2 -saturated maximum photosynthetic rate (A max ), maximum rate of carboxylation (V Cmax ), and the rate of ribulose-1,5bisphosphate (RuBP) regeneration (J max ) showed no difference between well-watered and drought-stressed plants ( Figure 2A). In addition, maximum and operating efficiencies of photosystemII (Fv/Fm, Fv9/Fm9, and Fq9/Fm9;Baker, 2008) showed no change during the drought period ( Figure 2B; Supplemental Figure 1D), suggesting that the overall primary metabolic capacity was maintained as drought conditions progressed. The sequential changes in photosynthetic physiology, relative water content, and leaf water potential suggest that these conditions allowed us to capture the transition between early physiological changes and later stress responses.

Metabolite Profiling Indicates the Stable Nature of Primary Metabolism during Drought Stress
The decline in stomatal conductance and carbon assimilation led us to perform metabolite analysis to evaluate changes in primary and secondary metabolism (Figure 3). Untargeted liquid chromatographymass spectrometry metabolite profiling was performed on samples harvested at early (day 2, ;80% rSWC), mid (day 7, ;45% rSWC), and late (day 13, 17% rSWC) stages of the drought stress. This analysis showed that the majority of the metabolome was unchanged throughout most of the drought treatment, and distinct clustering between well-watered and drought-stressed samples emerged only by the final day of drought stress (17% rSWC) ( Figure  3A). Leaf development was a major factor for sample separation, with days clustered more closely together than treatments ( Figure 3A).
Targeted metabolite analysis was performed to determine the foliar levels of 102 stress-associated compounds (Supplemental Data Sets 1 to 3), which also revealed a mainly late response for many of these stress-associated metabolites (Supplemental Data Sets 1 to 3). Often these changes were limited to the last two to three time points (between ;30 and 17% rSWC; Supplemental  Figures 2 and 3). For example, metabolites indicative of drought stress increased only during the late stages of the dehydration period ( Figures 3B and 3C). There was a significant increase in abscisic acid (ABA) levels during the last four time points (Xiong et al., 2002; Figure 3B), while proline, a drought stress-responsive compatible solute in vascular plants (Sperdouli and Moustakas, 2012), accumulated to significant levels only during the last two time points ( Figure 3B). Additionally, the accumulation of secondary metabolites commonly associated with stress responses such as anthocyanins and flavonols were altered during the late stages of the drought response (Sperdouli and Moustakas, 2012; Supplemental Data Sets 1 and 3 and Supplemental Figure 3). Oligosaccharides/disaccharides associated with osmotic protection during drought and osmotic stresses (galactinol and raffinose; Taji et al., 2002) significantly accumulated during the last 4 days of the drought response ( Figure 3C; Supplemental Data Set 1). In conclusion, leaf metabolism remained largely stable during the first 9 d of the experiment, while changes previously associated with drought stress (Taji et al., 2002;Xiong et al., 2002;Sperdouli and Moustakas, 2012) only became evident during the last three to four time points.

Transcriptomics Analysis on a Single Leaf Identifies 1815 DEGs during Progressive Drought Stress
Transcriptome profiling was performed on leaf 7 to integrate the complex physiological and metabolic responses with changes at the gene expression level. Leaf 7 was fully expanded at the time of the experiment (Supplemental Figure 1C) and was chosen because a detailed temporal transcriptome analysis of leaf development was available (Breeze et al., 2011). Single leaf 7 samples for transcriptome analysis were taken each day at the midpoint of the light period. RNA from four leaf samples per treatment and time point was hybridized on CATMA v4 arrays (Sclep et al., 2007;see Methods). An adapted MAANOVA (microarray analysis of variance) method was used to analyze the data for each comparison (Wu et al., 2003;Churchill, 2004;Breeze et al., 2011;Windram et al., 2012). This generated a single normalized expression value for each gene. A Gaussian process two-sample test (GP2S; Stegle et al., 2010) was used to identify DEGs. Choosing a Bayes factor value (likelihood of differential expression) of >6 resulted in a total of 1815 DEGs (Supplemental Data Set 4). The upregulated group of genes showed on overrepresentation of Gene Ontology (GO) terms related to carbohydrate biosynthesis, flavonoid, and secondary metabolic processes, while downregulated genes were enriched in protein translation, cell wall-associated processes, pigment biosynthesis, and chloroplast associated processes (Supplemental Data Set 5).
GO terms related to stress, dehydration, and hormonal regulation, including ABA, were not enriched in the complete data set. This result suggested that the overall progressive drought experiment was not a severe dehydration stress response, as indicated by the maintenance of primary metabolic capacity ( Figure  2) as well as the late responses of stress-associated metabolites ( Figure 3). Leaf water potential has been used as a measure of the progression and effect of drought stress on plants (Zhang et al., 2014). We estimated the cumulative number of DEGs at each time point by determining the time of first differential expression for each gene (Supplemental Data Set 6; Figure 4A). In our experiment, leaf water potential correlated significantly with rSWC ( Figure 4B) as well as the number of DEGs ( Figure 4C) and showed a weaker correlation with carbon assimilation ( Figure 4D). The biggest drop in leaf water potential occurred between 40% (day 8) and ;30% (day 9) rSWC, which coincided with the biggest increase in DEGs ( Figures 4A and 4C), potentially indicating a shift from mild to severe drought stress.
To evaluate the transcriptome data set in the context of other drought experiments, we also compared our data set to two soilbased drought studies in Arabidopsis. The first experiment performed by Harb et al. (2010) comprised a microarray comparison of leaf samples under moderate drought stress (maintaining soil water content at 30% of field capacity) and progressive drought stress at the prewilting (;15% field capacity) and wilting (;10% field capacity) stages. In the progressive drought stress treatment, 3005 genes responded >2-and <0.5-fold, in comparison to 441 genes for the moderate drought treatments (Supplemental Data Set 7). The second study analyzed samples at a loss of 25% soil water of field capacity measured at 6-h intervals across a 24-h period (Wilkins et al., 2010) and identified 570 genes that responded across a 24-h interval (hereafter called diel; Supplemental Data Set 7). A general overview of overlapping genes between the different experiments showed that only 30 genes were common to all 4 treatments ( Figure  5A). Among the overlapping genes were known stress-responsive, ABA-responsive, and secondary metabolism genes, which predominantly responded during the latter half of the drought (C) Relative concentrations of galactinol (black bars) and raffinose (gray bars). The data represent the mean of the ratio (n = 4; 6SE; *P < 0.05 or **P < 0.01). experiment (Supplemental Data Set 8). While there were common elements in all four treatments, 63% of the DEGs in our data set were unique to the time series ( Figure 5A) and were potentially related to the adjustments to early and moderate drought stress.

Slow Soil Drying Induces a Senescence Response in Leaf 7
To assess developmental changes in leaf 7 during drought stress, we compared the drought time series to an Arabidopsis leaf 7 senescence time-series data set (Breeze et al., 2011;Supplemental Data Set 7). In all, 842 genes overlapped with the senescence data set (p-hyper, 4.4E-135; Figure 5B), of which 83% responded between days 8 and 13 (40 to 17% rSWC; Supplemental Data Set 9). The overlap contained genes associated with oxidation/reduction-related processes, pigment biosynthesis, and primary metabolism, which were predominantly downregulated during drought stress ( Figures 5C and 5D; Supplemental Data Set 10). The induction of senescence-related processes during drought stress is a known phenomenon (Munné-Bosch and Alegre, 2004) and was further confirmed by a significant overlap of genes between the published drought and senescence data sets (Supplemental Data Set 7 and Supplemental Figure 4A).
In addition, changes in the expression of secondary metabolism genes were observed in the drought time series ( Figure 5E), including CHALCONE SYNTHASE, FLAVONOL SYNTHASE1, LEUCOANTHOCYANIDIN DIOXYGENASE, PRODUCTION OF ANTHOCYANIN PIGMENT1, and ANTHOCYANINLESS2 (Supplemental Figure 5 and Supplemental Data Set 9). This coincided with increased accumulation of flavonol and anthocyanin (from day 11; Supplemental Data Set 1 and Supplemental Figure 3) and suggested that the plants had entered a severe stress phase (Vanderauwera et al., 2005). Therefore, it was concluded that slow soil drying induces senescence in leaf 7 but only at the point of severe drought stress. By contrast, early responses to soil drying (days 1 to 7) were mostly unique to the drought time series. Importantly for later considerations, the induction of leaf senescence in response to drought did not affect flowering time (Supplemental Figure 4B).

Temporal Clustering Reveals Coregulated Groups of Genes, but Does Not Reveal Specific Regulatory Mechanisms
The cumulative number of DEGs at each time point (Supplemental Data Set 6) confirmed that major gene expression changes occurred late during the drought experiment as the number of genes that showed first differential expression at each time point was highest between days 8 and 11 ( Figure 4A). A total of 336 genes responded during the first half of the experiment, while the majority of genes (1479) showed first differential expression during the latter half of the experiment. A Euclidean distance matrix of the average expression values of the four biological replicates for each data set was generated and used in hierarchical cluster analysis. The resulting dendrogram showed that samples clustered in relation to treatments and within the drought treatment into early and late stage responses ( Figure 6A). Dividing the data set into early (days 1 to 7, ;95 to 45% rSWC) and late responses (days 8 to 13, ;40 to 17% rSWC) also revealed functional groups of genes responding at different times throughout the progressive drought. Early upregulated genes were associated with carbohydrate and glycoside biosynthetic processes, general carbohydrate metabolic processes, and inorganic cation transporter activities (Table 1). The late upregulated genes encompassed flavonoid and secondary metabolite biosynthesis, while the late downregulated genes were involved in translation, pigment biosynthesis, photosynthesisrelated processes, and oxidation/reduction processes (Table 1). To gain insight into the molecular factors underlying this temporal (A) Comparative meta-analysis of the 1815 DEGs with publicly available drought data sets. The Venn diagram shows the overlap of time series DEGs with those responsive to moderate (mDr; Harb et al., 2010) or progressive drought (pDr; Harb et al., 2010) and a moderate drought at different times of day (diel; Wilkins et al., 2010). (B) Comparative meta-analysis of the 1815 DEGs with a publicly available leaf 7 senescence time-series data set (Breeze et al., 2011). The Venn diagram shows the overlap of drought and senescence DEGs. (C) Overview of antioxidant, photosynthesis. and photorespiration-related gene expression at two different time points (95% rSWC and 17% rSWC). (D) Overview of oxidation/reduction-related gene expression at two different time points (95% rSWC and 17% rSWC). (E) Overview of secondary metabolism-related gene expression at two different time points (95% rSWC and 17% rSWC). All MapMan diagrams show gene expression data in leaf 7, where blue indicates increased and yellow indicates decreased gene expression according to the scale. Each square represents a single gene within the pathways. separation of samples, hierarchical cluster analysis of the DEGs was performed using SplineCluster (Heard et al., 2005) on the basis of gene expression patterns in the drought stressed leaf only. Using a prior precision value of 0.01, the 1815 genes were divided into 28 clusters ( Figures 6B and 7; Supplemental Data Set 11). The first 14 clusters showed an overall upregulation and contained 1149 genes, while the last 14 downregulated clusters contained 667 genes ( Figure 7; Supplemental Data Set 11). The 28 clusters were hypothesized to represent groups of genes that are coregulated during the drought experiment. To explore potential regulatory mechanisms of genes clustered in specific temporal expression profiles, we analyzed each individual cluster for over-and underrepresentation of GO terms in the Biological Process and Molecular Function and for overrepresentation of known TF binding motifs in promoters (Supplemental Data Sets 12 and 13 and Supplemental Figure 6). A few clusters showed enrichment of expected GO terms in response to drought, such as flavonoid biosynthesis, photosynthesis, pigment biosynthesis, and response to stress (Supplemental Data Set 12; Figure 7). However, most clusters did not show any enrichment or underrepresentation of GO terms (Supplemental Data Set 12; Figure  7), and only two clusters (cluster 1 and 9) contained the ABRE binding motif, known to perceive ABA-mediated drought and osmotic stress signals (Supplemental Figure 6 and Supplemental Data Set 13; Kim et al., 2011).

Bayesian State Space Modeling Identifies Genes That Link Drought Responses to Plant Development through Regulating Transcriptional Networks
The data presented so far did not allow us to draw conclusions about a specific regulatory mechanism across the whole time series but suggested that early and late responses to a decline in soil water content are regulated differently. Large-scale transcriptional reprogramming and metabolic adjustment did not play a dominant role during the early phases of the dehydration response. Nevertheless, among the 337 DEGs responding between 95 and 40% rSWC were 33 TF genes (Supplemental Data Set 14), of which 25% had a functional annotation of development (GO:0032502; Supplemental Data Set 14). This suggested that a reprogramming of developmental processes during drought stress may have occurred in response to the observed early physiological changes (closure of stomata and reduction in leaf water potential). We reasoned that the expression of early TF genes must therefore play a role in orchestrating this acclimation to drought stress.
Metropolis Variational Bayesian State Space Modeling (M-VBSSM; Penfold, University of Warwick; Supplemental Methods) was initially performed on 176 differentially expressed TFs in the data set (Supplemental Data Set 14). This approach selects subsets of TFs to generate a network, which is continually updated by probabilistic replacement of TFs to generate a series of networks and provides a consensus model based on the marginal likelihood (Supplemental Methods). From the consensus model, we could calculate the occurrence of each TF (the number of times a particular TF appeared over all models) and a count of the number of downstream connections each TF had across all models at a particular z-score, in this case indicating a 95% confidence threshold (Supplemental Data Set 15). TFs that scored well in both rankings were deemed highly connected hubs with many significant edges. The M-VBSSM consensus model indicated that developmental genes played an important role in the regulation of drought, as four out of the top 10 highly connected TFs were associated with the regulation of plant development (Supplemental Data Set 15). In addition, two of the top 10 TFs were among the group of early responding TFs (Supplemental Data Sets 14 and 15). While the consensus model was useful for the initial ranking of genes, it did not represent a causal model, instead representing a type of averaging of many different network models. For this reason, we opted to model a smaller selection of genes, which was advantageous for two reasons. First, the final model was smaller and sparser, and therefore more interpretable, and second, the resultant interactions could be interpreted causally.
Therefore, the top 10 "hub" TFs with the highest frequency of occurrence and 90 random transcription factors were chosen for analysis with the VBSSM package (Beal et al., 2005; Supplemental Data Set 16). As part of this selection, we also chose early and late responding TFs. In total, 19 early responding TFs (days 1 to 7) were included in the model to establish a potential transcriptional link between early and late responses (Supplemental Data Set 16).
The resultant model placed the early-responding TF gene AGL22 at the center of a 25-TF gene network ( Figure 8A). AGL22 was differentially expressed in the drought experiment beginning at day 5 (Supplemental Figure 7A). Fifteen TF genes that were part of the GRN were initially analyzed by qPCR to check for differential expression under drought in wild-type plants (20% rSWC). All genes were differentially expressed in line with the levels observed in the microarray experiment, except for PACLOBUTRAZOL RESISTANCE1, BASIC HELIX-LOOP-HELIX038, and AUXIN RESPONSE FACTOR1 ( Figure 8B; Supplemental Figure 7B).
AGL22 is known to affect flowering time and plant development (Supplemental Figure 7C; Gregis et al., 2013;Méndez-Vigo et al., 2013); however, it was not regulated during leaf senescence (Supplemental Data Set 9), suggesting that AGL22 uniquely regulated a transcriptional network during drought stress. This was further explored by performing VBSSM using the control timeseries data set of the same transcription factor genes. If AGL22 was a hub gene in the control time series, it would suggest a role in developmental reprogramming over the 13-d experimental period regardless of drought stress. The Bayesian modeling resulted in a number of fragmented connections of a small number of genes (Supplemental Figure 8A), suggesting these genes were not part of a gene regulatory network under well-watered conditions. The highly connected genes in the drought model, AGL22 and RAP2.12, did not feature, not even as peripheral genes (Supplemental Figure 8A). Therefore, the early responding TF AGL22 was chosen for further analysis to establish how far an unbiased modeling approach can be used to identify genes capable of influencing plant drought phenotypes and downstream network connections. Two independent T-DNA insertion lines were isolated (see Methods), both of which were confirmed knockout mutants for AGL22 ( Figure 8C; Supplemental Figure 8B). The mutants were subsequently analyzed for their effect on the AGL22-centered network interactions after drought stress. Eight out of 15 TF genes differentially expressed under drought conditions exhibited altered gene expression in at least one of the agl22 mutants compared with the wild type ( Figure 8D). This implied that ;50% of the network connections were regulated at least partially through AGL22, but also suggested the possibility of redundancy within the network ( Figures 8A and 8D). Four late TFs (WRKY20, GIS, DREB1A, and FBH3) were substantially downregulated in both agl22 mutants after drought, suggesting that these TFs were primarily regulated through AGL22 ( Figure 8D).

Both agl22 Mutants Had Early-Flowering, Fast-Drying Drought Escape Phenotypes
It is important to note that due to the early-flowering phenotype of the agl22 mutants (Supplemental Figure 7C), drought stress was begun at day 22 after sowing, when there was no visible differences in rosette leaf number (Supplemental Figures 9A to 9C), but with a significant increase in rosette area ( Figure 9A).
We observed an increased drying rate in both agl22 mutants, suggesting increased water use ( Figure 9B). To determine if this was due to developmental or metabolic changes, we performed light response curve measurements of photosynthesis at specific times throughout the drying period (Supplemental Figure 9D). Due to the small leaf and rosette size, light curves were measured in whole plant chambers at 90, 74, and 25% rSWC (see Methods). The increased water loss was primarily driven by a greater rosette area ( Figure 9A), despite a significant reduction in stomatal conductance ( Figure 9C). Accordingly, light saturated carbon assimilation (A sat ) was significantly reduced throughout the drying period already under wellwatered conditions ( Figure 9D), leading to a significant reduction in total aboveground biomass (Supplemental Figure 10A). Flowering time remained constant between well-watered and drought treatments in both agl22 mutants, indicating that drought stress conditions did not affect flowering time in the agl22 mutants ( Figure 9E).   Zhang et al., 2014). The slow steady drought experiment performed in this study allowed us to investigate the full range of temporal physiological, transcriptional, and metabolic responses in a single fully expanded Arabidopsis leaf. Additionally, by measuring leaf water potential ( Figure 1B) and RWC ( Figure 1A), we were able to monitor the progression and degree of drought stress in relation to the physiological, transcriptome, and metabolome changes. The decline in carbon assimilation at ;45% rSWC was primarily driven by reduced stomatal conductance limiting CO 2 diffusion. There were no underlying metabolic constraints, as photosynthetic capacity was unaffected by the drought treatment ( Figure 2; Supplemental Figure  1D), and metabolite profiles remained unchanged throughout the majority of the drying period ( Figure 3A; Supplemental Figures 2 and 3), suggesting that Arabidopsis Col-0 is a drought-tolerant ecotype.
Stomatal limitation as the primary factor in reducing photosynthesis under mild drought conditions has been observed in other studies; however, severe dehydration stress is believed to lead to metabolic constraints, associated with RuBP availability (Flexas and Medrano, 2002). We did not observe a reduction in the maximum capacity for carbon assimilation (A max ), Rubisco carboxylation (V Cmax ), and RuBP regeneration (J max ; Figure 2A). In addition, maximum and operating efficiencies of photosystem II photochemistry ( Figure 2B), photochemical and nonphotochemical quenching (Supplemental Figure 1D), were maintained throughout the drying period despite a decline in carbon assimilation ( Figures 1D), indicating very little stress on photosystem II. This suggested that an alternative electron sink, most likely photorespiration (reviewed in Chaves et al., 2003;Lawson et al., 2014), must have been operating under drought stress, and increased gene expression in the photorespiratory pathway supports this notion ( Figure 5C).
In general, two distinct phases in response to progressive soil drying could be discerned (Figures 6 and 7). Early responses were predominantly adjustments to stomatal conductance leading to restricted CO 2 diffusion for photosynthetic carbon assimilation (Boyer, 1970;Passioura, 1996; Figure 1C) with some associated transcriptional changes accounting for 17% of the 1815 DEGs (Supplemental Data Sets 1 to 12; Table 1). By contrast, late responses (from 40% rSWC) encompassed hormonal (ABA), transcriptional, and major metabolic changes associated with senescence (Supplemental Data Sets 1 to 12; Table 1). These later responses corresponded with the many different phenological and physiological changes observed in other studies, including impaired photosynthesis, increased solute accumulation, and growth arrest (Boyer, 1970;Passioura, 1996). At the cellular level, soluble sugars, oligosaccharides, antioxidants, and proline accumulation are known to enhance the tolerance to drought stress by acting as osmolytes or as reactive oxygen species scavengers, especially hydroxyl radicals (Smirnoff and Cumbes, 1989;Cuin and Shabala, 2007). The accumulation of these compounds during the latter stages of the drought period ( Figures 3B and 3C), together with the increase in secondary metabolites, such as flavonoids (Supplemental Data Set 1 and Supplemental Figures 3  and 4), suggests a role in the defense against severe drought stress (Tattini et al., 2004;Lei et al., 2006;Xiao et al., 2007;Xu et al., 2008;Harb et al., 2010;Fini et al., 2011;Page et al., 2012). In conclusion, this time series covers all phases during a progressive drought stress and therefore provides the opportunity to study different stages of stress responses in greater detail than has been previously possible (Supplemental Figure 10B).

Transcriptional Regulation of Drought Stress Responses
At the gene expression level, a slowly developing soil water deficit is different from rapid tissue dehydration. From this study and in  (Benjamini and Hochberg, 1995). "Category" indicates the timing and direction of change in gene expression across the time series. "Fold" indicates the fold enrichment.
comparison with two other drought experiments (Harb et al., 2010;Wilkins et al., 2010), it is clear that different genes respond depending on the nature of the drought stress applied ( Figure 5A). Only few genes overlapped, with the majority of genes, responding in the final 4 d of the experiment (Supplemental Data Set 8).
Previous studies have often focused on the identification of genes coding for TF classes responding to terminal or severe drought stress, including BASIC LEUCINE ZIPPER (bZIPs, e.g., ABA-responsive element binding protein/ABRE binding factor), AP2/EREBP (e.g., DREB/CBF), NAC transcription factors (NAM, The blue line indicates the mean expression profile for each of the 28 clusters. Individual genes present in each cluster are available in Supplemental Data Set 11. The red line indicates the switch from early (95 to 45% rSW; days 1 to 7) to late (40 to 17% rSWC; days 8 to 13). Selected enriched GO terms (Supplemental Data Set 12) are indicated on each cluster. ATAF1-2, CUP-SHAPED COTYLEDON2), CCAAT binding (e.g., NUCLEAR FACTOR Y), and ZINC-FINGER (e.g., C2H2 zinc finger protein) families (Umezawa et al., 2004;Bartels and Sukar, 2005;Karaba et al., 2007;Li et al., 2008;Licausi et al., 2010;Jensen et al., 2013). The majority of TF genes in our study also responded relatively late in the drought period (from 40% rSWC), especially those associated with ABA, dehydration, and oxidative stress responses (Supplemental Data Set 14). Similar classes of TF genes have also been shown to respond to drought stress in Medicago truncatula where 8% of the responding genes coding for TFs responded late throughout the drying period (Zhang et al., 2014). By contrast, among the early-responding TF genes (95 to 45% rSWC; Supplemental Data Set 14), ;25% were linked to plant development, indicating that early physiological changes (A) Gene regulatory network generated using VBSSM with the drought time-series data (threshold z-score = 1.65). The nodes highlighted in red were upregulated during drought stress including the central hub gene, AGL22. Nodes highlighted in green were genes downregulated during drought stress. Blue nodes signify genes that were not regulated by AGL22 as predicted from the model (see [D]). All red, green, and blue nodes were selected for evaluation after drought stress; (B) Relative gene expression of selected genes under drought conditions (17% rSWC). Gene expression was analyzed by qPCR. The numbers are expressed as fold changes of drought over control (n = 5 6 SE). Significance of the fold changes are indicated by either *P < 0.05 or **P < 0.01. For gene and primer list, see Supplemental Data Set 16 and Supplemental Table 1. (C) Relative expression levels of AGL22 in two knockout lines, agl22-3 and agl22-4, compared with the wild type determined by qPCR. Significance of the fold changes are indicated **P < 0.01. (D) Relative gene expression profiles of 16 genes predicted to be regulated by AGL22 under drought stress in agl22-3 (black bars) and agl22-4 (gray bars) compared with the wild type. The data represent the mean (n = 7; 6SE), and significance of the fold changes are indicated by either *P < 0.05 or **P < 0.01. Asterisks located centrally indicate both mutants are significantly different to the wild type, while asterisks located over one mutant indicate significance for the specific mutant. Due to the early-flowering phenotype of both agl22 mutant alleles, drought stress was begun at 22 d after sowing. (A) Rosette area (cm 2 ) of Col-0 (light gray), agl22-3 (black), and agl22-4 (dark gray) plants at different soil water contents (n = 5). The asterisk indicates significant difference compared with the wild type at P < 0.05. (B) Rate of water loss in agl22-3 and agl22-4 plants compared with the wild type averaged over 13 d of water withdrawal (n = 10). The asterisk indicates significant difference compared with the wild type at P < 0.05. (C) Stomatal conductance (Gs) at different soil water contents in Co-0 (light gray), agl22-3 (black), and agl22-4 (dark gray) (n = 5). The asterisk indicates significant difference compared with the wild type at P < 0.05. (D) Light-saturated carbon assimilation (Asat) at different soil water contents in Col-0 (light gray), agl22-3 (black), and agl22-4 (dark gray; n = 5). The asterisk indicates significant difference compared with the wild type at P < 0.05. (E) Days to flowering in well-watered (light gray) and drought-stressed (black) plants (n = 10). Plants were grown under short-day conditions as described in Methods. At 5 weeks, plants were subjected to progressive drought stress. When 17% rSWC was reached, plants were rewatered and flowering time was recorded as days after sowing. Control plants were maintained well watered. The asterisk indicates significant difference compared with the wild type at P < 0.05. may influence lifetime traits before the initiation of acute stress defense and senescence responses, highlighting the balancing act between the need to grow and to induce effective stress tolerance mechanisms (Claeys and Inzé, 2013).

Dynamic Bayesian Network Modeling Identifies Genes That Regulate Plant Development
Analysis of promoter binding sites (Supplemental Data Set 13 and Supplemental Figure 6) did not indicate specific regulatory networks or mechanisms during the early events. We therefore assessed the use of a high-throughput gene expression approach coupled with dynamic Bayesian network modeling to identify genes associated with the regulation of early drought responses. To make sense of large high-throughput data sets, network inference algorithms were developed, which are capable of establishing regulatory interactions among genes (Bansal et al., 2007). A number of different inference algorithms were used to successfully reconstruct known gene regulatory networks to validate these approaches (Cantone et al., 2009;Penfold and Wild, 2011). VBSSM is such an algorithm developed specifically for highly resolved temporal gene expression data sets, with the aim of identifying genes that are the key regulators in a given system (Beal et al., 2005). A recent comparison of modeling algorithms used to infer GRNs has shown that VBSSM is competitive with network reconstructions based on experimental data (Penfold and Wild, 2011;Windram et al., 2014;Penfold and Buchanan-Wollaston, 2014). Due to the limited number of experimental observations compared with the much greater number of differentially expressed genes, the system is inherently underdetermined (Penfold and Buchanan-Wollaston, 2014), and previous experience suggested that for these kinds of data sets VBSSM can model around 100 genes (Breeze et al., 2011). However, this type of approach can introduce a bias during the process of gene selection, while the M-VBSSM approach (see Methods) generally avoids this gene selection bias (Supplemental Methods). We therefore opted to use both M-VBSSM and VBSSM to identify key drought-regulatory genes and drought phenotypes associated with those genes.
The initial emergence of several development-associated transcription factors from the M-VBSSM approach (Supplemental Data Set 15) and the subsequent selection of developmental and nondevelopmental TFs for VBSSM (Supplemental Data Set 16) confirmed the flowering time regulator AGL22 as a hub gene in the drought response. Importantly, we did not observe any involvement of AGL22 in leaf senescence (Supplemental Data Set 9 and Supplemental Figure 8A), suggesting that the gene plays a unique role in drought stress responses.
The connection between flowering time and drought in Arabidopsis has been established independently in a number of studies, in which carbon isotope discrimination (Farquhar et al., 1982(Farquhar et al., , 1989 quantitative trait loci colocated with known developmental/flowering time loci Masle et al., 2005;Juenger et al., 2005;McKay et al., 2008). The identification of a flowering time gene and subsequent verification of its influence on plant water use, photosynthesis, and phenology (discussed below) suggest that this type of dynamic modeling can provide an important means of discovering genes that will produce phenotypes associated with lifetime water use and plant development. However, unlike quantitative trait locus mapping, VBSSM also managed to establish some valid network interactions from time-series transcriptomics data, potentially allowing for a temporal reconstruction of events and biological processes occurring during progressive drought stress.

A Flowering Time Gene Influences Water Use and Photosynthesis under Well-Watered and Drought Stress Conditions
Both agl22 mutants exhibited elevated water loss and rapid development already under nonstress conditions (Figures 9A and  9B). This could imply a trade-off between drought avoidance and escape in environments where drought shortens the growing season (Franks, 2011). Selecting for early flowering may be beneficial for plant survival but not necessarily for achieving high biomass (Supplemental Figure 10A), which suggests that drought survival and the ability to maintain biomass under sustained waterlimiting conditions depend on different mechanisms (Skirycz et al., 2011).
The agl22 mutants exhibited 36 and 46% reductions in the steady state light-saturated photosynthetic rate under wellwatered conditions ( Figure 9D), which appeared to be partly associated with reduced stomatal conductance ( Figure 9C). However, during drought stress, the photosynthetic rate in both agl22 mutants was reduced by only 11 and 13%, suggesting that both agl22 mutants were able to maintain substantial photosynthetic rates ( Figure 9D). This is supported by the fact that agl22 mutants also maintained rosette growth throughout the drying period in comparison to wild-type plants ( Figure 9A), and although total aboveground biomass was significantly reduced in both mutant alleles (Supplemental Figure 10A), biomass distribution shifted from vegetative growth to reproductive growth (Supplemental Figure 10C). The complex links between plant growth, primary metabolism, and flowering time in Arabidopsis are highlighted in a recent article where increased plant growth was positively associated with early-flowering phenotypes (El-Lithy et al., 2010), which may explain the larger rosette area observed prior to flowering in 30-d-old agl22 mutant plants compared with the wild type ( Figure 9A; Supplemental Figure 9B). In addition, starch/ carbohydrate status and metabolite levels have been linked to rosette growth Sulpice et al., 2009), as well as development and flowering time (Zhou et al., 1998;Moore et al., 2003;Funck et al., 2012). Both photosynthesis and flowering are regulated by the light environment and are clearly linked via the carbohydrate status (Zhou et al., 1998;Moore et al., 2003;Funck et al., 2012), connecting primary metabolism with plant growth and development. However, little is known about the link between photosynthetic performance and flowering time especially in flowering time mutants, and at this point, it is unclear why the agl22 mutants exhibited a substantial reduction in photosynthesis already under well-watered conditions. Furthermore, the observations regarding AGL22 reinforces that water use in relation to overall plant productivity requires a balance of developmental and physiological processes to successfully complete a lifecycle in the prevailing climatic conditions. AGL22 was identified from moderately drought-stressed plants, which also suggests that not all drought-responsive genes may work in all water deficit scenarios. This may especially be the case for those genes that have been selected under terminal or severe drought conditions (Hu et al., 2006;Nelson et al., 2007;Xiao et al., 2009). It is important to note that none of the TF genes selected as key regulators of the drought response in earlier studies (Hu et al., 2006;Nelson et al., 2007;Xiao et al., 2009) was a hub in the TF GRN ( Figure 8A), although many of these TF genes were differentially expressed during drought stress (Supplemental Data Set 14). This is supported by the notion that many genes identified with a role in stress tolerance under severe stress conditions seem to have little effect on plant growth in mild drought conditions (Skirycz et al., 2011).
Interestingly, two targets of AGL22, DREB1A and FBH3 ( Figure  8D), have previously been shown to be involved in the regulation of abiotic stress responses. Overexpression of DREB1A leads to drought, salt, and freezing tolerance (Kasuga et al., 1999(Kasuga et al., , 2004, while FBH3 has been shown to regulate stomatal opening (Takahashi et al., 2013) and functions in ABA signaling in response to osmotic stress (Yoshida et al., 2015). Both genes were lateresponding targets with few network connections ( Figure 8A; Supplemental Data Set 15). The model therefore may allow us to predict the role of AGL22 during drought stress and provide a potential link between mild/moderate and severe drought responses.
This study demonstrates that network inference incorporating highly resolved time-series transcriptomics data is able to predict TF networks and identify genes with regulatory importance during drought stress. Moreover, by focusing on the transition from early physiological changes to drought stress responses, we were able to identify AGL22 as a gene associated with lifetime water use. Consequently, VBSSM as a gene discovery tool promotes the selection of unknown, yet highly connected genes for further phenotypic evaluation.

Plant Material, Plant Growth, and Drought Stress
Arabidopsis thaliana plants (Col-0, agl22-3 [SALK_141674], and agl22-4 [SAIL_583_C08]) were obtained from the European Arabidopsis Stock Centre and were grown under a 8:16-h light:dark cycle at 23°C, 60% relative humidity, and light intensity of 150 µmol m 22 s 21 , using a mixture of cold and warm white fluorescent tubes. Arabidopsis seed was stratified for 3 d in 0.1% agarose at 4°C before individual seeds were sown onto a soil mix (Scotts Levington's F2+S compost:fine grade vermiculite in a ratio of 6:1). For the drought time-course experiment, pots (7 3 7 3 9 cm) were filled with the same amount of soil mix. Control pots, to determine 100 and 0% soil water content, were set up at the same time. Plants were transferred into individual pots 2 weeks after the sowing date and were kept well watered until the beginning of the drying episode at 5 weeks after sowing. Half the plants were maintained under well-watered conditions, while for the remaining half, water was withdrawn and pot weight was determined daily. Relative soil water content was calculated for each day and pots were left to dry until 17% rSWC was reached. Five-week-old plants were saturated in water to reach 95% rSWC, and watering was stopped in the treatment plants until ;17% rSWC was reached. The control plants were maintained under well-watered conditions at ;95% rSWC. Due to the early-flowering phenotype of both agl22 mutant alleles, drought experiments were performed on 22-d-old plants to ensure the experiments were performed at similar rosette developmental stages and prior to the onset of flowering, as indicated in Supplemental Figures 9B and 9C. The drying rate was determined as the slope of the decline in relative soil water content, measured daily throughout the drying period.

RNA Extractions, Labeling, Microarray Hybridization, and Analysis
Total RNA was extracted, labeled, and hybridized to CATMA v4 arrays (Sclep et al., 2007) as previously described (Breeze et al., 2011;Windram et al., 2012). The experimental design for the drought time-series hybridization is shown in Supplemental Figure 11.
Arrays were hybridized and washed as described (Windram et al., 2012). Arrays were scanned on a 428 Affymetrix scanner at wavelengths of 532 nm for Cy3 and 635 nm for Cy5. Cy3 and Cy5 scans for each slide were combined and processed in ImaGene version 8.0 (BioDiscovery) to extract raw intensity and background corrected data values for each spot on the array. The data have been deposited in the Gene Expression Omnibus under accession number GSE65046.
An adaptation of the MAANOVA package (Wu et al., 2003) was used to analyze the extracted microarray data as described by Breeze et al. (2011) and Windram et al. (2012), using a mixed-model analysis. The MAANOVA fitted model considered dye and array slide as random variables, and time point, treatment, and biological replicate as fixed variables. The model allowed assessment of the main effect of treatment, the main effect of time point, the interaction between these factors, and the nested effect of biological replicate. Predicted means were calculated for each gene in each of the 112 combinations of treatment, time point, and biological replicate and for each of the 28 combinations of treatment and time point, averaged across biological replicates.

Differential Gene Expression Analysis
Genes were ranked based on their Gaussian process two-sample (GP2S) Bayes factor differentially expressed score, a cutoff of $6 gave 2496 differentially expressed genes. Genes identified in the F-test as being differentially expressed that were not in the GP2S list of 2496 genes were added manually. The expression profiles of the genes ranked 1800 to 3150 were then plotted and assigned visually as differentially expressed or not. This resulted in a false positive rate of 23.3% for this group, and a final cutoff of 6 for the GP2S was chosen, duplicates were removed, and a list of 1934 differentially expressed genes was produced. Removal of probes and genes with no annotation in TAIR9 left a list of 1815 unique differentially expressed genes. The time at which genes first became differentially expressed (TOFDE) was subsequently determined using the GP2S timelocal method (Stegle et al., 2010).

Promoter Analysis
Publically available position-specific scoring matrices (PSSMs) were collected from the PLACE and JASPAR databases (Higo et al., 1999;Sandelin et al., 2004). PSSMs were clustered by similarity, and a representative of each cluster was chosen for screening. Promoter regions corresponding to 200 bp upstream of the transcription start site were retrieved from the Ensembl Plants sequence database (release 50). For any given PSSM and promoter, we scanned the sequence and computed a matrix similarity score (Kel et al., 2003) at each position on both strands. P values for each score were computed from a score distribution obtained by applying the PSSM to randomly generated sequences. We took the top k nonoverlapping hits and performed the binomial test (pbinom function in R Stats package) for the occurrence of k sites with observed P values within a sequence of length 200 bp. The parameter k is optimized within the range 1 to 5 for minimum binomial P value to allow detection of binding sites without a fixed threshold per binding site. To determine the presence or absence of a PSSM in a promoter, in each case the promoters were sorted by binomial P value, and we applied a cutoff to select the top 2000. For each PSSM, its frequency in promoters of each cluster was compared with its occurrence in all promoters in the genome. Motif enrichment was calculated using the hypergeometric distribution (see statistical analysis). For motif enrichment analyses P values # 1e-5 were considered significant, to allow for multiple testing.

RT-PCR and qPCR
Leaves were harvested and frozen in liquid nitrogen. Total RNA was extracted from a pool of three plants using Tri-reagent (Sigma-Aldrich) according to the manufacturer's instructions. A minimum of five replicates of whole plants from separate experiments was performed for mutant analysis and drought treatments. For cDNA synthesis for real-time qPCR and RT-PCR, 1 µg total RNA was treated with RNase-free DNase (Ambion) according to the manufacturer's instructions and reverse transcribed as previously described (Ball et al., 2004). qPCR-PCR was performed using a SYBR green fluorescence-based assay as described previously (Bechtold et al., 2010(Bechtold et al., , 2013. Gene-specific cDNA amounts were calculated from threshold cycle (Ct) values and expressed relative to controls and normalized with respect to ACTIN and CYCLOPHILIN cDNA according to Gruber et al. (2001). RT-PCR was performed to amplify the full-length AGL22 gene on both agl22 mutant alleles. The primers used for qPCR and RT-PCR are given in Supplemental Table 1.

Relative Water Content
Whole rosettes of five plants were harvested each day throughout the drying period. The RWC of the leaf was calculated using the formula: rLWC (%) = (FW 2 DW)/(SW 2 DW) 3 100, where FW is the actual rosette weight at the day of harvest, SW is the fully saturated rosette weight, and DW is the dry weight of the rosette.

Leaf Water Potential
The leaf water potential was measured via the Scholander pressure bomb technique (Scholander et al., 1964) using a SKPM1400 plant moisture system (Skye Instruments). Leaf water potential was measured daily throughout the drying period on both control and drought-stressed plants according to the manufacturer's instructions.

Plant Development and Biomass Measurements
Arabidopsis development was assessed using the scale developed by Boyes et al. (2001). Once the final flower had opened, watering was ceased and plants were bagged and left to dry out before harvesting. At harvesting, rosettes, stalks, and seeds were separated. The seed weight and dry weight of rosettes and stalks/pods were determined (Bechtold et al., 2010). At least 10 plants per line and watering regime were measured.

Photosynthetic Rate (Snapshot Measurements)
Instantaneous measurements of net CO 2 uptake rate (A) and stomatal conductance to water (gs) were made on leaf 7, using an open gas exchange system (PP Systems). Leaves were placed in the cuvette at ambient CO 2 concentration (C a ) of 400 µmol mol 21 , leaf temperature was maintained at 22 6 2°C, vapor pressure deficit was ;1 kPa, and irradiance was set to growth conditions (150 µmol m 22 s 21 ). A reading was recorded every 3 min when the IRGA conditions had stabilized (;1.5 min), but before the leaf had a response to the new environment (Parsons et al., 1997).

A/C i Curves (Maximum Photosynthetic Rates)
Five weeks after emergence, (A) and (gs) were measured on leaf 7, using an infrared gas exchange system (PP Systems). The response of A to changes in the intercellular CO 2 concentration (C i ) was measured under a saturating PFD, provided by a combination of red and white LEDs (PP Systems). In addition, the response of A to changes in PFD from saturating to subsaturating levels was measured using the same light source at the current atmospheric CO 2 concentration (390 µmol mol 21 ). All gas analysis was made at a leaf temperature of 20 (61)°C and a vapor pressure deficit of 1 (60.2) kPa. Plants were sampled between 1 and 4 h after the beginning of the photoperiod. For each leaf, steady state rates of A and g s at current atmospheric [CO 2 ] were recorded at the beginning of each measurement. The A/C i parameters, V Cmax (maximum RuBP-saturated rate of carboxylation in vivo), A max (light and CO 2 saturated rate of carbon assimilation in vivo), and J max (maximum in vivo rate of electron transport contributing to RuBP regeneration) were calculated by fitting equations described by Farquhar et al. (1980) with subsequent modifications described by McMurtrie and Wang (1993).

Light Response Curves Using Whole-Plant Chambers
A/Q response curves were measured using whole-plant gas exchange system developed at the University of Essex, with a heliospectra LED light source (Heliospectra). Input air was maintained at a relative humidity of 50 to 60%, air temperature of 22°C, and CO 2 concentration of 400 mmol mol 21 , matching that of the growth conditions. Plants were initially stabilized for 30 min at saturating irradiance 800 µmol m 22 s 21 , after which PPFD was reduced in nine steps (Supplemental Figure 9D), with assimilation (A) and stomatal conductance (g s ) being recorded at each new PPFD level.

Metabolite and Hormone Analysis
During the experimental period, two leaves per plant were harvested every day, for a total of six plants per treatment. Samples were frozen in liquid nitrogen, freeze-dried overnight, and stored at room temperature in darkness until extraction. Primary metabolites were extracted from frozen tissue with chloroform-methanol as described by Lunn et al. (2006). T6P, other phosphorylated intermediates, and organic acids were measured by highperformance anion-exchange chromatography coupled to tandem mass spectrometry as described by Lunn et al. (2006). Trehalose was measured enzymatically with fluorometric detection as described by Carillo et al. (2013).
For sugars, amino acids, hormones, and secondary metabolites, freezedried leaf powder (10 mg) was extracted in 0.8 mL methanol containing 1% acetic acid. After centrifugation (10 min at 16,100g, 4°C), the samples were filtered through a 0.2-µm PVDF syringe filter (Chromacol). For nontargeted Liquid chromatography quadrupole time-of-flight mass spectrometry metabolite profiling, 5 mL extract was injected onto a Zorbax StableBond C18 1.8 mm, 2.1 3 100 mm (QToF) reversed-phase analytical column (Agilent Technologies). Chromatography and mass spectrometry conditions were described by Page et al. (2012). Peaks were extracted and aligned using XCMS (Smith et al., 2006), and statistical analysis and data visualization were performed with MetaboAnalyst 2.0 (Xia et al., 2015). For liquid chromatography-tandem mass spectrometry analysis of hormones, 10 mg freeze-dried leaf powder was extracted in 0.8 mL 10% methanol + 1% acetic acid containing deuterated standards (Forcat et al., 2008). Secondary metabolites and hormones were analyzed with an Agilent 6420B triple quadrupole (QQQ) mass spectrometer (Agilent Technologies) coupled to a 1200 series Rapid Resolution HPLC system. Two microliters of sample extract was loaded onto a Zorbax Eclipse Plus C18 (3.5 mm, 2.1 3 150 mm) for amino acids, sugars, and hormones, respectively. The following gradient was used: 0 min to 0% B; 1 min to 0% B; 5 min to 20% B; 20 min to 100% B; 25 min to 100% B; 27 min to 0% B; 7 min post-time. QQQ source conditions were as follows: gas temperature 350°C, drying gas flow rate 9 liters min 21 , nebulizer pressure 35 psig, and capillary voltage 4 kV. The polarity, fragmentor voltage, and collision energies were optimized for each compound. The multiple reaction monitoring modes used for compound identification are shown in Supplemental Data Set 17, and data are reported as peak areas. Flavonoid identification was based on previous tandem mass spectrometry identification of flavonoids in Arabidopsis (Tohge et al., 2005;Stobiecki et al., 2006). For sugar and polyol analysis, 5 mL of sample extract was loaded onto an XBridge amide HILIC column (particlesize3.5 µm, 2.1 mm i.d. 3150mm; Waters)with a constant flow rate of 0.3 mL min 21 and a column temperature of 35°C for the duration. Mobile phases comprised water:acetonitrile with 0.1% ammonia (mobile phase A was 90% acetonitrile, and B was 10% acetonitrile with 5 mM ammonium formate). Sugars (5 µL) were separated using the following gradient: 0 to 17 min, 0 to 54% B; 17 to 19 min, 54% B; 19 to 20 min, 54 to 0% B, with a 10 min reequilibration time. The QQQ was operated in negative ion mode. Electrospray ionization source conditions were as follows: gas temperature, 350°C; drying gas flow rate, 9 liters min 21 ; nebulizer pressure, 35 psig; and capillary voltage, 4 kV. Data were acquired in selected ion monitoring mode with a dwell time of 50 ms. The fragmentor voltage was 50 V for all sugars. The sugars were quantified by reference to standards (Supplemental Data Set 17). Amino acids were separated with a ZIC-HILIC column (150 3 2.1 cm, 3.5-µm particle size; Merck SeQuant). Sample (2 µL) was injected into the column with a flow rate of 0.25 mL min 21 . Mobile phases comprised of water:acetonitrile with 0.1% formic acid (mobile phase A was 95% acetonitrile and B was 5% acetonitrile with 5 mM ammonium acetate). Compounds were separated using the following gradient: 0 to 10 min, 5 to 50% B, 10 to 15 min, 50-90% B, 15 to 20 min, 90%B, 20 to 25 min, 90 to 5% B, with an 11-min reequilibration time. The QQQ was operated in positive ion mode and electrospray ionization source conditions were as follows: gas temperature, 350°C; drying gas flow rate, 9 liters min 21 ; nebulizer pressure, 35 psig; and capillary voltage, 4 kV. Data were acquired in multiple reaction monitoring mode with a dwell time of 50 ms. Amino acids were quantified by multiple reaction monitoring (Supplemental Data Set 17) and data are reported as peak areas (Supplemental Data Set 18).

Analysis of Publicly Available Microarray Data Sets
Publicly available microarray data sets from different experiments in which Arabidopsis was subjected to drought and senescence (Harb et al., 2010;Wilkins et al., 2010;Breeze et al., 2011) were located in supplementary files of already published articles and were compared with the drought timeseries data sets using VENNY (Oliveros, 2007). Hypergeometric distributions were calculated for different overlaps using the phyper function in R version 3.0.2. A cutoff of -p(log) of 5 was chosen as highly a significant overlap between two or more data sets.

Analysis of GO
GO annotation analysis was performed using DAVID version 6.7 (Huang et al., 2009), BINGO (Maere et al., 2005), and Agrigo (Du et al., 2010) with the GO_Biological_Process category, as described by Ashburner et al. (2000). Overrepresented GO_Biological_Process and GO_Molecular_Function categories were identified using a hypergeometric test with a significance threshold of 0.05 after Benjamini-Hochberg correction (Benjamini and Hochberg, 1995) or Bonferroni correction (Holm, 1979) with the whole annotated genome as the reference set.

VBSSM and M-VBSSM
Significant numbers of genes can be differentially expressed in response to environmental stress, which, given the limited number of experimental measurements, means that network models are often unidentifiable (Penfold and Buchanan-Wollaston,2014;Windram et al., 2014, and references therein). Furthermore, the interpretation of large, densely connected networks can often be difficult, and any hypothesis we extract from them can therefore be ambiguous. One solution is to select a more limited number of genes to model, either based upon prior knowledge, heuristic approaches, or random selection. The reduced number of genes means network inference approaches can be applied such as the VBSSM of Beal et al. (2005). Within the VBSSM the expression of the genes can be written in the form: where the term ½BC þ D ij captures all information about how gene j regulates gene i. Rather than infer a point estimate for each interaction ½BC þ D ij , Beal et al. (2005) infer a posterior distribution and use standard Z-statistics to assess the statistical significance. However, the preselection of genes described above may result in some bias. Here, we chose to additionally build a network model around a particular gene of interest, using random selection via a Metropolis algorithm. At each step in the Metropolis algorithm, a Bayesian state space model is fitted to the timeseries gene expression profiles for the selected genes, and the marginal likelihood or "model evidence" used as the selection criteria. In this way, we can infer small network models around each gene that we are interested in (for full details, see Supplemental Methods and Supplemental Figure 12).
For the drought data, a total of 176 transcription factors were differentially expressed (Supplemental Data Set 14). We therefore used the Metropolis model selection to systematically build a network of 88 genes around each of the 176 genes in turn. Within the Metropolis selection, each of the 176 network models was run for 2000 iterations in the MCMC chain, by which point the marginal likelihood was seen to be plateau and the algorithm was terminated. The 176 networks at step 2000 were then combined to create a meta-network, which was used to compile summary statistics, such as the number of times a particular gene was found in each of those 176 network models or the number of downstream connections a particular gene had over those 176 models. Ranked lists of genes can be found in Supplemental Data Set 13. Of the top 10 genes with the highest number of downstream connections, four were annotated as being developmental in nature, suggesting a link between drought response and developmental programs. A selection of 99 random differentially expressed transcription factors including the top 10 highly connected genes were selected from the larger pool of 176 TFs that were DE during the drought stress (Supplemental Data Set 16). Both the control and drought time series for these genes were normalized to have zero-mean and unit variance and subsequently modeled using the VBSSM of Beal et al. (2005) using a z-score of 1.65 to select the control and drought-specific networks, and a final VBSSM model (Beal et al., 2005) was fitted to the gene expression for AGL22.

Statistical Analysis
Statistical analyses were performed using SPSS version 19.0. Parameter differences between the wild type and agl22 mutants were determined using one-way ANOVA with appropriate post-hoc analysis. Tukey's HSD test was used if variances of means were homogenous, and Games Howell test, if variances were not homogenous. The SE of the calculated ratios of fold differences for metabolite and gene expression data and errors of individual means were combined "in quadrature" as the final ratio was a combination of the error of the two different means of the control and drought stress samples. Correlations were estimated among drying rate and flowering time as the standard Pearson product-moment correlation between the genotype means. Hypergeometric distributions were analyzed using the phyper function in the R stats package. Generally, P values # 1e-5 were considered significant to allow for multiple testing. The metabolite data were analyzed by between subjects two-way ANOVA for time series data and probability values with false discovery rate multiple testing correction are tabulated. For secondary compounds, amino acids, and sugars, compounds not detected in >50% of samples were discarded and remaining missing data imputed using KNN (Supplemental Data Set 3). Abundance data were normalized to total signal in each sample, log 2 transformed, mean-centered, and divided by SD of each variable using Metabolonalyst 3.0 (Xia et al., 2015).

Accession Numbers
Sequence data from this article can be found in the Gene Expression Omnibus data library under accession number GSE65046. Supplemental Figure 8. Validation of the knockout phenotype in agl22 insertion mutants and the specific role of AGL22 during drought stress.
Supplemental Figure 9. Growth and photosynthetic phenotype of Col-0 and agl22 mutants during drought stress.
Supplemental Figure 10. Biomass production in agl22 mutants compared with the wild type and timeline of events.
Supplemental Figure 11. Schematic overview of array hybridizations across the 13 time points and two different treatments.
Supplemental Figure 12. Validation of the M-VBSSM approach.
Supplemental Table 1. Primers for qPCR, mutant screen, and RT-PCR analyses. Supplemental Data Set 18. Output from XCMS alignment of peaks from LC-QToF metabolite profiling.