Inference of the Arabidopsis lateral root gene regulatory network suggests a bifurcation mechanism that defines primordia flanking and central zones

A large number of genes involved in lateral root (LR) organogenesis have been identi ﬁ ed over the last decade using forward and reverse genetic approaches in Arabidopsis thaliana . Nevertheless, how these genes interact to form a LR regulatory network largely remains to be elucidated. In this study, we developed a time-delay correlation algorithm (TDCor) to infer the gene regulatory network (GRN) controlling LR primordium initiation and patterning in Arabidopsis from a time-series transcriptomic data set. The predicted network topology links the very early-activated genes involved in LR initiation to later expressed cell identity markers through a multistep genetic cascade exhibiting both positive and negative feedback loops. The predictions were tested for the key transcriptional regulator AUXIN RESPONSE FACTOR7 node, and over 70% of its targets were validated experimentally. Intriguingly, the predicted GRN revealed a mutual inhibition between the ARF7 and ARF5 modules that would control an early bifurcation between two cell fates. Analyses of the expression pattern of ARF7 and ARF5 targets suggest that this patterning mechanism controls ﬂ anking and central zone speci ﬁ cation in Arabidopsis LR primordia.


INTRODUCTION
The root system of land plants is essential for anchorage plus water and nutrient acquisition. The ability of the root system to fulfill these functions is highly dependent on its architecture (Den Herder et al., 2010). Root system architecture is defined by three parameters: root growth rate, root branching rate, and root tropisms. Root system architecture plasticity determines the ability of plants to adapt to soil heterogeneity and various abiotic stresses, such as nutrient limitation or drought (Malamy, 2005). Lateral root formation in particular is a highly plastic and environmentally responsive trait that plays a central role in plant adaptation to different soils (Desnos, 2008;Lavenus et al., 2013).
In Arabidopsis thaliana, lateral roots form from pairs of pericycle cells called lateral root founder cells located in front of the xylem poles (Malamy and Benfey, 1997). Lateral root founder cells acquire their identity near the root tip in a region called the basal meristem or oscillation zone (De Smet et al., 2007;Moreno-Risueno et al., 2010;Van Norman et al., 2013). Priming is followed by the coordinated migration of founder cell nuclei toward the common cell walls and asymmetric anticlinal division (De Smet et al., 2007;De Rybel et al., 2010;Goh et al., 2012) leading to a stage I lateral root primordium (LRP). A two-layered LRP is then formed by a switch in the orientation of division planes from anticlinal to periclinal (Stage II). Additional divisions generate a multilayered dome-shaped primordium (Stage III to VII) that protrudes into the external tissues (endodermis, cortex, and epidermis) and shows a recognizable root meristem organization from stage VI onwards (Malamy and Benfey, 1997;Lucas et al., 2013). After emergence, the meristem is activated and the lateral root (LR) starts elongating. Many aspects of lateral root formation from priming to emergence are under the control of the phytohormone auxin (Lavenus et al., 2013). Two auxin-signaling pathways have been described. The ABP1 pathway has been proposed to mediate nontranscriptional auxin response through Rho of Plant GTPase signaling and transmembrane receptor-like kinases of the TMK family (Robert et al., 2010;Chen et al., 2012;Lin et al., 2012;Xu et al., 2014), while the nuclear TIR1/AFB pathway alters gene transcription by regulating the degradation of the Aux/IAA transcriptional repressors (Dharmasiri et al., 2005;Kepinski and Leyser, 2005). Auxin signaling modules have been defined as groups of coexpressed and strongly interacting Aux/IAA and ARF proteins that regulate subsets of auxin responsive genes (De Rybel et al., 2010). Priming was shown to be under the control of an IAA28dependent auxin signaling module, the only known target of which is GATA23, a transcription factor encoding gene triggering lateral root founder cell identity acquisition (De Rybel et al., 2010). The polarization of the founder cells prior to the first asymmetric division is regulated by another auxin signaling module featuring IAA14/SLR, ARF7, and ARF19 (Fukaki et al., 2005;Okushima et al., 2005aOkushima et al., , 2007Wilmoth et al., 2005;Lee et al., 2009;Goh et al., 2012). The known targets of this module encode LBD transcription factors (LBD16, 18, 29, and 33) as well as ARF19 itself (Okushima et al., 2005a(Okushima et al., , 2007. Besides its role in LR founder cell polarization, the IAA14/SLR module is also involved in LR initiation and LRP patterning together with a third module composed of IAA12/BDL and ARF5/MP (Vanneste et al., 2005;. One of the downstream targets of the IAA14/SLR module involved in LRP patterning is the PUCHI transcription factor (Okushima et al., 2005a;Vanneste et al., 2005;Hirota et al., 2007). So far, no target has been reported for the IAA12/BDL module.
While several regulators of lateral root formation and some of their targets have been identified in Arabidopsis, our knowledge of how these and other genes interact within a complex gene regulatory network (GRN) to control root branching is very limited. A powerful way to generate such knowledge is to make use of GRN inference techniques that reconstruct network topologies from transcriptomic data sets (De Smet and Marchal, 2010;Marbach et al., 2010). These techniques are based on a fairly simple principle: If two (or more) genes interact with each other (i.e., one is regulating the other), their expression levels should be linked to each other statistically. Therefore, it is theoretically possible to identify pairs or groups of interacting genes by analyzing statistical dependencies of their respective expression levels. Arguably, time-series expression data sets represent the best type of data for GRN inference modeling, as they permit the detection of delays between gene profiles, thereby revealing the directionality of the interactions. In the case of positive regulatory interaction, if gene X directly regulates gene Y, the expression profile of Y should be highly similar to the expression profile of X but shifted forward in time by m time units ( Figure  1C). This interval reflects the time needed to synthesize protein X and for protein X to bind to the promoter of gene Y. In the case of a negative regulatory interaction, the expression profile of Y is expected to be highly similar to the inverted expression profile of X shifted forward in time by m time units ( Figure 1C). Hence, by analyzing statistical linkage between gene expression profiles, it should be possible in principle to infer regulatory relationships, and it has already been successfully done in diverse biological systems (De Smet and Marchal, 2010;Wang and Huang, 2014).
Many GRN inference algorithms have been developed over the last decade, based on different formalisms and having different specificities and abilities (De Smet and Marchal, 2010;Marbach et al., 2010;Wang and Huang, 2014). In their comparative study as part of the DREAM3 challenge, Marbach et al. (2010) pointed out important limitations of the GRN inference approach. The authors reported that all algorithms, including the best performers, made systematic errors to various extents. These were classified into three categories: cascade errors, resulting in the introduction of shortcuts in genetic cascades; fanout errors, leading to the misinterpretation of coregulation as regulation; and fan-in errors, resulting from a failure to predict combinatorial interactions. Cascade errors and fan-out errors create false positives, whereas the fan-in errors generate false negatives (Marbach et al., 2010).
In this study, we created an algorithm called TDCor (for timedelay correlation) to infer the GRN that controls lateral root formation in Arabidopsis. We developed an approach based on GRN motif modeling to limit the number of the systematic errors made by the algorithm. TDCor was run on a subset of ;100 genes known to be involved in auxin biosynthesis, transport, or signaling pathways, plus lateral root initiation, patterning, or root apical meristem (RAM) maintenance. The predicted network topology links the very early activated genes involved in LR initiation down to the late activated cell identity markers through a multistep genetic cascade exhibiting both positive and negative feedback loops. We experimentally validated TDCor predictions for the key regulator ARF7 using an ARF7-GR transcriptomic data set and chromatin immunoprecipitation-PCR. Interestingly, TDCor predicts that the selected genes cluster in two major groups. The first group is predicted to be positively regulated by ARF7 and ARF19 and to repress the second group. The second group is predicted to be positively regulated by MP, ARF6, and ARF8 and to repress the first group. Using various reporter lines, we observed a spatial separation of two nonoverlapping spatial expression domains matching the two groups of genes, supporting the existence of two mutually exclusive sets of genes generating two zones with different identities in the LRP.

Statistical Analyses Can Predict Linear Interactions between Two Genes from Time-Series Transcriptome Data in a Complex Tissue
A GRN inference approach using time-series transcriptome data might be used to infer the GRN regulating LR development in the model plant Arabidopsis. However, the models underlying the currently available GRN inference algorithm make the implicit assumption that the data on which they are based were generated from homogenous cell populations. Because roots are made of heterogeneous cell types, we first analyzed whether this heterogeneity could pose a problem for GRN inference. Using an analytical approach (Supplemental Methods 1), we showed that only linear interactions can be inferred from timeseries data set generated from heterogeneous cell populations. We established that in a single cell or in a homogeneous cell population, under the conditions listed in Supplemental Table 1, if gene X regulates gene Y, the amount of gene Y transcript present at time t [which we call YðtÞ], would be linked to the amount of gene X transcript present at time t 2 d with d> 0 [which we call X(t 2 d]) by the differential equation (E). In this equation, the function f models the effect of protein X on the transcription rate of gene Y.
∃f : ℝ þ ↦ ℝ and a; g; d > 0; ∀t∈ ℝ; The different terms of the equation are explained in detail in Supplemental Table 2. When the g parameter is high, meaning that the degradation rate of the Y transcript is high enough, the system reaches quasi-steady state. Hence, ∀t ∈ ℝ; 0 a þ fðXðt 2 dÞÞ 2 g: YðtÞ⇒YðtÞ 1 g :ða þ fðXðt 2 dÞÞÞ Therefore, a strong relationship exists between the amount of transcript of a regulator gene X at time t minus a delay and the amount of transcript of the target gene Y. This implies that in this case, one can infer the existence of the interaction between X and Y by looking for a strong statistical dependency (linear or nonlinear depending on f) between the two variables Xðt 2 dÞ and YðtÞ. More rigorously, one can say that in the absence of (A) Changing the gravity vector (turning the Petri dish by 90°= applying the gravistimulus) causes the root to reorient downwards and creates a bend where a lateral root systematically initiates.
(B) The initiation and development of the lateral root in the bend follows a tightly reproducible timing. Each line on the plot shows the distribution of lateral root primordia stages (Malamy and Benfey, 1997) at a given time after the gravistimulus was applied. At 12 hag, >90% of LRPs have not initiated yet (preinitiation [PI] stage). Initiation occurs between 15 and 18 hag (stage 1), and then the LRP passes from one stage to the next one, approximately every 3 h. The E0, E1,., E4 stages correspond to successive stages after emergence. At E0, the LRP tip has reached the surface of the parental root; at E1, it is approximately one epidermal cell width from the surface of the parental root; while at E2, it is two epidermal cell widths, etc. The FE stage corresponds to a fully emerged and elongating lateral root.
(C) Gene regulatory network inference algorithm often make the assumption that if a gene X regulates a target Y positively, Y expression profile is similar to X expression profile but shifted forward in time. Conversely, if X regulates Z negatively, the profile of Z is expected to resemble the inverted profile of X. Here, we show that this principle stands in the case of heterogeneous cell populations when the interaction is linear (Supplemental Methods 1).
(D) Two expression profiles from the LR data set (of the SHR transcription factor and its direct target SCR), which follows the rule illustrated in (C).
statistical dependency between Xðt 2 dÞ and YðtÞ, the null hypothesis that X regulates Y following the assumptions shown in Supplemental Table 1 can be rejected. We next considered a heterogeneous cell population made of two homogeneous groups of cells in constant proportions, which could be for instance a very simple nongrowing and nondeveloping organ made of only two cell types. We showed that in this case, an ordinary differential equation (ODE) identical to (E) links the expression level of genes X and Y at the organ level when f is linear but not when f is nonlinear. No differential equation exists that can link X and Y transcript levels when f is nonlinear. This means that when X regulates Y in a linear manner, the possible existence of the interaction can be inferred (H0 not rejected) from transcriptomic data obtained from the whole organ. In contrast, when the interaction is not linear, the expression level of the target is not linked to the expression level of the regulator at the organ scale. Hence, it is not possible to infer the existence of nonlinear interactions by using transcriptomic data generated from a complex organ.
To illustrate this point, we simulated an interaction between two genes in an organ made of two homogeneous groups of cells (Figures 2A and 2B) and chose two different forms for f (one linear and one nonlinear). As predicted, we observed a strong relationship in the individual groups of cells in both cases ( Figures 2C and  2D, blue curves), although this relationship was linear in the first case while non-linear in the second case. However, at the wholeorgan scale, a relationship existed between transcripts level only when f was linear ( Figures 2C and 2D, blue curves). These results were then generalized to the case where the organ is composed of n homogeneous groups of cells in changing proportions (Supplemental Methods 1). In brief, GRN inference from timeseries transcriptomic data obtained from heterogeneous cell populations like roots is possible, but limited to linear interactions, or more realistically nearly linear interactions.
TDCor, an Algorithm Based on Time-Delay Correlation to Infer Linear Regulatory Interactions from Time-Series Transcriptomic Data Taking advantage of our ability to synchronize Arabidopsis LR development following a root gravitropic stimulus (Ditengou et al., 2008;Laskowski et al., 2008;Lucas et al., 2008;Péret et al., 2012), a time-series transcriptomic data set (herein after referred to as the LR data set) was generated encompassing every LR developmental stage from preinitiation to postemergence (Voß et al., 2015;Figures 1A and 1B). In our experimental conditions, stage I LR primordia are detected in Columbia-0 (Col-0) ;15 h after gravistimulation (hag) and LRs emerged around 42 hag. Root bends were microdissected every 3 h from 6 to 54 hag and RNA extracted to cover the entire process of LR development. At time point 0, a mature root segment located between the bend and the shoot was harvested at 9 h after an inductive gravitropic stimulation.
An initial study of well-characterized regulators and their targets within the LR data set revealed that, as expected, the expression profile of the target is often highly similar to the expression profile of its regulator and shifted in time. For example, the target of the SHORT-ROOT (SHR) transcriptional regulator termed SCARECROW (SCR) (Di Laurenzio et al., 1996;Helariutta et al., 2000;Sozzani et al., 2010) closely followed the dynamic changes in SHR mRNA abundance within the LR data set with a time delay estimated around 2 h ( Figure 1D). Similar relationships were identified for other well-known regulators and their direct targets, such as ARF19 and LBD16 (Okushima et al., 2005a(Okushima et al., , 2007. Hence, we concluded that a GRN inference approach based on the LR data set could potentially be employed to reconstruct the GRN controlling LR formation. We developed the TDCor algorithm specifically to identify linear interactions, as only linear (as opposed to nonlinear) interactions may be inferred from the LR data set (as described above). A detailed description of TDCor can be found in Supplemental Methods 2.
TDCor reconstructs the network topology from a subset of expression profiles (typically around 120 genes) in four main steps ( Figure 3). First, TDCor requires data preprocessing from the LR data set to obtain expression profiles of the genes selected. For each gene, an exact interpolating cubic spline function of the normalized expression profile is generated. Second, TDCor builds a list of plausible interactions. This is achieved by performing an extensive pairwise and bidirectional testing of whether a sufficient linear model with time delay exists between the expression (A) Schematic representation of the system that was modeled. We considered a simple noncombinatorial interaction between a regulator X and its positive target Y occurring in an organ made of two homogeneous groups of cells in equal proportion. The interaction was modeled using the ODE in equation (E). (B) Two profiles were taken from the LR data set to be used as template for gene X expression in the two groups of cells. The transcript accumulation profile of gene X in the whole organ was obtained by averaging the transcript accumulation profiles in the two groups of cells.
(C) The X-Y interaction was first modeled using a linear model (linear regulation function shown on the left of the panel). The profile of gene Y in each of the two homogeneous tissues was obtained from the corresponding profiles of gene X in these two tissues (shown in [B]) by solving equation (E) while taking the linear regulation function as function f. The profile of gene Y in the whole organ was obtained by averaging the computed profiles of gene Y in the two groups of cells. The three plots at the bottom show the relationship between the expression level of Y at time t and the expression level of its regulator X at time t 2 m, with m the time shift between the expression profiles of X and Y. m is the sum of the d parameter of the ODE and the relaxation time of the system (0.87 h), which depends on the g parameter. (D) The X-Y interaction was next modeled using a nonlinear model (nonlinear regulation function of Mikaelis and Menten type shown on the left). The profile of gene Y in the two homogeneous tissues was obtained from the profiles of gene X in these two tissues (B) by solving equation (E) and taking the nonlinear regulation function as function f shown on the left of the panel. The profile of gene Y in the whole organ was obtained by averaging the computed profiles of gene Y in the two groups of cells. The three plots at the bottom show the relationship between the expression level of Y at time t and the expression level of its regulator X at time t 2 m, with m the time shift between the expression profiles of X and Y. m is the sum of the d parameter of the ODE and the relaxation time of the system (0.87 h), which depends on the g parameter. profiles of the selected genes. Because it is critical to calculate precisely the delay between the expression profiles for each pair, four different methods were combined (Supplemental Methods 2). Furthermore, to reach high precision in the calculations (i.e., below the sampling interval), TDCor uses exact interpolating cubic spline functions. To build the list of plausible interactions, TDCor makes use of prior knowledge concerning the nature of the genes used to reconstruct the network. Genes are classified into four categories: transcriptional activators, transcriptional repressors, transcriptional regulators (including activators or repressors or proteins with unknown effect on transcription), and nonregulators (e.g., cell wall remodeling enzyme). Using this minimal prior knowledge, some interactions can be ruled out. For instance, a nonregulator cannot directly regulate target gene expression. Hence, any interactions involving this type of gene as a regulator can be eliminated. Similarly, a transcriptional activator cannot negatively regulate its primary target. Hence, any negative regulatory interaction with a short delay involving a transcriptional activator is regarded as unlikely and subsequently removed.
In its third step, TDCor tests whether the linear statistical dependencies detected at step 1 between pairs of expression profiles are more likely to be due to true direct regulation or rather to be a consequence of coregulation or of indirect regulation. The preliminary network topology established at step 1 is pruned at step 2 using several filters. All the decisions made by the algorithm for whether to retain or not retain an interaction involve the comparison of a numerical predictor with one or two user-defined threshold parameters. Because the choice of the threshold values influences the final result and because the result of the pruning can depend on the order of node analysis, TDCor uses a double bootstrap protocol in order to estimate the sensitivity of the various edges in the network to changes in the parameter values and in the order of node analysis. The two steps described above are repeated N times (typically N>1000) with different parameters set independently and randomly chosen at the beginning of each iteration within a user-defined parameter space, and each pruning procedure is repeated n times (typically n 5 10) with the same parameters. This produces N:n network topologies that are then merged. The bootstrap value for each edge is obtained by counting the number of predicted network topologies in which it occurs. At step four, the n R best regulators (typically n R = 4) showing the highest bootstrap scores are selected for each node.

Comparison of TDCor Performance versus Other GRN Inference Algorithms
We inferred a LR GRN from the LR transcriptomic data set using TDCor and two existing GRN inference algorithms capable of reconstructing network topologies from time-series transcriptomic data, namely, VBSSM and TDARACNE (Beal et al., 2005;Zoppoli et al., 2010). A set of 120 genes was selected on their known or suspected function in LR development. This includes genes from the auxin pathway (metabolism, transport, or signaling pathways), genes involved in root apical meristem patterning and maintenance, or genes that play a role in lateral root initiation, patterning, or emergence (Supplemental Table 3). The broad span of functions for these genes enables one to cover major aspects of lateral root formation and therefore obtain a network that is representative of the entire LR GRN.
GRNs predicted by VBSSM with default parameters and with three different cutoff levels of significance (Z = 1.65, Z = 2.33, and Z = 3) are shown in Supplemental Table 4. The network obtained using the highest significant threshold (Z = 3) features 56 nodes and 82 edges, whereas the one obtained with the lowest stringency has 120 nodes and 651 edges. The GRNs predicted by TDARACNE with low and high stringency parameter sets are shown in Supplemental Table 5. The two TDARACNE networks feature 291 (89) and 582 (88) edges (nodes), respectively. The TDCor algorithm was run with low stringency parameters (Supplemental Table 6), and only interactions with a bootstrap value greater than 10% were kept in the network. The network topology predicted by TDCor with these parameters features 122 nodes and 358 edges corresponding to 206 positive interactions and 152 negative interactions (Supplemental Figure 1).
To assess the relative performance of the three algorithms on the lateral root data set, we determined how many of 23 already TDCor uses a double bootstrap protocol. The parameter values are randomly chosen within a user-defined parameter space at the beginning of each big iterative loop. Node order analysis is randomized before each pruning filter and the last part of the pruning procedure is itself repeated several times with the same parameter values (but different order of nodes analysis). MRST, master regulator and signal transducer, i.e., a gene that could be posttranslationally activated in response to the stimulus and could start the whole genetic cascade. TPI, triangle pruning index; DPI, diamond pruning index. These two indices are numerical topology predictors for respectively some specific three-and four-gene motifs (Supplemental Methods 2). experimentally validated regulatory interactions were predicted (Supplemental Table 7). A surprising low number of interactions was predicted by VBSSM and TDARACNE. VBSSM was able to correctly predict only one interaction (BBM self-regulation) out of 23. In most cases, these direct interactions were misinterpreted by VBSSM as coregulation of low significance. TDARACNE consistently predicted five interactions out of 23 (22%), while seven others were misinterpreted as coregulation in the topology obtained with low stringency parameters (Supplemental Table 7). This suggests that the LR data set was not suited for these algorithms. On the other hand, TDCor was able to correctly predict 14 out of the 23 interactions (60%) described in the literature, flagged another one as coregulation, and failed to predict the remaining eight (Supplemental Table 7). We concluded that TDCor is suitable for analyzing the LR transcriptomic data set to predict new regulatory interactions with a high level of confidence.

Experimental Validation of TDCor Predictions for ARF7 Targets
In order to further evaluate the performance of TDCor, we experimentally tested its predictions for ARF7 targets. ARF7 encodes a key transcriptional activator required to trigger lateral root initiation (Okushima et al., 2005a(Okushima et al., , 2007Moreno-Risueno et al., 2010). To determine ARF7 primary and secondary targets in roots, the arf7 arf19 double mutant was complemented by expressing a dexamethasone (DEX) inducible ARF7-GLUTICORTICOID RECEPTOR (GR) fusion protein under the native ARF7 promoter. The PRO ARF7 :ARF7-GR arf7 arf19 line was treated for 4 h either with DMSO (mock control), the synthetic auxin naphthalene acetic acid (NAA), cycloheximide (CHX), DEX, or with one of the four possible combined treatments. The whole seedling root was then harvested for transcriptomic analyses. The global response of the plant to the various treatments was formally decomposed into five independent response components that may or may not be activated depending on the treatments (Supplemental Table 8). These components include the ARF7-dependent primary and secondary auxin responses and the ARF7-independent primary and secondary auxin responses. Because CHX can have a strong impact by itself on transcript accumulation (not directly related to it blocking translation), a CHX response control was taken into account.
In order to quantitatively estimate the ARF7-dependent primary and secondary responses from the data, a model in the form of a system of seven equations with seven unknowns was solved for all genes (Supplemental Table 9 and Supplemental Methods 3). In these equations, the two components of interest were represented respectively by the a and b unknowns, which are interpreted as fold changes. Primary (respectively secondary) positive targets of ARF7 should show a > 1 ( b > 1), while primary (respectively secondary) negative targets of ARF7 should show 0 < a < 1 ( 0 < b < 1). Conversely, genes that are not primarily regulated by ARF7 are expected to show an a parameter close to 1. For our analysis, we regarded a gene as regulated through a given pathway if the parameter corresponding to this response component was either greater than 1.5 (upregulation) or lower than 1/1.5 (downregulation).
TDCor uses an index of directness in order to predict whether edges in the network are likely to correspond to direct or indirect interactions. In the network topology, 10 genes were located one edge downstream of ARF7, five of which are predicted to correspond to direct interactions (Figures 4A; Supplemental Figure 2). All five were validated as ARF7 primary targets using the ARF7-GR data set, although ARF4 is potentially both a primary and secondary target of ARF7 ( Figure 4B). Consistently these genes display between one and four canonical ARF motifs termed AuxRE in their promoters (Supplemental Table 10). Furthermore, binding of ARF7 to the target promoters was confirmed by chromatin immunoprecipitation-quantitative PCR (ChIP-qPCR) for ARF19, LBD16, LBD29, and IAA19 (Supplemental Table 10; Figure 4C). TDCor also correctly predicted the previously documented double regulation of LBD16 by ARF7 and ARF19 (feed-forward loop). GATA23, LBD17, PIN7, IAA11, and AFB3 are located one edge downstream of ARF7, but their index of directness predicts they are likely to be indirect targets of ARF7 ( Figure 4A). GATA23, LBD17, PIN7, and IAA11 are validated by the ARF7-GR data as positive secondary targets of ARF7 ( Figure 4B). AFB3 was not validated as negative secondary target of ARF7 ( Figure 4B). The estimated delay between ARF7 and AFB3 is 6 h, which would be too long to test using the ARF7-GR data (which was harvested 4 h after treatments). In total, eight out of the 10 genes located one edge downstream of (A) Topology of the network predicted by TDCor corresponding to genes located up to two edges downstream of ARF7 (bootstrap > 10%). The TMO5 gene is predicted to be a negative target of ARF4 and to downregulate several of ARF7 targets. For clarity, this node was removed from the network. TDCor also predicted that the activity of the ARF7 protein could be regulated posttranslationally upon gravistimulus application to start a transcriptional cascade (denoted by the blue "Stimulus" node). (B) ARF7-GR transcriptomic data analysis for the genes located one edge downstream of ARF7 in the predicted network topology. The a parameter corresponds to the ARF7-dependent primary auxin response. The b parameter corresponds to the ARF7-dependent secondary auxin response. a:b is the total ARF7-dependent auxin response. These parameters are equivalent to fold changes. The bar plots show the log 2 values of the parameters estimated from the ARF7-GR transcriptomic data for the various genes located one edge downstream of ARF7 in the network topology (information about the estimation method is provided in Supplemental Methods 3).
(C) Validation of the predictions of ARF7 primary targets by ChIP-qPCR. The schematic on the left indicates the position of the putative ARF binding sites termed auxin response element in the promoters of the predicted ARF7 primary target genes and the location of the analyzed amplicons (see Supplemental Table 10 for more details).
Eighteen genes are found two edges downstream of ARF7 in the network topology predicted by TDCor, 14 of which are predicted as positive targets of ARF7, while the remaining four are predicted as negative targets of ARF7 ( Figure 4A). The ARF7-GR data validated as positive secondary ARF7 targets nine of the 14 predicted positive targets, including two LBDs (LBD18/33), two auxin conjugating enzymes (GH3.1 /3.5), three Aux/IAAs (IAA29, SLR, and SHY2), and two auxin carriers (PIN3 and AUX1). The AP2 transcription factor PUCHI and the auxin conjugating enzyme GH3.3 are predicted as a positive target of LBD16 and/or LBD29 and consequently an indirect target of ARF7. However, the ARF7-GR suggests that these two genes are primary targets of ARF7. The four predicted positive targets that were not validated were AFB2, SHP1, KRP2, and ARF1. None of the four predicted negative secondary targets of ARF7 could be validated by the ARF7-GR data. ARF9 is an indirect (A) Whole-network topology predicted by TDCor. The layout shows two groups of genes featuring different sets of positive ARFs (yellow nodes) and negative ARFs (green nodes). Within each group, genes regulate each other positively (red edges), while each group regulates the other negatively (blue edges). SHP1 is the only exception to the rule as it is predicted to be positively regulated by members of the two groups. (B) Expression profile of PLT5 (black) and PHV (red) in the LR data set. PLT5 is predicted to be part of the group featuring ARF7-ARF19, while PHV is predicted to be part of the other group, which features ARF6-ARF8-MP. The high similarity that can be observed between the two profiles during the first 30 h suggests that PLT5 might be indirectly positively regulating PHV. (C) Expression profiles of MYB52 (black) and MYB56 (red) in the LR data set. These profiles could be compatible with these two genes acting as intermediary links between PLT5 and PHV.
Hence, comparison of the ARF7-GR data with TDCor predictions from the LR data set validated most of the predicted ARF7 targets. The genes that are the more strongly and specifically regulated by ARF7 in the ARF7-GR data set were correctly predicted as ARF7 targets. TDCor was also able to successfully differentiate between primary and secondary targets. Ninety percent (9/10) of the genes located one edge downstream of ARF7 could be validated as such, and for all of them but one, the predicted directness was correct. The validation rates of the genes located up to two edges downstream, which would correspond to two regulatory steps downstream or more, is very good though lower (20/28; 71%). It needs to should be stressed that the ARF7-GR data set was designed to detect targets regulated only 4 h after treatment, which may be sufficient for two regulatory steps. Overall, TDCor recovered >65% (21/31) of ARF7 targets.
Besides the ARF7 node, a number of predictions from TDCor were consistent with published data. For instance, most ARF7 targets are predicted as repressed by the chromatin-remodeling factor PICKLE (PKL) (Fukaki et al., 2006). Furthermore, most of the genes encoding cell cycle regulators (CYCA2;4, DPA, CYCB1;1, CDKB2.1, E2FA, CYCD6;1, and CYCD3;1) included in the analysis were predicted by TDCor as direct or indirect targets of transcription factors of the PLETHORA (PLT) family, namely, PLT1, PLT2, PLT3, PLT4/BABYBOOM (BBM), and AINTEGUMENT (ANT), as well as of the LONESOME HIGHWAY (LHW) regulator. This is fully consistent with the reported function of the PLT genes in the control of mitotic activity in the RAM, with the role of LHW in the regulation of the size of the vascular stem cells pool in the RAM, and also with the ascribed function of ANT in the control of shoot lateral organ size (Mizukami and Fischer, 2000;Aida et al., 2004;Nole-Wilson et al., 2005;Galinha et al., 2007;Ohashi-Ito and Bergmann, 2007;Krizek, 2009;Horstman et al., 2014;Ohashi-Ito et al., 2013). The negative ARF, ARF2, is predicted to repress indirectly cell cycle genes through indirect repression of PLT genes and other regulators. Consistently ARF2 is known as a repressor of cell division and organ growth (Okushima et al., 2005b;Schruff et al., 2006). Altogether, the LR development GRN predicted by TDCor could be validated by our analysis of ARF7 targets and is consistent with our current knowledge of the function of the genes involved.

TDCor Reveals That the LR Development GRN Is Organized in Two Regulatory Modules
Strikingly, the network predicted by TDCor was organized into two groups of genes ( Figure 5A). Each of these groups features different sets of positive and negative ARFs. On one hand, around 50 genes are predicted by TDCor to be downstream of ARF7 and ARF19 and form a first module, while the remaining genes are associated with MP, ARF6, and ARF8 to form a second module. The genes within the same module are predicted to positively regulate one another ( Figure 5A, red edges) to form a genetic cascade, although the cascade organization is less clear in the second group due to the presence of multiple positive feedback loops. In contrast, the two modules are predicted to interact through negative interactions ( Figure 5A, blue edges). We noted that the genes predicted to be at the bottom of the    ARF7/ARF19-activated cascade (e.g., PLT5 and PLT7) have relatively similar expression profiles to the genes predicted to be at the top of the second cascade (e.g., PHV and ANT), particularly during the first 35 h of the LR time series. However, the time shift between these genes (around 8 to 10 h) exceeds the TDCor upper limit ( Figure 5B). This suggested that the predicted network topology could be artificially cleaved into two parts simply because the genes establishing the link between the two blocks could be missing from our list. For this reason, we searched for genes that could potentially be acting as missing links between PLT5-PLT7 and ANT-PHV. We searched in a list of ;1060 transcription factors from the AtTFDB database (Davuluri et al., 2003;Palaniswamy et al., 2006) for genes displaying an expression pattern with high correlation (>0.8) and +3 h time delay with PLT5/PLT7 and high correlation (>0.8) with 23 h time delay with ANT/PHV in the LR data set (Table 1, Figure 5C). We also noted that PLT1 has many predicted targets but no regulators. However, TDCor delay measurements indicated that PLT1 was ;7 h downstream of ARF8, MP, and BBM. Therefore, we searched and identified 13 genes potentially acting as missing link between ARF8 and PLT1 ( Table 2). All of these genes were added to our original list and TDCor was rerun. The new network topology contained a single very long genetic cascade featuring ARF7 at the top and meristematic and patterning genes at the bottom ( Figure 6A). The "full cascade" topology is supported by the observations that (1) MP and PLT3 are now located several steps downstream of ARF7 (De Hofhuis et al., 2013), and (2) the MYB56 gene, one of the predicted links, has recently been shown to control seed size by regulating cell division in the integument of the seed coat, exactly as ANT, its predicted target (Mizukami and Fischer, 2000;Zhang et al., 2013). Interestingly, at the very top of the cascade, the gravistimulus is predicted by TDCor to trigger a nontranscriptional response leading to posttranslational activation of the ARF7 protein ( Figures 4A and 6A). The most likely hypothesis is that upon gravistimulation, auxin accumulates in the root bend, triggering the degradation of Aux/IAAs and subsequently the activation of the ARF7 protein.

The LR Development GRN Topology Creates Successive Waves of Transcription
Topological features of the inferred LR network can explain many aspects of the observed dynamics of the network. First, TDCor predicts that the genes upregulated several regulatory steps downstream of ARF7 would act as repressors of the early ARF7-activated genes, creating late negative feedback loops. This topology could explain why ARF7 targets and downstream target genes are transiently expressed, hence creating successive transcriptional waves ( Figure 6B). Two groups can be distinguished among the ARF7 targets. Genes such as ARF4, IAA11, GATA23, and LBD17 are very quickly repressed after initiation, while others such as LBD16 and PUCHI remain expressed in the developing LR primordium and are repressed more slowly (Supplemental Figure 2). The rapid repression of early ARF7 targets is predicted to be mediated by two negative ARFs (ARF16 and ARF18), PHV and MYB56, as soon as initiation has taken place around 15 hag ( Figure 7A). The group of late-repressed ARF7 targets is predicted to be repressed by genes that are activated a few hours later such as the HDIII-Zip PHB (starting around 20 h), the KANADIs (starting around 30 h), and PLT1, PLT2, PLT3, and PLT4/BBM (starting between 30 and 40 h). Most noticeably, the repression of ARF7 targets is predicted to involve epigenomic factors. In particular, TDCor predicts that the chromatin remodeler PKL/SUPPRESSOR OF SLR2 and the Polycomb-group protein CURLY LEAF (CLF) are upregulated by meristematic and patterning genes to repress ARF7 targets ( Figure 7C). Supporting this prediction, we found that the Gene Ontology categories related to DNA methylation (P values < 1e-67), histone methylation (and in particular H3K9 methylation; P values < 1e-90), and negative epigenetic regulation of gene expression (P values < 1e-61) are very significantly overrepresented among the ;3000 genes positively correlated with PKL, CLF, and other meristematic genes in the lateral root data set (Supplemental Figure 3). Many positive feedback loops are also predicted between meristematic genes ( Figure 7D). This could explain how, once activated, the meristematic genes remain expressed even though their regulators get repressed ( Figure 6B). In the absence of positive feedback loops, the negative feedbacks alone would generate oscillations of genes expression.

Mutual Inhibition between ARF7 Targets and Meristematic Genes Specifies LRP Central and Flanking Domains
One characteristic of the predicted network topology is the presence of incoherent feed-forward loops (Mangan and Alon, 2003;Alon, 2007). On the one hand, ARF7 and its targets are predicted to upregulate meristematic and patterning genes through a long genetic cascade ( Figure 6A). On the other hand, ARF7 targets are also predicted to repress these very same genes through a much shorter pathway ( Figure 7B). In particular, TDCor predicts that ARF7 quickly activates the expression of three negative ARFs, ARF4, ARF2, and ARF9, which would repress meristematic genes, such as PLT3, SHR, MP, and ARF8, and patterning genes, such as GLABRA2 (GL2), ENHANCER OF GLABRA3, or KAN4 ( Figure 7B). Hence, ARF7 targets and meristematic/patterning genes are predicted to mutually inhibit (A) A few genes identified as potential missing links between the two groups of genes shown in Figure 5A were added to the list of genes used to build the network and TDCor was rerun. As expected, TDCor predicted that these genes could be acting as intermediates between the two groups of genes, producing a gigantic genetic cascade headed by ARF7. For clarity, the negative interactions are not represented. Refer to Figure 2A for the key.
(B) Normalized expression profiles of some of the genes in the cascade. The red lines indicate the position of a primary peak of expression, which is noticeable mostly for the genes located upstream in the network. This peak moves as wave downstream along the network. The pink lines indicate the position of a secondary peak of expression that is mostly present in the profiles of the genes located downstream in the network. each other ( Figure 5A). The existence of the incoherent feedforward loop could explain why the genes downstream of PHV ( Figure 6A) are more slowly upregulated than the genes located more upstream in the cascade ( Figure 6B). Indeed, at the time when they get upregulated around 20 to 30 h after gravistimulation, the ARF2/4/9 repressors are expressed at a relatively high level. This mutual inhibition should lead to a bifurcation between two different fate decisions and the two groups of genes should be expressed in different nonoverlapping spatial domains.
In order to validate this prediction, we analyzed the spatial expression pattern of several genes predicted to be part of one or the other of the two groups. For the ARF7 targets, we analyzed the expression pattern of LBD16, PUCHI, and PLT5 genes, while MP and PLT1/2/3/4 genes were selected as meristematic genes. Using transcriptional and translational fusions with fluorescent proteins, we observed that the two groups of genes are expressed in separate spatial domain. PLT5, LBD16, and PUCHI were expressed in the flanks and at the base of LR primordia and only weakly expressed in the center of the primordium (Figure 8). After emergence, PLT5 is upregulated in the prospective central vasculature domain (Figure 8). In contrast, MP and the PLT1/2/3/4 genes were observed to be expressed in the center of the primordia and to be excluded from the flanks (Figure 9). Taken together, this suggests that the mutual inhibition between ARF7 targets and meristematic genes lead to the specification of two spatial domains in the LR primordium: the flank (where division stops) and the center (which will produce the new meristem with the new stem cell niche).

DISCUSSION
Over the last decade, rapid progress has been made in the genetic analysis of LR development (reviewed in Lavenus et al., 2013). Forward and reverse genetic screens have led to the identification of a number of key regulators of LR formation. Interestingly, LR development is not dependent on a stereotypical cell division pattern  but is characterized by conserved events of symmetry breaking that generate different cellular domains. Such complex coordinated cellular behaviors result from interactions between multiple genes within a GRN. Hence, in order to understand how a pair of pericycle cells generate a new root meristem showing a complex and yet robustly reproducible organization, it is necessary to (1) identify the topology of GRN that controls LR formation and (2) understand how the topology can generate the observed behavior.
Here, we performed GRN inference on a LR time-series transcriptomic data set (Voß et al., 2015). The LR data set was generated from tissues composed of multiple cell types rather than a homogenous cell population. Our analysis revealed that only linear interactions or more realistically close-to-linear interactions could be inferred from such data. We developed an algorithm called TDCor specifically designed to infer close-tolinear interactions from a time-series transcriptomic data such as the LR data set. While designing TDCor, we paid careful attention to two points. First, in order to obtain accurate estimates of the delay between expression profiles, which is crucial information for network reconstruction, we combined four methods of delay estimation. Second, to limit the number of systematic errors made by our algorithm, we developed novel filtering methods based on modeling of small network motifs. These filters use indices, the distributions of which depend on the topology of the network. TDCor pruning methods depend on the unknown distributions of biological/technical parameters, such as transcript degradation rates or the time needed for protein synthesis and maturation. However, the indices were made in such way that their interesting and useful properties should not depend on the shape of the parameter distribution but only on the topology. Moreover, we assumed that the biological parameters were distributed uniformly in a large interval. In particular, the time shift between expression profiles that could be obtained under our parameters distribution ranged from 40 min to 4 h. Thus, we chose the most extreme case producing the widest index distributions and, therefore, the softest pruning, which limits the risk of data overinterpretation (i.e., excessive pruning).
From the combined analysis of the TDCor predictions, the transcript profiles in the LR data set, and our ARF7-GR data, the following model can be proposed for the functioning of the LR GRN. Before the gravistimulus (t = 0 h), according to our transcriptomic data, very few genes related to LR organogenesis are expressed at high level in the root. One of them is ARF7. TDCor predicts that at t = 0, the ARF7 protein is in an inactive state and subsequently is activated during the first 6 h after gravistimulation. The most likely hypothesis is that ARF7 is bound to Aux/IAAs repressors. Upon gravistimulation (0 to 6 h), auxin redistribution (Ditengou et al., 2008) triggers Aux/IAAs degradation, enabling ARF7 to activate its targets. Among the latter targets are the auxin transporters PIN7, PIN3, PIN1, and AUX1 (according to TDCor predictions, our ARF7-GR data, and previous studies; Laskowski et al., 2006;Péret et al., 2013) whose induction would serve to reinforce auxin flux toward the future initiation site (positive feedback). Subsequently (6 to 12 h), major targets of ARF7 including LBD16, LBD29, or ARF19 reach their highest expression level and trigger LR initiation (Okushima (A) A late negative feedback involving noticeably the two negative ARFs, ARF16 and ARF18, could be involved in the repression of some of ARF7 targets. (B) ARF7 is predicted to upregulate a number of genes among which three negative ARFs (ARF2, ARF4, and ARF9) that in turn would repress the "meristematic and patterning genes" located more downstream in the cascade ( Figure 6A). (C) ARF7 targets are predicted to be repressed through epigenomic mechanisms involving the chromatin remodeling factor PKL and the Polycombgroup protein CLF. (D) Patterning and meristematic genes are predicted to form many positive feedback loops.  (Laplaze et al., 2005) at initiation, preemergence stage, and postemergence stage. (C) Expression pattern of PUCHI monitored by genomic PRO PUCHI :PUCHI-GFP (Hirota et al., 2007) at initiation, several preemergence stages, and postemergence stage.  predicted to repress their own repressors and our analysis suggests chromatin silencing mechanisms may play a central role in this process. The involvement of chromatin silencing as well as the presence of numerous predicted and known positive feedback loops between meristematic genes (Passarinho et al., 2008;Lau et al., 2011) could explain why, once activated, the meristematic genes would self-sustain and durably repress ARF7 targets. At the very bottom of the cascade TDCor predicts that PLT upregulates cell cycle genes as well as radial patterning genes such as genes involved in epidermis stem cells activity (FEZ and SMB) and in their daughter cell differentiation (KAN1, ATML1, WER, and GL2), hence completing the building of the new meristem. Mathematical models will certainly be most helpful for understanding in more detail the importance of the timing and the role of the successive transcriptional waves and predicted feedbacks in the LR GRN, as done previously in Escherichia coli, for instance (Zaslaver et al., 2004).
Another notable feature of the predicted network topology (beside its bimodularity) is the presence of three positive ARFs at the head of the meristematic genes: ARF6, ARF8, and MP. This suggests the existence of an auxin-dependent developmental checkpoint. It was previously demonstrated (1) that although most of the Aux/IAAs and positive ARF interact with each other in yeast (Vernoux et al., 2011), some Aux/IAA and ARFs are likely to interact more strongly with each other, thus defining auxin signaling modules (Weijers et al., 2005); and (2) that Aux/IAA proteins do not show the same sensitivity to auxin (Calderón Villalobos et al., 2012). In particular, IAA14/SLR, the Aux/IAA protein that associates with ARF7, is more sensitive to auxin than IAA12/BDL (Calderón Villalobos et al., 2012), the Aux/IAA protein associated with MP. In the apical part of the primordium where auxin concentration is higher, BDL would get degraded and MP-ARF6-ARF8 would activate the genetic program leading to meristematic identity acquisition and to ARF7 target repression in this region. However, cells in the flanks would not have enough auxin to activate the MP-ARF6-ARF8 module and therefore would not express the meristematic genes. This auxindependent switch could control the bifurcation between meristematic cell and flank cells based on the positional information provided by auxin concentration. This would explain why meristematic genes are expressed only in the center of the primordium where the concentration of auxin is the highest (Benková et al., 2003) and not in the flanks. Overall, we report a GRN topology that is able to (1) causally link the earliest molecular events to the latest events of LR formation and (2) could explain how cells acquire different identity (flanks versus meristematic) depending on their position in the LR primordium. The proposed model provides a novel conceptual framework for understanding lateral root primordium development and stem cell niche establishment.

Plant Material and Growth Conditions
Arabidopsis thaliana seeds were surface-sterilized for 6 min in a solution of 0.88% sodium dichloroisocyanurate (v/v) and 90% ethanol (v/v) and then washed three times in 90% ethanol (v/v) and dried under sterile air flow. Seeds were plated on 0.53 Murashige and Skoog solid medium containing 0.7% (v/v) plant agar and supplemented with B5 vitamins. Plates were cold treated at 4°C for 48 h to synchronize germination and then incubated in continuous light (22°C) in a nearly vertical position. The LR time-series transcriptomic data set was generated using wild type plants of Col-0 ecotype (Voß et al., 2015). The ARF7-GR transcriptomic data set was generated using PRO ARF7 :ARF7:GR nph4-1 arf19-1 plants (Okushima et al., 2007). The PRO LBD16 :GUS line was previously described (Laplaze et al., 2005). The PRO PUCHI :GFP-PUCHI line was published by Hirota et al. (2007). The PRO PLT :PLT-YFP lines were kindly provided by Ben Scheres (Hofhuis et al., 2013). For the ChIP-qPCR experiment, Col-0 and arf7-1 plants were used (Okushima et al., 2005a).

Root Culture Treatments and Fixation for Chromatin Immunoprecipitation
Root tissue was generated by adding 10 to 20 surface sterilized Arabidopsis seeds to flasks containing 100 mL root growth media (3.2 g/L Gamborg's B5 basal salt, 1 g/L MES hydrate, 20 g/L sucrose, and 1 mL/L Gamborg's B5 vitamin mix [10003], pH 5.8) and gently shaking in the dark for several weeks. Flasks containing ;5 g of root culture tissue were pretreated were with 1 mM NAA. Cross-linking of DNA and protein complexes was performed by submersing the tissue in 40 mL of fixation buffer (0.1 M sucrose, 50 mM NaCl, 10 mM KH 2 PO 4 , pH 7.0, 1% formaldehyde, and 10 mM MG132) and placing under vacuum for 20 min. Cross-linking was stopped by adding 1.25 M cold glycine to a final concentration of 0.125 M. The tissue was rinsed with water and flashfrozen in liquid nitrogen.

Lateral Root Data Set
To generate this data set, 3-d-old Col-0 plants were gravistimulated at t = 0 as described (Lucas et al., 2008). The root bends subsequently formed by gravity-induced downwards reorientation of the root apices were cut every 3 h from 6 to 54 hag and used for microarray analysis with the Affymetrix ATH1 array. For the first time point corresponding to t = 0, young unbent mature root material from the 9-h time point seedlings was used. For each time point, four independent biological replicates were generated (Voß et al., 2015).

VBSSM
The VBSSM algorithm was run in Matlab with the parameters indicated in Supplemental Table 11. Network predictions were visualized and analyzed using Cyctoscape software (Shannon et al., 2003;Cline et al., 2007;Smoot et al., 2011).

TDARACNE
The TDARACNE was obtained from Bioconductor and was run in R with the parameters indicated in Supplemental Table 12. Network predictions were visualized and analyzed using Cyctoscape software (Shannon et al., 2003;Cline et al., 2007;Smoot et al., 2011).

TDCor
The TDCor algorithm was run in R with the parameters indicated in Supplemental Table 13. This set of parameters, which leads to a "soft" pruning of the network, was determined by manual optimization. The performance of the algorithm was estimated based on interactions reported in the literature and plausibility of the predictions given the available knowledge (Supplemental Table 7). Network predictions were visualized and analyzed using Cyctoscape software (Shannon et al., 2003;Cline et al., 2007;Smoot et al., 2011). The code has been made available in the form of an R package called "TDCor" that can be downloaded from the CRAN repository (http://cran.r-project.org/) and directly installed in R version $ 3.1.2. In complement to the package documentation a step-bystep tutorial is provided in Supplemental Methods 4.

Missing Link Identification
The description of the method for identifying candidate genes potentially acting as missing links between PLT5/7 and ANT/PHV, and between ARF8 and PLT1, is provided in Supplemental Methods 5.

ARF7-GR Data Set
We used PRO ARF7 :ARF7:GR nph4-1 arf19-1 for DNA microarray analysis of downstream target of ARF7 (Okushima et al., 2007). The 5-d-old seedlings were transferred to the medium with or without the indicated chemicals (1 µM NAA, 2 µM DEX, and 10 µM CHX), and roots were harvested after 4 h of treatment. Microarray analyses were independently performed two times with the Arabidopsis ATH1 Genome Array (Affymetrix). Data analysis was performed by R software and Microsoft Excel (Supplemental Methods 3). The processed transcriptomic data for all genes on the chip is provided in Supplemental Data Set 1 in the tab entitled "Whole data set." The "TDCor validation" tab contains the processed data for the genes that were included in the network reconstruction. The cells containing a or b parameters (e.g., direct or indirect regulation by ARF7) above the validation threshold for positive regulation by ARF7 (1.5) are highlighted in red, while the cells containing values smaller than the validation threshold for negative regulation by ARF7 (1/1.5) are highlighted in blue. The mathematical model used to process the data is included as a reminder in the tab entitled "The model."

Nuclei Preps and Chromatin Shearing
Nuclei preps were made from ;5 g of root tissue using the method described by Bowler et al. (2004), with plant protease inhibitor cocktail (Sigma-Aldrich) added fresh to all buffers to a final concentration of 13. Nuclei were resuspended in 1 mL sonication buffer (10 mM potassium phosphate, pH 7, 0.1 M NaCl, 0.3% sarkosyl, 10 mM EDTA, 0.1 mM PMSF, and 13 Sigma-Aldrich plant protease inhibitor cocktail). To shear DNA to fragments of 200 to 600 bp, samples were sonicated four times for 15-s bursts with a Soniprep 150 Exponential Microprobe set at 10 mm of amplitude. The samples were kept on ice between sonications. Samples were centrifuged at 16,000 rpm for 10 min and the supernatant transferred to a fresh Eppendorf tube. At this stage, a fraction of the sonicated chromatin (50 mL) was reverse cross-linked by heating to 95°C in 0.5 M NaCl for 20 min and incubating overnight with 4 mL of 20 mg/mL proteinase K. DNA was extracted with phenol (buffered to pH 8):chloroform: isoamyl alcohol, 25:24:1, and ran on a 1% agarose gel to verify DNA shearing to an average of 400 bp in all samples.

Chromatin Immunoprecipitation
For the nonimmunoprecipitated "input" sample, 50 mL of chromatin was put aside. For antibody precipitated samples, 200 mL of sonicated chromatin was added to 1 mL immunoprecipitation buffer (50 mM HEPES, pH 7.5, 150 mM KCl, 5 mM MgCl 2 , and 0.1% Triton X-100) and incubated along with 3 mg of anti-ARF7 at 4°C on a slow-moving rotator for 4 h. Storage buffer was removed from 50 mL of Protein A Dynabeads (Invitrogen), which were immediately resuspended using the chromatin and antibody mix and further incubated at 4°C on a slow-moving rotator overnight. The magnetic beads were washed four times for 1 h with immunoprecipitation buffer and twice with water. Elution and reverse cross-linking of the immunoprecipitated and input samples were performed by heating to 95°C in 0.5 M NaCl solution for 20 min. Upon cooling, 4 mL of 20 mg/mL Proteinase K was added and the samples incubated at 55°C overnight, then at 65°C for 6 h. The magnetic beads were removed and DNA extracted from the remaining solution using phenol (buffered to pH 8):chloroform:isoamyl alcohol, 25:24:1. Precipitated DNA from the immunoprecipitated samples was resuspended in 30 mL and the input DNA in 300 mL water.

Quantitative PCR Analysis
ChIP input and immunoprecipitated samples were diluted 1 in 10 and 5 mL of this used in a 12-mL volume reaction of SYBR green master mix and 1 mM each of forward and reverse oligonucleotides (Supplemental Table  14). All qPCR reactions were performed as quadruplicate technical replicates using a Light Cycler 480 qPCR machine. ChIP-qPCR experiments are representative of at least three biological replicates.

Antibody Generation
The anti-ARF7 antibody was generated in rabbits using an Escherichia coli-expressed antigenic region from amino acids 795 to 1039 that was also used to affinity purify the resulting antiserum.

Gene Ontology Analysis
The Gene Ontology analysis was performed using BiNGO! software (Maere et al., 2005), which works as a Cytoscape plug-in (Shannon et al., 2003;Cline et al., 2007;Smoot et al., 2011). The analysis was performed using the latest Arabidopsis annotation (October 11, 2014).

Supplemental Data
Supplemental Figure 1. The whole network topology predicted by TDCor.
Supplemental Figure 2. Expression profiles in the lateral root transcriptomic data set of the ARF7 regulator and its predicted targets by TDCor.
Supplemental Figure 3. Gene Ontology analysis on genes correlated with PKL in the LR data set.
Supplemental Table 1. Assumptions for the dynamics of the regulation of gene X by gene Y in a single root cell.
Supplemental Table 2. Definition of terms for the equation describing the regulation of gene Y by gene X in a single cell.
Supplemental Table 3. List of genes used to run the algorithms (for TDCor).
Supplemental Table 4. List of interactions predicted by VBSSM.
Supplemental Table 5. List of interactions predicted by TDARACNE.
Supplemental Table 6. List of interactions predicted by TDCor. Table 7. Comparison of the performance of the algorithm based on the interactions well established in the literature. Table 8. Response components activated in different treatment of the ARF7-GR line.

Supplemental
Supplemental Table 9. Parameters of the model for ARF7-GR data analysis.
Supplemental Table 10. The position and sequences of auxin response elements.
Supplemental Table 13. Parameters used to run TDCor.
Supplemental Table 14. Oligo used in the ChIP-qPCR experiment.
Supplemental Methods 1. Analysis and proofs showing that nonlinear interactions cannot be predicted from the LR data set.
Supplemental Methods 2. Detailed description of the TDCor algorithm.
Supplemental Methods 3. Analysis of the ARF7-GR data.
Supplemental Methods 4. A step-by-step tutorial on the use of the TDCor algorithm to infer gene regulatory networks from time series transcriptomic data.
Supplemental Methods 5. Identification of potential missing links in the network.