Changes in Gene Expression in Space and Time Orchestrate Environmentally Mediated Shaping of Root Architecture[CC-BY]

Tracking cell-type transcriptomic responses to treatments that change root development at an unprecedented level of temporal detail reveals molecular mechanisms that shape root architecture. Shaping of root architecture is a quintessential developmental response that involves the concerted action of many different cell types, is highly dynamic, and underpins root plasticity. To determine to what extent the environmental regulation of lateral root development is a product of cell-type preferential activities, we tracked transcriptomic responses to two different treatments that both change root development in Arabidopsis thaliana at an unprecedented level of temporal detail. We found that individual transcripts are expressed with a very high degree of temporal and spatial specificity, yet biological processes are commonly regulated, in a mechanism we term response nonredundancy. Using causative gene network inference to compare the genes regulated in different cell types and during responses to nitrogen and a biotic interaction, we found that common transcriptional modules often regulate the same gene families but control different individual members of these families, specific to response and cell type. This reinforces that the activity of a gene cannot be defined simply as molecular function; rather, it is a consequence of spatial location, expression timing, and environmental responsiveness.


INTRODUCTION
As sessile organisms, plants are highly responsive to the environment, at both the molecular and phenotypic levels. A facet of multicellularity is that it enables organisms to have the scope for differentiation between different cell types and within a cell type as it develops. To maintain compartmentalization of specialized functions, different cell types maintain distinct transcriptional programs. This discovery was facilitated by the use of fluorescence-activated cell sorting (FACS), which enables isolation and analysis of endogenous responses in specialized cell types (Birnbaum et al., 2003Brady et al., 2007). Cell-type specialization manifests not only in the development of different cell type functions, but also in how different cell types respond to their environment. For example, lateral root organogenesis is initiated specifically in the pericycle cell type through a program of asymmetric, anticlinal cell division and is a root developmental process that is highly responsive to environmental change (Malamy and Benfey, 1997). This enables the architecture of the root to be tuned to particular environments or remodeled over time, a manifestation of plant plasticity. In response to nitrogen influx, a vast array of changes occurs at the molecular level in plants, and these changes influence root development. However, it is unclear how short-term and long-term changes are linked temporally. Furthermore, although many environmental factors affect root development, we do not know the extent to which lateral root development regulation in different environments is controlled by the same genes.
FACS-based approaches to cell-type-specific profiling have previously been used to study a range of processes and environmental responses in the root. Previous studies have addressed the responses of distinct root layers to stresses, including salinity (Dinneny et al., 2008;Geng et al., 2013), low pH, and sulfur (Iyer-Pascuzzi et al., 2011). FACS has also been used to quantify endogenous auxin levels at the tissue-specific level (Petersson et al., 2009) and to identify auxin-responsive genes (Bargmann et al., 2013). Other studies have employed FACS to isolate cells of the root quiescent center (Nawy et al., 2005) and generate metabolic profiles of distinct root cell types (Moussaieff et al., 2013).
Complementary to spatial analysis of responses is the dissection at a temporal level. Collection of sequential time points of transcriptomic data enables the use of network inference methods to reconstruct causative gene networks Windram et al., 2012;Lewis et al., 2015). Using network inference methods, it is thus possible to make predictions of regulatory interactions and how these interactions are influenced by different treatments or under different conditions (Hecker et al., 2009;Krouk et al., 2010;Emmert-Streib et al., 2012). Furthermore, networks can identify control hubs and enable prediction of the effect of genetic perturbations . In the case of root environmental responses, nitrogen responses have been studied at the whole root level over the first minutes of nitrate influx  at 2 h in isolated root cell types (Gifford et al., 2008) and at the phenotypic level to look at longer-term effects of nitrogen on root development (Gifford et al., 2013). A number of previous network inference studies have investigated the response of Arabidopsis thaliana to nitrogen and have helped to identify major nitrogen-responsive sensors and hubs. The best understood of these is the nitrate transporter/ sensor NRT1.1 which also plays a key role in signal transduction during nitrogen responses (Ho et al., 2009). During the early nitrate response, NRT1.1 is responsible for the activation of numerous genes involved in nitrate assimilation. There is also strong evidence that NRT1.1 represses lateral root branching in response to depleted nitrate levels by diverting accumulating auxin from lateral root primordia (Bouguyon et al., 2015). Two bZIP family transcription factors that are downstream of NRT1.1, TGA1 and TGA4, have been demonstrated to regulate expression of the nitrate transporters NRT2.1 and NRT2.2 (Alvarez et al., 2014). It has also been shown that the key circadian clock gene CCA1 is regulated by organic nitrogen signaling and in turn regulates genes involved in nitrogen assimilation itself (Gutiérrez et al., 2008).
However, each of these studies captures a specific snapshot of nitrogen responsiveness and is insufficient for reconstruction of the temporal link between molecular and developmental responses. In order to address this important question, we tracked gene expression responses in roots over a 48-h time series in Arabidopsis. This allowed us to compare temporal changes in the pericycle to those in the cortex, which has potential for postdifferentiation development (e.g., nodule development in leguminous species), but does not divide during lateral root development. We used two environmental perturbations (replete nitrate and rhizobia) that each affect lateral root development and asked to what extent similar changes in gene expression are involved in similar developmental responses. This unique high-resolution time series analysis enabled an investigation of the links between molecular and developmental dynamic behaviors. Our data suggest that the responsiveness of different cell types to the environment is a fundamental functional aspect of a multicellular system such as the root.

Cell-Type Profiling of Arabidopsis Roots Reveals Dynamic Gene Expression over Time
To characterize the transcriptional changes that underlie changes in lateral root development in response to environmental changes, we subjected Arabidopsis root cell-type-specific GFP marker lines to tissue-specific transcriptional analysis after either treatment with nitrogen, inoculation with rhizobia, or no treatment (as a baseline). To investigate transcriptional changes underlying or preceding developmental regulation, we analyzed the pericycle cells using a GFP marker specific to xylem pole pericycle cells from where lateral root primordia originate via reactivation of cell division in differentiated cells (Gifford et al., 2008). We also analyzed a GFP marker line specific to cortical cells (Brady et al., 2007), since this is a more outer cell type in the root, providing the opportunity to examine spatial variation in responses across root cell types. A time series of sampling was performed over a 48-h growth period starting at dawn (= time point 0) on day 9 postgermination, using protoplast generation and FACS (Gifford et al., 2008) (Figures 1A to 1E), and the whole experiment was performed in at least biological triplicate. Each biological replicate sample (for each time point in each time series and cell type) consisted of total mRNA isolated from ;10,000 FACS-isolated cells from pooled harvested roots from ;200 GFP-expressing seedlings grown on the same plate. Each (destructive) sample was considered as an independent replicate within subsequent gene expression and network inference analysis (see Methods).
We first investigated changes in transcriptional levels over the untreated 48-h time series starting with the first harvest of cells at dawn to establish how gene expression changes over time within the cortex or pericycle. We found that 6432 transcripts were differentially expressed (DE) over time in either untreated cortex (CU; 2041 DE transcripts) or pericycle (PU; 4682 DE transcripts) cells ( Figure 1F; Supplemental Data Set 1 and Supplemental Figures 1A to 1C). Twice as many transcripts change in the pericycle as the cortex. To investigate the range of differential expression in time profiles, we clustered DE transcript expression within each time series using SplineCluster (see Methods). This generated 95 clusters for the CU time series and 106 clusters for the PU time series (Supplemental Figures  1E and 1F). To help investigate timing changes, we grouped clusters chronologically into 4 classes of gene response speeds that we term for discussion purposes: "immediate" (by 1 h), "fast" (1-2 h), "downstream" (2-6 h), and "late" (6+ h), determined by the time of the first significant time point ( Figure 1G). We further divided clusters into those whose expression was generally increased (upregulated) or decreased (downregulated) compared with the 0-h time point (see Methods). Clusters are clearly demarcated over time. For instance, the immediate responders show a dominant activation/inhibition within 1 h, with gradually slower changes being evident through the other three classes (Supplemental Figures 1E and 1F). In the cortex, 29% of gene expression changes were immediate, 21% fast, 42% downstream, and 7% late ( Figure 1H; Supplemental Figure 1E). By contrast, in the pericycle, 97% of varying genes showed a substantial expression change within 2 h; 48% are immediate (E) Sample preparation: At each time point, roots are harvested and subjected to protoplast generation treatment, cells are sorted using FACS, and then RNA is extracted, labeled, and hybridized to microarrays (diagram adapted from Grønlund et al., 2012). responders, with 49% being fast, 2% being downstream, and 1% being late ( Figure 1I; Supplemental Figure 1F). For each time series' clusters, there were a wide range of differential expression patterns. These include transcripts whose expression increases over time (e.g., fast upregulated cluster CU62; Supplemental Data Set 2), decreases over time (e.g., pericycle immediate downregulated cluster PU37; Supplemental Data Set 2), as well as transcripts that cycle over time (e.g., putative circadian clusters CU37 and PU64; Supplemental Data Set 2). Clustering also delineates transcripts with the same expression timing but different expression amplitude (e.g., cortex immediate upregulated clusters CU47 and CU50 or pericycle immediate upregulated clusters PU87 and PU88; Supplemental Data Set 2). Overall then, the pericycle appears to have more frequent and earlier (in the day) perturbations in gene expression, potentially reflecting the fact that the pericycle is comprised of cells at a wide range of mitotic activities (Malamy and Benfey, 1997).

A Differentially Expressed Core of Genes Coordinates Changes in Cortex and Pericycle
Despite the large numbers of DE transcripts in the cortex and pericycle, there is a very high degree of cell specificity with only 291 "core" transcripts (transcripts that are DE in both cell types) differentially expressed (as a function of time) in both of the cell types ( Figure 1F). The expression responses of the core transcripts reflect well the range of patterns of all transcripts that are differentially expressed in cells; despite accounting for only 14% (cortex) and 6% (pericycle) of the DE genes, these core 291 transcripts were present in 63% of the clusters. A similar percentage of clusters were reflected from each cell type, invariant of cluster size and including clusters in all four response timing groups. This indicates that this core of DE transcripts is part of a wide range of responses across time in both cell types and could thus act to coordinate the function of these cell types in the context of the whole root, including housekeeping processes such as primary metabolism.
Among these 291 core dynamic transcripts, the Gene Ontology (GO) term "circadian rhythm" was overrepresented (P < 9.15E-4), including five previously characterized core circadian clock genes (Hsu and Harmer, 2014). LHY (LATE ELONGATED HYPOCOTYL), CCA1 (CIRCADIAN CLOCK-ASSOCIATED1), PRR5 (PSEUDO-RESPONSE REGULATOR5), PRR9, and EARLY FLOWERING4 are found to be expressed in accordance with their characterized dynamics and are invariant between cell types, supporting the timeseries nature of our data (Supplemental Figure 2 and Supplemental Data Set 2). The core group also includes transcripts important in a range of metabolic and cellular processes that could be considered to be important for root function, independent of cell type (Supplemental Data Set 2). For example, the nucleosome assembly protein NRP1 (NAP1-RELATED PROTEIN1), together with NRP2, has already been shown to be critical for cell cycle regulation in roots (Zhu et al., 2006). TET1/TRN2 (TETRASPANIN1/TORNADO2), a gene that is expressed throughout the root and plays a role in determining cell fate in the root apical meristem  is regulated by PRR5 and ERF115 (ETHYLENE-RESPONSIVE TRANSCRIPTION FACTOR115) . Consistent with this, we found TET1 and PRR5 to be coregulated (fast cluster CU62 for PRR5 and downstream cluster CU66 for TET1; fast cluster PU79 for PRR5 and downstream cluster PU76 for TET1; Supplemental Data Set 2).
To determine if the 291 core transcripts also share the same temporal dynamics between the pericycle and cortex, we compared their expression using Gaussian process modeling (see Methods). We found that 40% have identical expression patterns (invariant; including clock genes CCA1 and LHY; Figure 1J), 40% have an expression pattern that is similar (e.g., MICROTUBULE-ASSOCIATED PROTEIN 65-1; fast downregulated cluster CU14 and immediate downregulated cluster PU46; Supplemental Data Set 2), and only 20% have expression whose pattern varies significantly between the cortex and pericycle (e.g., RAB GTPase HOMOLOG G1; immediate downregulated cluster CU73 and immediate activated cluster PU87; Supplemental Data Set 2). Together, this suggests that if transcripts are differentially expressed in different cell types, they are often under the same directional control. This makes sense for critical functions such as the circadian clock, whose invariant expression between cell types enables other functions to be coordinated across the root as a whole.

Reprogramming of Cell-Type Responses after Nitrogen Influx
The plant root response to nitrogen has previously been found to be highly cell type specific; however, this conclusion was made from work at a single time point after nitrogen treatment (Gifford et al., 2008). Other work has analyzed an early time series of nitrogen responses to create a functional timeline, but at the whole-root level where signals from several cell types are mixed . To investigate how nitrogen responses are partitioned across a longer-term environmental exposure, we followed the same experimental design as above, but instead, immediately after the 0-h time point, seedlings were treated with replete nitrogen by transferring them to plates containing 5 mM (J) The core regulated transcripts include key clock components, e.g., LHY and CCA1, which are differentially expressed over time in an invariant fashion between cell types (late downregulated CU37, downstream downregulated PU64); error bars, SE with an average sample size of 2.95 replicates. NH 4 NO 3 . In this case, we used a Gaussian process method (GP2S) to determine which of the dynamic transcripts change upon N treatment relative to the untreated time series. We again found a very high level of cell type preference in the N response, with 2404 transcripts N-regulated only in the cortex (CN), 2382 transcripts N-regulated only in the pericycle (PN), and 450 transcripts (only 16% of total) N-regulated in both cell types (Figures 2A to 2C). All differentially expressed transcripts were clustered using SplineCluster, generating 42 clusters for CN and 29 for PN (Supplemental Figure 3).
Supporting the effect of our N treatment, we detected many N-response markers (Canales et al., 2014) in our data (13 only in cortex, 7 only in pericycle and 11 in both [core]; Supplemental Data Set 2). These N-response markers are typically found in whole roots; thus, their cell-type responses here suggest that they are strongly regulated in those cell types. Many of these transcripts are in clusters PU6/PU11 and CU14 (Supplemental Data Set 2), suggesting that these represent N-responsive modules. By analyzing the expression regulation of key components of nitrogen uptake (nitrogen and ammonia transporters) and metabolism (primary nitrogen metabolism), it is possible to analyze the temporal aspects of nitrogen uptake across the tissues ( Figure 2B). Both the high-affinity nitrate transporter NRT2.1 (NITRATE TRANSPORTER2.1) and the dual-affinity nitrate transporter NRT1.1 as well as nitrate reductase activity (NIR1) are immediately upregulated specifically in the cortex ( Figure 2B). The rapid upregulation of NRT2.1 in the cortex (cluster CN3, immediate; Supplemental Data Set 2) corroborates previous expression data (Feng et al., 2011) and is consistent with the role of NRT2.1 in signaling rescue from N starvation, preparing the plant to start absorbing N again (Yong et al., 2010). At the same time ("immediate"), nitrite reductase (NIA1) and other nitrogen transporters (NRT2.3 and NRT2.5) are downregulated in the pericycle, suggesting a tissue-specific effect. Ammonium transport (via AMT1;1) appears to be downregulated in both tissue types. Nitrogen signaling via the key N-responsive hub bZIP1 (Obertello et al., 2010) is immediately upregulated in both cortex and pericycle. Other N-responsive regulators that we found include the GARP transcription factor HRS1, which integrates N and P signaling in the root and is N-induced late in the pericycle (cluster PN8; Supplemental Data Set 2) (Medici et al., 2015). SAG21 (SENESCENCE-ASSOCIATED GENE21), a regulator of both primary and lateral root development as well as biotic responses (Salleh et al., 2012), is downregulated in the pericycle late after N treatment (cluster PN11; Supplemental Data Set 2).
Another line of support that the N responses we see are physiologically relevant and related to molecular changes that precede regulation of lateral root development comes from the presence of a NIN1-like putative transcription factor, NLP8, among the core N responses (clusters CN6 and PN11; Supplemental Data Set 2). NIN is a critical regulator that coordinates cell-type responses during nodulation in the legume Medicago truncatula (Vernié et al., 2015). NLP8 has recently been found to regulate N-responsive seed germination (Yan et al., 2016). We analyzed T-DNA lines and found that a nlp8 homozygous loss of expression knockout mutant (nlp8-2 as named in Yan et al., 2016) has a shorter primary root and fewer lateral roots on replete levels of nitrate ( Figure 2F; Supplemental Data Set 3) but not on depleted levels of nitrate. This is consistent with findings by Yan et al. (2016) who did not analyze root development but did observe N-dependent effects on germination.

Cell-Type Preferential Nitrogen Responses Coordinate
Nitrate Influx with Induction of Lateral Root Development and Remodeling of the Root Our temporal data uncovered the rapid speed with which the Arabidopsis transcriptome responds to a new environment. By examining the cell type preferential N response at the gene and process level (GO term analysis of gene clusters), we can determine the extent to which compartmentalization in different cell types underlies regulation of lateral root development ( Figures  2D and 2H). Within the pericycle, WRKY75, a gene that has previously been shown to affect lateral root numbers (Devaiah and Raghothama, 2007), is strongly and immediately upregulated after nitrogen treatment (cluster PN6; Supplemental Data Set 2). Next, there is N repression of genes regulating "multidimensional cell growth" (fast cluster PN23, P = 1.22E-4; Supplemental Data Sets 2 and 4). A combination of processes are regulated relating to modification of the cell wall during lateral root primordial development, including genes involved in cytoskeleton reorganization, cell elongation, and cell wall biogenesis (e.g., KOBITO, Actin-1, ATFUT5, CESA1, and Profilin-1; all cluster PN23; Supplemental Data Set 2), which can all be directly related to lateral root development (Vilches-Barro and Maizel, 2015) ( Figure 2H). Genes involved in "cellular respiration" (cluster PN25, P = 1.49E-8; Supplemental Data Sets 2 and 4) are induced early by N (at 2 h) and again at 16 h, in a diurnal fashion. This could provide energy during the lateral root growth phase. Late N repression of "RNA methylation" (cluster PN4, P = 2.52E-11; Supplemental Data Sets 2 and 4) could be linked to development of new lateral root primordial cells since it follows the first cell division in the pericycle. Also late N-induced are the auxin biosynthetic genes YUC3 (YUCCA3; PN15; Supplemental Data Set 2) and YUC7 (PN18; Supplemental Data Set 2) and the N-responsive transcription factor WRKY15 (PN6; Supplemental Data Set 2). Using mutant analysis, we found that a wrky15 homozygous loss of expression mutant, relative to the wild type, has a longer primary root and more lateral roots on deplete N but a similar root system on replete N ( Figure 2G; Supplemental Data Set 3). This confirms that our analysis is able to identify novel functionally important regulators of N-responsive lateral root development.
Our time-series data showed that the cortex also behaves dynamically in response to nitrogen. There is immediate upregulation of N-responsive genes including NRT2.1 (CN3; Supplemental Data Set 2), NRT1.1 (CN27; Supplemental Data Set 2), and NRT3.1 (CN27; Supplemental Data Set 2), and cluster CN28 is overrepresented for "nitrate transport" processes (P = 3.36E-5; Supplemental Data Set 4) ( Figure 2D). NRT1.1 directs lateral root growth to patches of high N in soil (Remans et al., 2006) and acts as a regulator for NRT2.1 (Muños et al., 2004), which itself forms a N-transporting complex with NRT3.1 (Yong et al., 2010;Kotur et al., 2012). We also found immediate N regulation of transcriptional regulators including the MYB-like transcription factor (TF) At1g25550 (cluster CN28; Figure 2D; Supplemental Data Set 2). At the same time, the cortex appears to be rapidly remodeled, possibly to accommodate lateral root emergence from underlying cell layers. Cytoskeleton reorganization and vesicle trafficking are  (B) Euler diagram showing the numbers of genes regulated by nitrogen in the cortex (purple font), pericycle (red font), or in both (core N-responding genes, blue font). Nitrogen assimilation genes that were found to be nitrogen regulated in the cortex, pericycle, or both (core) are indicated next to the Euler diagram sections; the arrows designate direction of regulation and timing (as in Figure 1G); NRT3.1 is differentially regulated in the cortex and pericycle. (D) Timeline of processes affected by nitrogen, with the log 2 expression of select indicative genes shown in nitrogen-treated and untreated cortical cells. Among cortex clusters that contain an overrepresentation of transcripts involved in particular processes are immediate N responses, including upregulation of nitrate responses (cluster CN28), downregulation then recovery of cytoskeleton organization components (cluster CN37), and biphasic upregulation of cell wall modification genes (cluster CN4). Solid black line, untreated (CU); dashed purple line, N-treated (CN); see Supplemental Data Set 4 for GO overrepresentation analysis corrected P values and Supplemental Data Set 2 for cluster membership. (E) N treatment induces lateral root development, as visible if seedlings are maintained in the treatment up to 4 d. Bar = 1 cm. (F) Perturbation of core NIN1-like putative transcription factor NLP8 in nlp8-2 (cluster CN6, PN11) leads to seedlings with shorter roots specifically on replete nitrogen. Bar = 1 cm. immediately N-repressed ( Figure 2D); this includes "vesicle coat" processes (P = 9.67E-3; Supplemental Data Set 4) and two genes involved in microtubule binding in cluster CN37 (At3g45850 and At1g65470; Supplemental Data Set 2), suggesting a change in cell growth processes. PIN-FORMED2, an auxin efflux carrier, was found to be fast N-induced (in cluster CN36; Supplemental Data Set 2). There is also an early biphasic (1 and 14 h) induction of "cell wall modification," overrepresented in cluster CN4 (P = 1.7E-2; Supplemental Data Sets 2 and 4), here involving N regulation of expansins (e.g., ATEXPB1; cluster CN39; Supplemental Data Set 2), together with an increase in cell adhesion processes ( Figure 2D).

Network Inference Reveals How Nitrogen Responses in the Cortex and Pericycle Are Connected via Regulation of Transcription Factor Activity across Cell Types
To identify candidate upstream regulators of the nitrogen and lateral root development response and understand how regulated processes are connected, we generated causal transcription networks using the Bayesian network inference tool GRENITS (Morrissey et al., 2012;Wang et al., 2014) from DE genes in each experiment. We assessed the accuracy of connection predictions by comparing predicted edges in all networks with the interactions found in published data (O'Malley et al., 2016) (see Methods; Supplemental Data Set 5) and by asking if the promoters of target genes contained a binding site motif for the predicted regulatory transcription factor using a test set of edges between nodes that were common between all interaction networks in a cell type (see Methods; Supplemental Data Set 6D). On average, 43% of DE genes were connected in networks and all temporal patterns of expression were represented (see Supplemental Files 1 to 14 and Supplemental Table 1). To gain a higher-order "backbone" view of regulatory connections, we created a view of the network by forming nodes from transcription factors together with their targets (each being a "TF module"), based on the fact that one TF was predicted to regulate another TF and thus the TF module ( Figure 3; Supplemental Files 7 to 14). This enabled us to query each module to ask if any functions or processes were significantly overrepresented. If a module, or connected modules, shared similar overrepresented GO terms, we shaded them to illustrate potential common function ( Figures 3A to  3C). The cortex nitrogen regulatory backbone network contains several such sets of connected modules that are involved in similar biological processes ( Figure 3A). These include abscisic acid (ABA), chitin, circadian, JA, gibberellin, and ethylene responses. Together, this suggests that a wide range of processes are controlled by nitrogen within the cortex. There is considerable evidence linking hormone signaling to N status in the plant. For example, ABAinsensitive and -deficient mutants are less sensitive to the inhibitory effect of high nitrate concentrations on root growth (Signora et al., 2001).
The finding of a decrease in expression levels in clusters of genes enriched for JA GO terms agrees with findings that JA genes are upregulated during nutrient starvation and quickly downregulated after resupply (Armengaud et al., 2004). Most of the CN network nitrogen-regulated transcriptional regulators were in the "immediate" response class (30/45), with a similar ratio to the number of differentially expressed genes within CN ( Figure 2A). Module size (number of edges) was found to be largest in the "fast" class (Supplemental Data Set 6). A test of the ranking of the fast modules when ordered by size indicates that larger modules significantly correlate with the fast category (P < 3.8E-2 using a Monte Carlo simulation for the null distribution). The largest module in the CN network with 34 targets is PIF7 (fast cluster CN10; Supplemental Data Set 2), a member of a basic helix-loop-helix-type TF family that has been linked to transcriptional activation leading to a rapid growth response (Li et al., 2012).
In the pericycle, the regulatory backbone network is more sparsely connected, with some significant modules that could be related to formation of new lateral roots, including N responses, ribosome biogenesis, organ structure development, and establishment of cell localization ( Figure 3C). As with the CN network, we again see that the "fast" modules are the largest (P < 3.3E-3 using a Monte Carlo simulation for the null distribution). The largest regulatory module in the PN network, with 77 regulatory targets, is At1g31760 (fast cluster PN5; Supplemental Data Set 2), an unknown member of the SWIB/MDM2 domain superfamily of chromatin remodeling proteins that facilitate transcription activation. The second largest module is SUVH4 (fast cluster PN5; Supplemental Data Set 2), a histone methyltransferase involved in the maintenance of DNA methylation, suggesting that posttranscriptional modification could play a crucial role in the regulation of the response to nitrogen in the pericycle. SUVH4 inhibits RAB HOMOLOG1 (late cluster PN11; Supplemental Data Set 2), a gene involved in the general response to auxin in roots, with mutants exhibiting defective gravitropism and auxin physiology as well as fewer lateral roots (Fortunati et al., 2008).
Although some genes with the same GO terms or molecular processes are N-responsive in both the cortex and pericycle, the lack of shared regulatory links between the cortex and pericycle networks suggests that individual N-regulated transcripts and  processes in each cell type are different. However, there are 11 N-regulated transcription factor modules that are present in both cortex and pericycle backbone networks and that thus could enable coregulation of N responses between cell types ( Figure  3C). These are the transcription factors that are regulated in common and could be suggested to act as "core" nitrogen regulators. A predominant cortex then pericycle timing of the N response is evident when comparing the timing of expression regulation of these genes; all but two of these 11 transcription factors are N-regulated in the cortex earlier than in the pericycle.
Downstream of these core transcription factors, there are putative targets with the same molecular function, based on the annotation of the target genes. For example, NF-YA10 regulates F-box proteins in both cortex (At1g62270) and pericycle (At5g39450), and AtbZIP63 (At5g28770) induces a LRR protein in the cortex (At3g15410) and represses a LRR kinase in the pericycle (At1g49100) ( Figure 3D; Supplemental Files 2, 5, 13, and 14). Individually observing these transcription factors regulating LRR and F-box family genes is not significant on the inferred CN/PN network (P = 0.056 for LRR and P = 0.085 for F-box; (A) and (B) Inferred causal transcription factor (backbone) network on CN (A) and PN (B) expression data. Node size is scaled by the outdegree of the module (number of targets) and colored according to response category: immediate (red), fast (orange), downstream (yellow), and late (blue) as in Figures 1 and 2 (see Supplemental Files 1 to 9). The edges denote that one TF regulates the other TF (and thus TF module), with regulation type (induction; arrowheaded lines or repression; bar-headed lines) shown by the arrowhead shape. Processes that are overrepresented among genes within module(s) are highlighted for each network. (C) Transcription factors that are N-regulated in both cell types (core) have been located between the cortex and pericycle networks; within the solid box, each pair of nodes is the same transcription factor with the regulation timing and node size indicating interactions in the cortex (left) and pericycle (right). Directionality of the N response is evident in the earlier response of these in the cortex compared with the pericycle (see colors designating timing, as in Figure 1G). (D) Targets of bZIP63, transcription factor in dotted line box section in (C). Expression profiles within time series of a transcription factor (At5g28770, bZIP63) that is core N-induced in both the cortex (fast) and pericycle (late), and predicted to regulate genes of similar function in each cell type, LRR protein At3g15410 in cortex (immediate) and LRR protein kinase At1g49100 in pericycle (late). Solid line, untreated (CU/PU); dashed line, N-treated (CN/PN); error bars indicate SE with average sample size of 2.95 replicates. permutation test). However, observing both in the CN/PN networks is significant at P = 0.0043 (permutation test). Therefore, our time-series analysis has revealed that a combination of variation in the timing, direction and location of core common responses can potentially regulate and enable cell type preferential outputs. To enable this common signaling, it is possible that transcription factors move between cell types (as SHR/ SCR protein moves from the stele to the endodermis; Nakajima et al., 2001), but it is also possible that TF regulation is independent in both cell types.

Arabidopsis Biotic Responses Can Be Partitioned into Defense Responses in the Cortex, and Subsequent Developmental Responses in the Pericycle
Having discovered that the cortex and pericycle coordinate changes in response to the abiotic change of nitrogen influx, we asked whether coordination was also evident for biotic responses that alter root architecture. We analyzed Arabidopsis responses to the M. truncatula (legume) N-fixing bacterium symbiont Sinorhizobium meliloti, a biotic stress that we discovered can alter root architecture. S. meliloti is not a symbiont of Arabidopsis, but S. melilotitreated seedlings have more lateral roots (P = 1.09E-2), although these are shorter (P = 1.77E-2) (Figures 4E to 4G; see data in Supplemental Data Set 3), indicating a possible stress/immunity response. We performed a time-series analysis of S. meliloti treatment and identified 2546 transcripts that respond specifically in cortex, 3350 that respond specifically in pericycle, and 202 that respond in both cell types ( Figure 4A). CR transcripts were clustered into 32 clusters and PR transcripts into 25 clusters (Supplemental Figure 4). As with the nitrogen responses, we used network inference to construct networks to identify transcription factors that are associated with the regulation of rhizobia responses (Supplemental Files 3 and 6).
We found putative defense responses in the cortex to be immediately activated then maintained, with three immediate/ fast regulated clusters enriched for genes annotated as chitin responses ( Figure 4B, clusters CN19, 28, and 30; Supplemental Data Set 2; see P values in Supplemental Data Set 4). Rhizobia contain microbe-associated molecular patterns (MAMPs) such as peptidoglycan that activate innate immunity via a conserved LysM signaling pathway that also responds to the MAMP chitin (Gust, 2015). Among the "chitin-responsive" genes are several cortex-specific rhizobial-responsive TFs including ERF11 (cluster CR19, Supplemental Data Set 2); WRKY48 (cluster CR30, Supplemental Data Set 2); a known repressor of plant basal defense genes (Xing et al., 2008); and MYB44 (cluster CR30, Supplemental Data Set 2), a positive regulator of SA-associated defense responses and a negative regulator of JA-associated defense response . Within the inferred cortex rhizobia-responsive regulatory network (Supplemental File 9), highly connected predicted defense modules included WRKY3 (cluster CR7; Supplemental Data Set 2) and WRKY28 (cluster CR28; Supplemental Data Set 2), with CESA1 (cluster CR26; Supplemental Data Set 2), a cellulose synthase critical for cell wall formation (Burn et al., 2002) predicted to control developmental responses (see Supplemental File 9). Modification of the responses of rhizobia-responsive genes (e.g., via mutation) could alter basal defense responses to bacteria such as rhizobia, as already seen for responses to pathogenic microbes. For example, wrky48 mutants exhibit stronger defense responses (Xing et al., 2008) and myb44 mutants exhibit higher pathogen resistance and are also more sensitive to JA-mediated root growth inhibition .
Later in the cortex there is upregulation of stress signaling, including "response to jasmonic acid stimulus" (cluster CR10, P = 4.27E-4; Supplemental Data Set 2), and induction of genes known to be involved in remobilizing nitrogen (e.g., NAXT1, cluster CR26; Figure 4B; Supplemental Data Set 2), which could indicate repartitioning of resources. Similar to in the nitrogen response, cell wall development is regulated, e.g., the methyltransferase GXMT (cluster CR29; Supplemental Data Set 2) that plays a role in the production of cell wall xylan (Lee et al., 2012), CESA3 (cluster CR27; Supplemental Data Set 2), and CSI1 (cluster CR28; Supplemental Data Set 2), which are both involved in cellulose production and for which mutation affects root development (Daras et al., 2009;Gu et al., 2010). "Plant-type cell wall modification" related transcripts are also regulated in the pericycle, e.g., the pectin methylesterase inhibitor At2g10970 (cluster PR6, P = 8.3E-4; Figure 4B; Supplemental Data Sets 2 and 4), which could be related to the induction of lateral root development or cell wall fortification as part of immunity.
In summary, we found that S. meliloti induces defense responses at early time points (immediate/fast) and stress signaling (late). In addition, other processes including cell wall modification and nitrate remobilization are rhizobia responsive. We hypothesize that combined, these lead to effects on root architecture as a form of stress adaptation ( Figures 4C to 4G).

Cell-Type Identity and Cell-Type Responses to the Environment
Around half of the transcripts that are differentially expressed in response to an abiotic (nitrogen) or biotic (S. meliloti) treatment are also differentially expressed in untreated cells over time. This suggests that many genes that respond to the environment are dynamically expressed within cell types, rather than only induced or repressed upon treatment. Cell types are known to have varied responses to the environment, at least in part due to the different patterns of gene expression within those cell types (Gifford et al., 2008). At the same time, just as it is well known that plant root cell position confers cell fate (Kidner et al., 2000), cell "state," meaning the functions and responses of a cell, are likely to be influenced or determined by signals from outside the environment of a cell type. To understand the extent to which cell-type identity determines the scale of cell-type responses to the environment (cell-type state), we evaluated the expression of 194 cortical and 104 pericycle transcripts previously defined as cell-type-specific enriched genes (CTSEs) (Birnbaum and Kussell, 2011;Bargmann et al., 2013) in our data. We identified CTSEs appropriate to each cell type and found that they were enriched in the cortex and pericycle, validating the cell-type identity of our sorted cells ( Figure 5A). However, when analyzing CTSE gene expression, we also found that CTSEs were dynamically expressed, either changing within untreated time series or in response to nitrogen and/or rhizobia ( Figures 5B and 5C; Supplemental Data Set 2). This suggested that the CTSEs were indeed cell-type-specific enriched but that they were also responsive to the environment.
To determine what controls the expression of the CTSEs that respond to the environment, we queried the PN causative network for N-responsive pericycle CTSEs ( Figure 5E) and the CR network for S. meliloti-responsive cortex CTSEs ( Figure 5D). Two SCL transcription factors are predicted to repress the pericycle response of the BR-6-ox1 N-responsive gene, potentially repressing brassinosteroid synthesis under replete N ( Figure 5E). Four transcriptional regulators were found to coregulate a group of three cortex rhizobia-responsive CTSEs and the expression of three out of four of these regulators was found to be similarly enriched in the cortex ( Figure 5F). Overall, we found that responsiveness to the environment is enabled by the transcriptome of a cell type and that this response also contributes to the state of that cell type. Using network analysis, we were able to identify novel cell-type-enriched responsive genes that could control responses of particular cell types.

Individual Genes Are Expressed with a High Degree of Temporal and Spatial Specificity, yet the Regulation of Processes Has a High Degree of Overlap
Although most differentially expressed transcripts were DE in one cell type and treatment only (Figure 6), our gene and network analysis suggests that similar processes are regulated across multiple conditions. To confirm this, we determined the degree of functional overlap between lists of differentially expressed transcripts. We assigned each transcript to its deepest nonoverlapping path molecular function GO term (see Methods). There is a significantly higher proportion of transcripts that are unique to a cell type and treatment than the proportion of unique GO terms (P < 10E-116), while the proportion of GO terms shared among all cell types and treatments is significantly higher than the proportion of shared transcripts (P < 10E-84) (z-tests) ( Figures 6A and 6B); transcript level overlap is low (62% DE in a single cell type and treatment; 1% DE in all cell types and treatments) compared with process level overlap (31% DE in a single cell type and treatment; 23% DE in all cell types and treatments). We term this mechanism "response nonredundancy," whereby individual transcripts within gene functional groups are expressed with a very high degree of temporal and spatial specificity, yet biological processes are commonly regulated.
Although most (84%) N responses are cell-type preferential, those that are shared between cell types (core N responses) could shed light on the coordination of the N response. We analyzed the timing of the responses of the 450 core N-responsive genes. We found that 24% (108) of these transcripts display invariant expression in cortex and pericycle. Assessing the 342 genes with timing directionality, a significantly larger number of the core N-responsive transcripts (295) are first regulated in the cortex, then in the pericycle ( Figure 6C; P < 10E-20, binomial test). This suggests that N responses occur first in outer root cell types, including the cortex, before then occurring in inner cell types including the pericycle ( Figure 6C). Previous work suggests that nitrate is rapidly transported throughout the root to regulate core N responses, with organic assimilated N typically playing a role in regulating cell-type responses including lateral root development (Gifford et al., 2008). The signal to regulate tissue directional responses could thus be assimilated N or a molecule downstream. This would likely involve signaling across the endodermis, also an N-responsive cell type (Gifford et al., 2008).
The N response of 135 of the 295 cortex-then-pericycle transcripts are in the same regulation direction (induced or repressed; see Methods). This timing of gene regulation suggests that an N-induced signal regulates gene expression in the cortex, before passing through and inducing a N-mediated response in the pericycle. This category is overrepresented for genes involved in responses to abiotic stimulus (P = 0.001), which includes many N-and ABA-responsive genes ( Figure 6C; Supplemental Data Set 4). Also present is the bZIP family transcription factor TGA1, which has been demonstrated to act downstream of NRT1.1 and has been implicated in the regulation of NRT2.1 and NRT2.2 (Alvarez et al., 2014). The remaining 160/295 genes whose N-response occurs later in the pericycle are regulated in opposite directions in cortex and pericycle. This includes two genes linked to auxin responses: EXTENSIN (clusters CN38 and PN23; Supplemental Data Set 2), an auxin-responsive member of environment-driven cell wall modifying enzymes (Eklöf and Brumer, 2010); and ZINC INDUCED FACILITATOR-LIKE1 (clusters CN6 and PN11; Supplemental Data Set 2) that has been shown to modulate auxin transport (Remy et al., 2013). These two genes could be part of auxin-responsive lateral root regulation, responding to auxin flow from the vasculature.
Systemic, non-cell-autonomous signaling or a mix of both could potentially control 108 of the 450 common DE transcripts, since they are regulated at the same time in both cell types ( Figure 6C). Seventy-three of these transcripts have identical expression, including bZIP1 that is immediately N-induced in both cortex and pericycle (clusters CN31 and PN6; Supplemental Data Set 2), and the key ammonia transporter AMT1;1 (AMMONIUM TRANSPORTER1;1) that is N-downregulated fast in both cell types (fast clusters CN10 and PN24; Supplemental Data Set 2). Our finding that bZIP1 plays a role in rapid and early N signaling is supported by studies of bZIP1 perturbation in studies of isolated root cells  (Para et al., 2014). The remaining 47/450 common DE transcripts are regulated first in the pericycle and second in the cortex, which could indicate cell type independence rather than core responses ( Figure 6C).
We observed widespread response nonredundancy regulation of genes involved in the same gene ontology molecular functions ( Figures 6A and 6B), and it is also evident when assessing the regulation of gene family members ( Figure 6D). We analyzed the regulation of 55 gene families and found that 47/55 gene families had members differentially expressed in each time series. The remaining eight had family members differentially expressed in at least three time series, but these had fewer members (an average of eight DE transcripts rather than an average of 35 DE transcripts; Supplemental Data Set 8). As an example of the typical widespread differential expression within gene families, we found that 29/51 transcripts from the glutathione S-transferase (GST) family are regulated in one or more time series and show diverse patterns of gene regulation ( Figure  6D; Supplemental Data Set 7). Although variation in expression patterns between members of a gene family could reflect the inherent range of functions, it is also likely that regulation of the same molecular function in a spatially and temporally diverse fashion will alter the outcome of that molecular activity. This represents the same phenomenon as discussed earlier for common nitrogen-responsive network modules: These were found to regulate the same gene families but control different individual members of these gene families, specific to response and location ( Figure 3C). Together, this shows how variations in temporal and spatial regulation in gene expression can enable specificity of responses and highlights the importance of large gene families for enabling complex responses to the environment.

DISCUSSION
Our novel time-series data have enabled us to show that environmental responses exhibit a very high degree of temporal and tissue specificity. By tracking the expression of transcripts over time, we were able to reconstruct a timeline of regulation of molecular processes that link short-term and long-term consequential developmental changes in the root system. We found that gene expression in the cortex and pericycle of the root is highly dynamic over day/night cycles, with thousands of genes changing in a wide range of different temporal patterns. Many of these genes also respond to the environment, potentially enabling both "enriched," no differential expression within a time series; "dynamic," differential expression within the untreated (U) time series; "responsive," differential expression within the nitrogen-treated (N) or rhizobia-treated (R) time series; purple shades, differential expression within cortex; red shades, differential expression within the pericycle. (D) and (E) Network inference-predicted upstream regulators of enriched and responsive CTSEs in the cortex (D) and pericycle (E); diamonds, genes annotated as transcription factors, circles, genes with non-transcription factor functions. Color denotes time of first significant differential expression (as in Figure 1G).  temporal-and spatial-specific responses. We hypothesize an interplay between cell-type identity and how these cell types respond to their environment, whereby responses are programmed by the cell type, but in turn reinforce the unique functions of those cell types (the cell type is programmed by its responses). The validity of this hypothesis could be investigated by expressing genes that are typically responsive in one cell type (e.g., pericycle), under the control of a promoter typically expressed in a different cell type (e.g., the promoter of a cortex-specific gene) and/or typically responsive in a different cell type (e.g., the promoter of a cortex-responsive gene). Such "promoter swapping" to compare responses would help to test our hypotheses regarding the connections between cell-type identity and cell-type responses.
We identified transcription factor modules that control modules of processes and also identified a novel root architecture role for NLP8 (Yan et al., 2016) and a new root architecture-regulating transcription factor, WRKY15, the functions of which are moderated by N abundance. The molecular regulation that we uncover underlies lateral root developmental responses under different environments, one component of root plasticity. The time series also enabled us to identify novel directionality in root responses between tissue layers to nitrogen.
Analysis of differentially expressed genes shows that there is a high level of specificity in the regulation of genes in response to treatment and within each cell type. Despite the fact that thousands of genes respond to treatments, each cell type before and after treatment still clusters together, suggesting that cell type identity is maintained during these responses. CTSE genes (Birnbaum and Kussell, 2011;Bargmann et al., 2013) are expressed in particular cell types but their expression is not static: They are highly responsive to the environment. This implies that cell type responses are controlled by the expressed transcriptome, but also that transcript responses contribute to modulating cell-type state. Within this analysis we identified new enriched genes using differential expression and network analysis. In our networks, we also found closely connected modules regulating similar processes, akin to "molecular network machines," a term for identifying common control of molecular processes, evident within network analysis (Gunsalus et al., 2005). Our analysis revealed a novel level of coordination between cell types, since several of the biological processes that are regulated by nitrogen treatment or rhizobia inoculation in cortical and pericycle cells were the same, but different genes are responsible for the responses in each cell type. For example, common transcriptional modules respond in both pericycle and cortex, and regulate similar gene families, but different members of these families act in each cell type and at different times during environmental responses. This suggests much less functional redundancy in gene families than is typically considered and that this is obscured at the whole-root level. This is an important finding since it shows that the scale and scope of gene activity is not simply the molecular function of a gene, but a consequence of its spatial location, timing of expression, and responsiveness to environmental conditions. It also helps to provide an explanation for the typically higher hit rate in identifying phenotypes in regulators identified in systems biology data set analysis using a reverse genetic approach (Ransbotyn et al., 2015); a higher level specificity in the measurement of expression levels, regulation, and location of genes means that investigators can determine more precise and accurate predictions of when perturbation of function might lead to a measurable phenotype.
Overall, these findings identify novel mechanisms of environmental responses and highlight the need to examine the expression of individual genes at the temporal and spatial levels to better understand how developmental plasticity is controlled and enabled in plants.
For GFP-expressing seeds, ;400 per plate were surface-sterilized and sown on 0.7-cm strips of autoclave-sterilized brown growth pouch paper (CYG germination pouch) on 0.8% agar plates containing a basal 13 Murashige and Skoog growth medium (Sigma-Aldrich M0529), supplemented with sucrose (30 mM), CaCl 2 (1.5 mM), MgSO 4 (0.75 mM), KH 2 PO 4 (0.625 mM), and 0.3 mM NH 4 NO 3 , pH 5.7. Plates were sealed with microporous tape, seeds stratified for 2 d at 4°C in the dark, and then plates placed in opaque (black polythene) covers (Bagman of Cantley) and grown (B) Euler diagram showing the number of molecular function GO terms from nonoverlapping paths represented within differentially expressed transcripts in (A) and illustrating the greater overlap at the function level compared with at the transcript level. The percentage of transcripts DE in one/two/three/all cell types/treatments is 62:33:4:1, respectively, and the percentage of GO terms DE in one/two/three/all cell types/treatments is 31:29:17:23, respectively. (C) Transcripts that are N-responsive in both the cortex and pericycle (core N responses) can be distinguished according to the timing of their regulation. The predominant pattern is cortex then pericycle regulation (295 transcripts), then regulated in both cell types at the same time (108 transcripts), followed by pericycle then cortex regulation (47 transcripts). Overrepresented processes or genes of note are indicated below indicative expression patterns for each response category; see Supplemental Data Set 4 for GO overrepresentation analysis corrected P values. (D) Euler diagram showing the number of GSTs significantly differentially expressed in the cortex or pericycle in response to nitrogen or rhizobial treatment (Supplemental Data Set 7), with the expression of selected transcripts across response categories shown to illustrate the diversity of responses within the same gene family. Log 2 -normalized expression values of all transcripts are shown on the same axes for comparative purposes; error bars indicate SE with average sample size of 2.95 replicates; colored boxes denote response time category as in Figure 1G; see key.
in a Sanyo growth chamber (MLR-351H; Sanyo, E&E Europe) with a 16/8-h photoperiod at 50 mmol m 22 s 21 and constant 22°C. For T-DNA lines and their wild-type siblings/Col-0, approximately eight seeds were surfacesterilized and sown directly, side-by-side on deplete N agar plates (as above) or replete N agar plates (as above but with 5 mM NH 4 NO 3 ) for 12 d, or seedlings were treated with replete N or rhizobia as above at 9 d, then imaged 4 d after treatment.

Plant Treatments
After 9 d in the growth chamber, treatments were performed at dawn (= time 0). To prepare the Sinorhizobium meliloti solution, 50 mL of TY/Ca 2+ medium (5 g/L Bacto-tryptone, 3 g/L yeast extract, and 6 mM CaCl 2 ) (Journet et al., 2001) was inoculated with S. meliloti and incubated for 26 h at 28°C and 220 rpm to an OD 600 of 1 to 2. Cells were then harvested by centrifugation (4000 rpm, 10 min, 4°C) and resuspended in 40 mL water to a final OD 600 of 0.01 (Ding et al., 2008). To carry out the treatments, seedling strips were removed from plates, briefly immersed (10 s) in liquid deplete N medium (as plates but with no agar), placed on a fresh plate (untreated or briefly immersed in liquid replete N medium, as plates but with no agar and with 5 mM NH 4 NO 3 ), and placed on a fresh 5 mM NH 4 NO 3 plate (N treatment), or briefly immersed (10 s) in a solution of S. meliloti and placed on a fresh deplete NH 4 NO 3 plate. Plates were then resealed with microporous tape, replaced in opaque (black polythene) covers, and returned to the Sanyo growth chamber. The same procedure was used to treat Col-0 plants for phenotypic analysis (Figures 2E, 4C, and 4E to 4G), but after growth of ;8 seedlings per agar plate for visualization purposes. At hourly time points 0 (immediately before transfer, at dawn), 1, 2, 4, 6, 8, 10, 12, 14, 16, 20, 24, 36, and 48 h after transfer, whole roots were harvested for protoplast generation and FACS by cutting roots just below the growth pouch paper strip. The whole experiment was performed in at least biological triplicate with seedlings for each biological replicate grown independently (in a nonoverlapping time period). Each replicate sample (for each time point in each time series) consisted of the pooled harvested roots from ;200 GFP-expressing seedlings (destructive sampling) grown on the same plate. Each replicate sample was considered independently within subsequent analysis.
Protoplasts were sorted using a BD Influx cell sorter (BD Biosciences) fitted with a 100-mm nozzle, using FACS-Flow (BD Biosciences) as sheath fluid. Pressure was maintained at 20 p.s.i. (sheath) and 21 to 21.5 p.s.i. (sample), and drop frequency was set to 39.5 kHZ, which yielded an event rate of <4000; these are optimal settings on a BD influx cell sorter for the type of protoplasts (Gifford et al., 2008). To optimize alignment of relevant lasers and detectors, Calibrite Beads (BD Biosciences) suspended in FACS-Flow were used, and to optimize sort settings, BD Accudrop Fluorescent Beads (BD Biosciences) suspended in FACS-Flow were used. GFP-positive protoplasts were identified using a 488-nm argon laser and plotting the output from the 580/30 band-pass filter (orange) versus the 530/ 40 band-pass filter (green). GFP/YFP-positive protoplasts were in the high 530/low 580 population, with non-GFP protoplasts in the low 530/low 580 population and dead/dying protoplasts and debris in the high 580 population (as in Grønlund et al., 2012). At least 10,000 GFP-positive protoplasts were sorted using methods shown previously to preserve endogenous gene expression levels (Birnbaum et al., 2003;Gifford et al., 2008). Sorted protoplasts were directly sorted into Qiagen RLT lysis buffer containing 1% (v:v) b-mercaptoethanol, mixed, and immediately frozen at -80°C for RNA extraction (methods according to Gifford et al., 2008).

RNA Isolation, qPCR, and Microarray Hybridization
RNA was extracted from sorted cells using the Qiagen RNeasy RNA kit and DNase treatment using the Qiagen DNase kit. The quantity and quality of RNA were checked with a Bioanalyzer 2100 RNA 6000 Pico Total RNA Kit (Agilent Technologies). cDNA was amplified from RNA using 1/2 reactions of the Ovation Pico WTA System V2 kit (NuGEN Technologies) according to the "quick" protocol. Post-amplification purification of cDNA was performed using the QIAquick PCR purification kit (Qiagen). Then, 0.5 mg NuGEN-amplified cDNA was labeled using the NimbleGen One-Color Cy3 labeling kit, and 4 mg of this was hybridized using the GeneChip hybridization kit on NimbleGen Arabidopsis 12 3 135k probe microarrays designed for the full TAIR-10 annotation Arabidopsis genome (design OID 37507; see GPL18735 and GSE91379; Roche Applied Science). OID has 131,524 probes designed against 27,143 genes in the Arabidopsis genome; 3675 genes had multiple probes along their length to determine expression of 31,524 gene transcripts (an average of 2.19 and up to four gene models each for these genes). Arrays were imaged using an MS200 microarray scanner using only the 480-nm laser using the autogain feature of the NimbleScan software, all according to the manufacturer's instructions. Grids were aligned automatically then manually verified. Raw probe (xys) and gene level unnormalized data were obtained using NimbleScan.
For qPCR analysis of gene expression in T-DNA and wild-type lines, root samples (consisting of ;16 roots per sample) were harvested at the end of the growth period on plates; three independent biological replicate samples were harvested. RNA was extracted using the Qiagen RNeasy RNA kit and used for cDNA synthesis (Primer Design). qPCR was performed using SYBR Green dye (Sigma-Aldrich) on a LightCycler 480 system (Roche) (primer sequences in Supplemental Data Set 3D). As a control, we used polyubiquitin UBQ10 using the geNorm REF gene kit (Primer Design) as described (Mestdagh et al., 2009).

Normalization and Quality Assessment of Microarray Data
The xys (NimbleGen) files, which store array coordinates and observed intensities, were imported with the Bioconductor "oligo" package (Carvalho and Irizarry, 2010) using the annotation package pd.120110.athal.mg.expr installed through the pdInfoBuilder package (Falcon and Carvalho, 2013). The RMA algorithm (Irizarry et al., 2003), which performs background subtraction, quantile normalization, and summarization via median polish, was applied to the raw data of expression arrays to obtain the log 2 -normalized gene expression levels. All arrays for pericycle and cortex were normalized together, and arrays were assessed for quality implemented by the Bioconductor package arrayQualityMetrics (Kauffmann et al., 2009). An object of class ExpressionSet, which was generated by the oligo package, was inputted to the arrayQualityMetrics package and then a utility report created. To generate a robust total data set, we removed outlier arrays based on the between array distances (using the sum of distances S a as the quality metric; with a cutoff of 337), box plot (using Kolmogorov-Smirnov statistic K a as the quality metric; with a cutoff of 0.063), and MA plot (using Hoeffding's statistic D a as the quality metric; with a cutoff of 0.15). This resulted in a final set containing 236 comparable quality arrays (110 for cortex and 126 for pericycle) over the 6 treatments and 14 time points, i.e., an average replication of 2.95 per time point with at least two replicates per time point.

Confirmation of Time-Series Data Robustness and Data Processing
The expression of routinely used housekeeping genes (CBP20, Actin-2/8, UBC, YLS8, TUB4, and GAPDH, EF-1a) was analyzed. All the housekeeping genes had constant expression and were invariant to treatment/ tissue type, changing less than 1-fold in either time series (Supplemental Figure 5). A full expression data set for all three of our time series is available to view at the Arabidopsis eFP browser: http://bar.utoronto.ca/ efp_arabidopsis/cgi-bin/efpWeb.cgi?dataSource=Lateral_Root_Initiation Initiation (Winter et al., 2007). For analysis of differential expression, the variation between gene model=transcript expression was investigated by examining the 1547 transcripts, representing 711 genes that were determined to be differentially expressed within (using the software BATS [Bayesian Analysis of Time Series]; see below) one or more time series presented in this article (14% of all DE transcripts). Expression was determined to be identical between transcripts for 34%, similar for 36%, and different for 30% (Supplemental Data Set 9). Identical expression could indicate redundancy between transcript expression levels, with divergence indicative of functional variation, to be explored in future work. Each transcript was considered separately for the purpose of analyzing expression patterns (in clustering) and for comparing the total size of transcriptional response. However, for analysis of functional categories/GO terms represented in responses, and for network inference, analysis was performed at the gene level, to avoid a bias introduced by including the same gene more than once. For network inference, this meant that the average expression value of a transcript set was used; we considered this to be justified since 70% of transcripts showed similar or identical DE.
For analysis of gene families, we determined which transcripts belonged to which gene families according to the same TAIR gene annotation as used for expression array analysis (file: TAIR_gene_families_sep_29_09_update), and we analyzed gene families with more than 20 members. Two-tailed t tests assuming equal variance were used to compare trait values for wild-type and mutant seedling roots and trait values for Col-0 with treatment compared with mock treatment.

Determination and Clustering of DE Genes
DE genes within each time series were obtained using the software BATS (Angelini et al., 2008) and determined to be robust using the independent method of gradient analysis (Breeze et al., 2011). The BATS input file for each time series contains the rescaled log 2 gene expression values such that the vector of log 2 expression values of each gene had mean zero and variance one. Genes were considered to be DE if their Bayes factors in the BATS output file were less than a threshold (log 2 Bayes factor >3), which was determined by the histogram of log 10 of Bayes factors and the regression plots of gene expression levels from BATS. DE genes between treatments in each cell type were obtained using Gaussian process modeling implemented in the software GP2S (Stegle et al., 2010). In the GP2S input file, the log 2 expression levels of each gene for the two time series were rescaled such that their mean was zero and their variance was one. GP2S assigned each gene a Bayes factor that equals the difference between the likelihood of the expression level of a gene in two treatments being sampled from different Gaussian processes and that from the same Gaussian process.
Immediate responding genes were defined as those with |x 1 -x 0 | > f 1 x sd(x t ) for gene transcription time series x t , i.e., the absolute value of the change in the first time point being larger than f 1 times the time series SD. Furthermore, fast and downstream responding genes were defined as genes that were not immediate responders but had a dramatic change within 2 or 6 h of treatment (sd(x 0..2 ) > f 2 x sd(x t ), sd(x 0..6 ) > f 6 x sd(x t )). The remaining genes not in any of these three classes were designated late. Clusters were further defined as upregulated or downregulated if their expression increased or decreased at this point of differential expression. The thresholds f 1 , f 2 , and f 6 were chosen and then clusters were visually inspected to confirm separation of the clusters into the appropriate classes.
Hierarchical clustering of the differentially expressed genes was performed using SplineCluster (Heard et al., 2006). The mean expression levels of biological replicates of a DE gene at each time point in a time series was used as input. A reallocation function was implemented to reallocate outliers of each cluster into other more appropriate clusters at each agglomerative step. A prior precision value was finally determined after trying different values and comparing their effects on clusters (Supplemental Figure 6). SE of the mean was plotted on all graphs of gene expression.

Statistical Analysis of GO Terms
GOrilla (Eden et al., 2009) with hypergeometric test correction for clusters or BioMaps  with Fisher exact test with FDR correction was used to determine which GO categories are statistically overrepresented for groups of DE genes. The whole TAIR-10 annotated Arabidopsis genome was used as the reference set. To consider a cluster to have a significantly overrepresented GO term, we considered clusters/ groups with $15 genes with the GO term in $3 genes and a corrected P # 0.05 (results are tabulated in Supplemental Data Set 4; generic GO terms were removed from the table during analysis). To assess the degree of functional overlap between lists of differentially expressed genes, we obtained GO terms for every gene from the GO website for Arabidopsis using the TAIR10 annotation (Berardini et al., 2004). Since every gene could be associated with more than one function, we assigned as few GO terms as possible for a gene using the R package Go.db and by removing terms that could be associated via a directed acyclic graph as an ancestor to other GO terms assigned for the same gene (Carlson, 2016). This approach determined GO terms from nonoverlapping paths that end with a GO term, which is deepest in its respective path. Sometimes there were GO terms with less significant branching but even these were considered unique for that particular gene. As a final step, we sorted the GO terms across a particular treatment condition in order to determine unique molecular functions.

Gene Regulatory Network Inference and Analysis
The Bioconductor package GRENITS (Morrissey et al., 2012;Wang et al., 2014) was used for gene regulatory network inference using the mean expression levels of biological replicates (as outlined above) for each DE gene at each time point from 0 to 16 h. We assigned 2460 genes (2946 transcripts) to be putative regulators due to known or putative transcription factor activity based on data from the NCBI Conserved Domains database indicating the presence of conserved protein domains indicative of DNA binding combined with gene annotations (GO and TAIR), indicative of regulatory effect (Supplemental Data Set 2). Genes that were DE within time series (as calculated using BATS) were used as input, then genes DE when comparing treatment and untreated control samples (GP2S) were identified within those networks (Supplemental Files 1 to 14 and Supplemental Table 1). GRENITS gave the posterior probability of each directed link between two genes. A link probability threshold was chosen based on the confidence required for assigning links, then a link in the gene regulatory network was assigned if its posterior probability was greater than or equal to the link probability threshold. To identify network modules, the out-degree for each transcription factor in the gene regulatory network was counted. For each module, the expression level of the most central transcription factor was plotted with the one time point left-shifted expression levels of the upregulated targets, and the mirror symmetry of the one time point leftshifted expression levels of the downregulated targets. We confirmed that this formed a tight cluster with high correlations compared with randomly selected genes. All networks connected a similar proportion of the DE genes (37-48%) and networks had 6.3 targets per TF on average (Supplemental Data Set 6). For "regulatory backbone" networks, the edge-to-node ratio was similar for all networks, with 1.2 to 1.4 edges per node (Supplemental Data Set 6).

Network Edge Testing Using Promoter Motif Analysis
Data on the presence of known cis-acting TF binding sites in the 3-kb upstream region (likely promoter) were obtained from the VirtualPlant platform  for a test set of edges between nodes that are common between all interaction networks in a cell type. Promoter searching based on the web tool MEME LaB (Brown et al., 2013) was used to match transcription factors to targets in clusters. To query a set of specific consensus motifs from Franco-Zorrilla et al. (2014), regions upstream and downstream relative to either the transcription start or termination site for 500 and 1000 bp for each gene locus in the TAIR10 Arabidopsis genome annotation (chr1-5, C and M) from TAIR10 were obtained (ftp://ftp.arabidopsis.org/home/tair/Sequences/blast_datasets/ TAIR10_blastsets/). The pattern-matching PatMatch algorithm from the TAIR10 website was used, following Franco-Zorrilla et al. (2014) with a search on the given strand only, but using both the given and the complementary sequence (to retain the orientation of the motif). Out of the 26 edges predicted by in cortex networks, 10 were confirmed as representing known TF binding sites. The same method was applied to pericycle networks, where 20 out of 66 edges (30%) were confirmed. The predicted networks consistently outperform randomized networks (with a similar topology and node/edge distribution) by a factor of 2 (2-fold) in terms of being validated using known TF cis-element binding site data. LHY and CCA1 both have a high number of regulator links in the networks in this study, and almost all of these links are confirmed in VirtualPlant (;80% for LHY and ;95% for CCA1).
Predicted edges were compared with the interactions found by O'Malley et al. (2016) using ampDAP-seq and DAP-seq where possible. Across all networks in this study, the targets of 86 TFs had been identified using ampDAPseq and the targets of 145 TFs had been identified using DAP-seq (a union of 155 TFs). Of 1025 edges from TFs tested using ampDAP-seq, and 1601 edges from TFs tested using DAP-seq, 26 and 22% were confirmed, respectively. As a union, 28% of 1668 edges were validated by ampDAP-seq and/or DAP-seq. The percentage of validated edges was similar in each of the six networks (CU/ CN/CR/PU/PN/PR); validation information is in Supplemental Data Set 5, and validated edges are colored green in the full network Cytoscape session files (Supplemental Files 1 to 6).

Accession Numbers
Sequence data from this article can be found in the GenBank/EMBL data libraries under accession numbers At2g43500 (NLP8) and At2g23320 (WRKY15). All raw and processed microarray data have been deposited in the Gene Expression Omnibus (GSE91379).  Supplemental Data Set 4. GO term overrepresentation analysis.

Supplemental
Supplemental Data Set 5. Network edge designations with validation information.
Supplemental Data Set 6. Network statistics and TF module annotation and characteristics.
Supplemental Data Set 7. Glutathione S-transferase gene family expression is highly variable.
Supplemental Data Set 8. Gene family regulation across cell types and treatments.
Supplemental Data Set 9. Gene model analysis.