2.1. Dynamics of nitrogen-responsive transcriptional networks show crosstalk with jasmonate signalling in Arabidopsis roots

As N is added to augment its limited bioavailability in the soil via fertilizer, we sought to characterize how a plant’s transcriptome changes upon a sudden increase in N concentration. Arabidopsis thaliana seedlings (Col-0) were grown for seven days on media with limiting nitrogen (1mM KNO₃) and then transferred to either limiting (1mM) or sufficient (10mM) KNO₃. Roots were collected immediately (0 minutes used exclusively for the 1mM KNO₃ condition as a control for mechanical transfer), 15, 45, 90 and 180 minutes after transfer (Figure 1a). To identify transcriptionally responsive genes to this change in KNO₃ concentration, we performed pairwise comparisons at each of these time points between N-sufficient and N-limiting conditions (q < 0.05). A total of 736 genes were detected as differentially expressed (DEGs), the majority of which showed distinctive and dynamic expression patterns. Specifically, 698 DEGs (94.84%) showed differential expression at only one time point (Figure 1b) and 38 DEGs (5.16%) were shared across at least two time points (Supplementary Figure 1, Supplementary Table 1). Such diversity in temporal transcriptional responses in response to a sudden increase in N suggests complex underlying regulatory networks.

Figure 1.

Transcriptional response upon nitrogen acclimation. (a) Experimental design of the time course. (b) Number of differentially expressed genes (DEGs) at each time point comparing nitrate sufficient with limiting conditions (q < 0.05). (c) Gene ontology enrichment of DEGs. (d) Heatmap of differentially expressed genes linked to N.

Figure 1. Long description

Panel A, top left, is a flowchart showing seedlings grown for 7 days on 1 m M K N O sub 3, then transferred to either 1 m M or 10 m M K N O sub 3. Roots are harvested at 0, 15, 45, 90, and 180 minutes for 1 m M, and at 15, 45, 90, and 180 minutes for 10 m M. Panel B, top right, is a Venn diagram with four overlapping ovals labeled 15 m, 45 m, 90 m, and 180 m, each with numbers indicating the count of differentially expressed genes at each time point and their intersections. For example, 353 genes are unique to 45 m, 149 to 15 m, 147 to 90 m, and 76 to 180 m, with smaller numbers in overlapping regions. Panel C, bottom left, is a dot plot with Y axis labeled G O biological processes and X axis labeled Fold enrichment. Dots represent enriched processes among differentially expressed genes, sized by number of genes and colored by F D R. Top processes include regulation of jasmonic acid mediated signaling pathway and cell wall macromolecule catabolic process. Panel D, bottom right, is a heatmap with genes listed on the Y axis (P T R 3, E R F 107, N R T 2.4, N R T 1.7, A B F 3, C C A 1, N A C 080, G S 2, A S N 2, N R T 2.2, U P M 1, T I P 2;3, N I T R 2;1, N I A 1) and time points on the X axis (15 min, 45 min, 90 min, 180 min). Color scale represents log 2 fold change, with red for upregulation and blue for downregulation.

To validate that our experimental design captured canonical aspects of the N response, we conducted enrichment analyses. When comparing our 736 DEGs against three diverse but complementary datasets that capture aspects of N-mediated transcriptional regulation, we found significant enrichment for N-responsive genes. Firstly, we found six DEGs in common with the fifty most broadly responsive N genes identified in a meta-analysis study (Canales et al., Reference Canales, Moyano, Villarroel and Gutiérrez2014) (Fisher’s exact p-value = 0.00211). In a comparable expression profiling dataset (Liu et al., Reference Liu, Niu, Konishi, Wu, Du, Sun Chung, Li, Boudsocq, McCormack, Maekawa, Ishida, Zhang, Shokat, Yanagisawa and Sheen2017), we found 27/474 genes in common (Fisher’s exact p-value = 0.00017). Finally, we found an enrichment of our DEGs in our previously published yeast one-hybrid (Y1H) transcriptional network (Gaudinier et al., Reference Gaudinier, Rodriguez-Medina, Zhang, Olson, Liseron-Monfils, Bågman, Foret, Abbitt, Tang, Li, Runcie, Kliebenstein, Shen, Frank, Ware and Brady2018) for N- and N-associated metabolism with 23/431 genes (Fisher’s exact p-value = 0.00126) (Supplementary Table 2A-C). Additionally, to further validate the robustness of our dataset, we performed a gene ontology (GO) enrichment and found that ‘response to nitrogen compound’ (GO:1901698) was significantly enriched within the DEGs (Figure 1c). When exploring the expression profiles for several nitrate transporters, nitrate assimilation genes and known nitrate-related transcription factors, we found that 11 of these genes showed a time point-specific induction, while 3 DEGs were repressed, highlighting the diverse temporal transcriptional nature of N-dependent responses (Figure 1d).

We then hypothesized that distinct transcriptional cascades result in a unique molecular response triggered at each time point upon differential nitrate availability. To identify transcriptional regulators enabling such a hierarchical response, we inferred gene regulatory networks (GRNs) with the previously identified 736 DEGs (FDR < 0.05), of which seventy-three are transcription factors (TFs) (Supplementary Table 3A). We inferred two networks using a Random Forest approach, one to explore regulatory interactions in nitrate-limiting conditions and the other in nitrate sufficiency. Changes in wiring (interactions between TFs and their downstream targets) were inferred using a machine learning approach (Spurney et al., Reference Spurney, Van den Broeck, Clark, Fisher, de Luis Balaguer and Sozzani2020) (Supplementary Table 3B,C) (see Materials and Methods ). The network in N-limiting conditions contained 277 genes and 382 interactions, while the N-sufficient network comprised 279 genes and 384 interactions (Supplementary Figure 2A,B). In 28% of these inferred interactions, we found that the cis-elements present in the promoter of the target gene were bound by their respective TF in vitro, providing validation of our approaches (O’Malley et al., Reference O’Malley, Huang, Song, Lewsey, Bartlett, Nery, Galli, Gallavotti and Ecker2016) (Supplementary Table 4). When comparing these two networks (Figure 2a), major changes in connectivity were observed between JA signalling genes.

Figure 2.

Network displaying the dynamic rewiring upon KNO₃ exposure. (a) The border colour represents the DyNet rewiring score (darker is a higher score). Yellow nodes are genes in 1 mM KNO₃ and 10 mM KNO₃ networks, blue nodes are in the 1 mM network only, red nodes are in the 10 mM network only. Activating and repressing interactions are represented by point and block arrows, respectively. (b) Top 10 TFs from network analyses with outdegree difference and Dynet rewiring score combined into a value to rank importance in network structure.

Figure 2. Long description

Panel A is a network diagram with nodes and directional edges. Three main modules are shaded: L B D 13 module at lower left, J A module at upper center, and E R F 107 module at lower right. Nodes are colored yellow if present in both 1 milli molar and 10 milli molar K N O 3 networks, blue if only in 1 milli molar, and red if only in 10 milli molar. Node borders vary in darkness, with darker borders indicating higher DyNet rewiring scores. Large diamond-shaped nodes labeled L B D 13, H S F 4, S C L 28, R A P 2.6, E R F 107, and C C A 1 are prominent. Arrows with points represent activating interactions, while block arrows indicate repressing interactions. Panel B is a table listing the top ten transcription factors: L B D 13, E R F 107, H S F C 1, C C A 1, F H A 3, T Z F 8, G A T A 27, S C L 28, H S F 4, and E R F 106. Columns are Outdegree Difference, DyNet Rewiring, and Combined, with L B D 13 scoring highest in all categories. The table is read left to right, with values: L B D 13 (37, 23.5, 2.76), E R F 107 (9, 24.5, 2.1), H S F C 1 (19, 11.5, 1.98), C C A 1 (10, 16.5, 1.76), F H A 3 (10, 11, 1.69), T Z F 8 (7, 8, 1.5), G A T A 27 (11, 10.5, 1.44), S C L 28 (8, 15.5, 1.4), H S F 4 (8, 19, 1.4), E R F 106 (13, 13.5, 1.35).

Within this set of 736 DEGs used to infer the two networks, JA signalling (‘regulation of jasmonic acid mediated signalling pathway’ [GO:2000022]) was the most significantly enriched GO category (Figure 1c, Supplementary Figure 3 [heatmap for all JA genes]); and 171 of these genes were previously identified to be transcriptionally responsive to JA (Zhang et al., Reference Zhang, Zhou, Yang and Fan2020) (Fisher’s exact p-value = 4.974e-13, Supplementary Table 2D). These genes include the JA precursor biosynthesis genes LIPOXYGENASE 3 (LOX3) (Bannenberg et al., Reference Bannenberg, Martínez, Hamberg and Castresana2009) and OXOPHYTODIENOATE-REDUCTASE 3 (OPR3) (Schaller et al., Reference Schaller, Biesgen, Müssig, Altmann and Weiler2000), the JA catabolic enzyme JASMONATE-INDUCED OXYGENASE4 (JOX4) (Caarls et al., Reference Caarls, Elberse, Awwanah, Ludwig, de Vries, Zeilmaker, Van Wees, Schuurink and Van den Ackerveken2017), as well as JA signalling repressors, including eight JASMONATE-ZIM-DOMAIN PROTEIN (JAZ) proteins (Chini et al., Reference Chini, Fonseca, Fernández, Adie, Chico, Lorenzo, García-Casado, López-Vidriero, Lozano, Ponce, Micol and Solano2007), and known JA-regulating TFs WRKY DNA-BINDING PROTEIN 18 (WRKY18) and WRKY40 (Pandey et al., Reference Pandey, Roccaro, Schön, Logemann and Somssich2010). Most of these JA signaling DEGs have reduced expression in sufficient nitrate, specifically at the 15-minute time point (Figure 3c, Supplementary Table 1). In addition to this over-representation, JA-associated genes form a highly interconnected module in both the limiting and sufficient nitrate inferred networks (Figure 2a). Despite their annotation as JA-associated genes, only 11 of the total 99 interactions across these modules are in common between the two networks (Figure 3a,b). Interactions with likely functional consequences include predicted targets of the JAZ transcriptional repressors (JAZ2, 5, 6, 7, 9, 13), as well as WRKY18 and 40, which are vastly rewired in the two N networks.

Figure 3.

Rewiring of JA-related subnetwork upon KNO₃ exposure. (a) 10 mM KNO₃ network. (b) 1 mM KNO₃ network. (c) Heatmap of differentially expressed genes linked to JA regulation and metabolism.

Figure 3. Long description

Panel A on the left shows a network diagram labeled 10 m M K N O sub 3, with circular and diamond-shaped nodes such as J A Z 10, J A Z 9, W R K Y 40, and E R F 13. Blue and red arrows connect nodes, indicating regulatory relationships. Panel B in the center is a similar network for 1 m M K N O sub 3, with more densely connected nodes including J A Z 10, J A Z 9, J A Z 7, W R K Y 40, and O P R 3. The pattern and density of arrows differ from panel A, with more red arrows. Panel C on the right is a heatmap labeled log 2 F C, with gene names listed vertically: W R K Y 40, J A Z 2, J A Z 11, J A Z 10, J A Z 9, J A Z 7, O P R 3, J A Z 6, L O X 3, W R K Y 18, J A Z 13. Time points are listed horizontally: 15 min, 45 min, 90 min, 180 min. Color intensity ranges from red (upregulation) to blue (downregulation), showing dynamic expression changes for each gene across time points.

Figure 4.

Rewiring of subnetwork linked to ERF107 and subnetwork linked to LBD13 upon KNO₃ exposure. (a) ERF107 linked 10 mM KNO₃ network. (b) ERF107 linked 1 mM KNO₃ network. (c) erf107-1 mutants have significantly shorter average lateral root lengths in 10 mM KNO₃ and 1 mM KNO₃. (d) LBD13 linked 10 mM KNO₃ network. (e) LBD13 linked 1 mM KNO₃ network. (f) LBD13 RNAi mutants have significantly more unemerged lateral roots per cm than Col-0 wild-type seedlings. ^*Indicates p-value < 0.05 using two-way ANOVA with post-hoc Tukey HSD. Box plots are centred at the data median and mark from the 25th to the 75th percentile. Individual measurements are plotted as dots. For plot in c, n = 21 (1 mM KNO₃), n = 18 (10 mM KNO₃). For plot in f, n = 14 (1 mM KNO₃), n = 13 (10 mM KNO₃)

Figure 4. Long description

Panel A at top left shows a network diagram for 10 mM K N O sub 3 with ERF107 as a central node, connected by blue and red edges to multiple genes, and CCA1 as a secondary hub. Panel B at top center shows the 1 mM K N O sub 3 network with ERF107, ERF106, and CCA1 as central nodes, with altered edge patterns compared to panel A. Panel C at top right is a box plot with x-axis labeled Col-0 and erf107-1 for both 10 mM and 1 mM K N O sub 3, y-axis labeled average lateral root length, showing significantly shorter roots in erf107-1 mutants, indicated by asterisks. Panel D at bottom left shows a network for 10 mM K N O sub 3 with LBD13 as the central node, connected to AT3G06590, LRL3, ABF3, EXT13, and others, with directional arrows. Panel E at bottom center shows the 1 mM K N O sub 3 network with LBD13 as the main hub, connected to more genes than in panel D, including HSF4. Panel F at bottom right is a box plot with x-axis labeled Col-0 and LBD13 RNAi for both 10 mM and 1 mM K N O sub 3, y-axis labeled unemerged L R per cm, showing significantly more unemerged lateral roots in LBD13 RNAi mutants, marked by asterisks. All box plots display medians, interquartile ranges, and individual data points.

2.2. Identification of key transcriptional players in the nitrogen-mediated network dynamics

Next, we delved into the identification of transcription factors (TFs) pivotal for the early transcriptional response to changes in N availability. Recognizing that TFs play a central role in the temporal rewiring of the network, we embarked on a comprehensive network analysis considering topology, regulatory rewiring and the impact of TFs on lateral root development upon perturbation (Brady et al., Reference Brady, Zhang, Megraw, Martinez, Jiang, Yi, Liu, Zeng, Taylor-Teeples, Kim, Ahnert, Ohler, Ware, Walhout and Benfey2011; Jeong et al., Reference Jeong, Mason, Barabási and Oltvai2001; Reich et al., Reference Reich, Mechler and Victor2001). To identify these key TFs, we ranked all the TFs taking into account network topology and rewiring. Specifically, we (1) quantified the outdegree of each TF, (2) computed the impact of each TF upon perturbation (Broeck et al., Reference Broeck, Spurney, Fisher, Schwartz, Clark, Nguyen, Madison, Gobble, Long and Sozzani2021), (3) used the DyNet plugin (Goenawan et al., Reference Goenawan, Bryan and Lynn2016) in Cytoscape to summarize the regulatory rewiring of each TF between the two networks (Figure 2b). We ranked the TFs by normalizing each metric and, subsequently, taking their sum. Among the most highly ranked genes, we found CCA1, which has previously been shown to link the circadian clock with N metabolism (Gutiérrez et al., Reference Gutiérrez, Stokes, Thum, Xu, Obertello, Katari, Tanurdzic, Dean, Nero and McClung2008). CCA1 also plays a role in root development: CCA1 overexpression results in altered lateral root architecture (Ruts et al., Reference Ruts, Matsubara, Wiese-Klinkenberg and Walter2012). Within our network, we found that CCA1 participated in 30 interactions in the N-limiting network and 20 interactions in the N-sufficient network, with 13 of these interactions found in both networks. The eighth-ranked gene, SCARECROW-LIKE 28 (SCL28), participated in 13 interactions in the N-limited network and 20 interactions in the N-sufficient network, with only one of these interactions found in both networks. SCL28 mutants have altered RSA with reduced primary root length (Goldy et al., Reference Goldy, Pedroza-Garcia, Breakfield, Cools, Vena, Benfey, De Veylder, Palatnik and Rodriguez2021).

Further exploration into TFs crucial for the early transcriptional response led us to ERF106 and 107. Accordingly, genes that were differentially expressed between limiting and sufficient N conditions at 45 and 90 minutes after transfer formed a module generally clustered around the top-ranked CCA1, as well as ERF106 and ERF107. Under N limiting conditions, we inferred extensive overlap in direct downstream targets for ERF106, ERF107 and CCA1 (10 shared targets), ERF107 and CCA1 (5 shared targets), and between ERF106 and CCA1 (5 shared targets) (Figure 4a,b). This redundant regulation is less prominent in the N-sufficient network where we found 2 common targets with all 3 TFs, 4 between ERF107 and CCA1, and 6 between ERF106 and CCA1. Overall, the topology of this module, specifically the redundant regulation, may indicate the importance of this module in N limiting conditions. As such, we hypothesized that ERF106, ERF107 and CCA1 played an essential role in N-mediated regulation of RSA. To support this hypothesis, we further investigated ERF107. Inferred downstream targets of ERF107 were important N-related genes such as NRT2.4 and CCA1 in the N-limiting network. Indeed, we found that erf107-1 mutants have reduced lateral root length in both limiting and sufficient nitrate (Figure 4c) (Gaudinier et al., Reference Gaudinier, Rodriguez-Medina, Zhang, Olson, Liseron-Monfils, Bågman, Foret, Abbitt, Tang, Li, Runcie, Kliebenstein, Shen, Frank, Ware and Brady2018).

2.3. Deciphering LBD13’s role in nitrogen-mediated root development dynamics

LBD13 was the top-ranked gene in terms of the number of outgoing targets and its rewiring score (Figure 2b). Within the inferred network, LBD13 was part of a network module, which consisted of 51 genes (first-degree neighbours) that were primarily responsive to changes in N at the 180-minute time point. Notably, a majority of the regulatory interactions involving LBD13 (40/43) were specific to N-limiting conditions (Figure 4e). In N sufficient conditions, regulatory interactions for some of these genes are rewired with connections to other top genes: HEAT SHOCK FACTOR C1 (12 genes, 4 genes with the same interaction type), HEAT SHOCK FACTOR 4 (9 genes, 3 genes with the same interaction type) and SCL28 (1 gene with the same interaction type) (Figure 4d,e).

While LBD13 has previously been shown to regulate lateral root development, its function in response to changes in N availability is unknown. We used an inducible RNAi line of LBD13 (Cho et al., Reference Cho, Jeon, Pandey, Ha and Kim2019) and analysed its RSA in limiting and sufficient N in the presence of dexamethasone. In limiting N, we observed a differential lateral root response: the number of unemerged lateral roots per cm primary root was increased in N-limiting conditions relative to Col-0 control (Figure 4f). However, we observed no other RSA-related mutant phenotypes related to primary and lateral roots in the LBD13 RNAi line.

Our transcriptomic time course data were obtained from RNA sequencing of bulk root samples. However, regulatory programmes comprising the initiation and elongation of lateral roots in different concentrations of available N are executed in a very small number of cells within the entire root. As the LBD13 RNAi line showed differential lateral root emergence in limiting vs sufficient N, we hypothesized that LBD13 may also play a role in the small number of cells that participate in early lateral root development. Indeed, LBD13 had a peak of expression in lateral root founder cells, which declined as the founder cells progressed through asymmetric cell division and auxin response (Cabrera et al., Reference Cabrera, Sanchez-Corrionero, Balaguer, Serrano-Ron, Barrio, Cubas, Perez-Garcia, Sozzani and Moreno-Risueno2024) (Supplementary Table 5). To further dissect the role of LBD13 in early lateral root development, we leveraged cell-type-specific single-cell expression profiling data that delineates a developmental trajectory associated with lateral root initiation in N-sufficient conditions and used dynamic Bayesian network inference to generate a TF regulatory network (TFRN) from the N-responsive TFs (Supplementary Table 3D, Supplementary Figure 5). Several biologically relevant network motifs were enriched within this early lateral root initiation network, including bi-fans, feed-forward loops (FFLs), feedback loops and bi-parallels (Milo et al., Reference Milo, Shen-Orr, Itzkovitz, Kashtan, Chklovskii and Alon2002) (Figure 5a, Supplementary Table 6). FFLs have been shown to be at the base of a plant’s adaptation to stimuli (Ma et al., Reference Ma, Trusina, El-Samad, Lim and Tang2009) and to play important roles throughout cellular development and growth, such as accelerating an output response and generating pulse-like dynamics (Joanito et al., Reference Joanito, Chu, Wu and Hsu2018). We therefore further focused on a type I coherent feedforward loop (Mangan & Alon, Reference Mangan and Alon2003) in which LBD13 participates and for which there is pre-existing evidence of the target genes regulating early lateral root development. LBD13 is predicted to regulate PDF2 and ERF107; and PDF2 to regulate ERF107. LBD13 regulates the number of emerged lateral roots per cm of primary root (Figure 4f); ERF107 regulates lateral root length (Figure 4c) and PDF2 has been shown to function in lateral root development prior to emergence (Nagata & Abe, Reference Nagata and Abe2023; Nagata et al., Reference Nagata, Ishikawa, Kawai-Yamada, Takahashi and Abe2021).

Figure 5.

Mathematical modelling of the LBD13 feed-forward loop in lateral root development. (a) The regulatory interactions in which LBD13, PDF2, ERF106 and ERF107 are involved within the N-responsive TF regulatory network during lateral root development. (b–e) Expression of LBD13, PDF2 and ERF107 modelled in (b). A wild-type condition, (c) nitrogen-limiting conditions, (d,e) and differential timing of LDB13 during lateral root development. Crosses on the graphs represent the experimental expression values at each development stage. Activating, repressing and undetermined interactions are represented by point, block and diamond arrows, respectively.

Figure 5. Long description

Panel A at the top left is a network diagram with L B D 13 at the center, connected by arrows to P D F 2, E R F 106, and E R F 107. Red arrows indicate activating, block arrows indicate repressing, and diamond arrows indicate undetermined interactions. Other nodes include R A P 2.6, J A Z 7, M Y B 45, M G P, A T 2 G 36480, H S F C 1, and M Y B 77, connected by gray lines. Panel B at the top right is a line graph with x-axis labeled time (hr) from 0 to 30 and y-axis labeled expression from 0 to 35. Three lines represent L B D 13 (red), P D F 2 (orange), and E R F 107 (blue). Key developmental stages are marked below the x-axis: Pericycle, Founder cells, Nuclear migration, Dividing Founder cells, and Auxin response. L B D 13 peaks early, P D F 2 and E R F 107 peak later. Panel C below B is labeled Nitrogen limiting conditions, with the same axes and color scheme. L B D 13 expression drops sharply after peaking, P D F 2 and E R F 107 show lower, delayed peaks. Percentage changes (+61%, +12%, +33%) are annotated. Panel D below C is labeled Differential timing L B D 13, showing a shifted peak for L B D 13 and corresponding changes in P D F 2 and E R F 107. Panel E at the bottom shows all three lines with similar axes and color scheme, summarizing expression dynamics over time. Crosses on the graphs indicate experimental data points at each developmental stage.

To analyse the expression dynamics of LBD13 and these two downstream targets in early lateral root initiation, and to link regulatory rewiring to a differential N response, we generated a quantitative model, using ordinary differential equations (ODEs); each of which included a production and degradation term that depended on its respective upstream regulators. We modelled LBD13 expression in a time-dependent manner based on cell-type expression data for single cell identities previously identified (Serrano-Ron et al., Reference Serrano-Ron, Perez-Garcia, Sanchez-Corrionero, Gude, Cabrera, Ip, Birnbaum and Moreno-Risueno2021). In addition, we included the homodimerization of PDF2 (Nagata & Abe, Reference Nagata and Abe2023), which led to a better fit of the model to the expression dynamics of its downstream target ERF107 as observed in the cell-type-specific single-cell expression data. When modelling a FFL, the two regulators, here LBD13 and PDF2, can act through an OR gate or an AND gate, where only LBD13 or PDF2 would need to be expressed to activate ERF107; or where both TFs need to be expressed to activate ERF107, respectively. The model fits the data better with an AND gate, suggesting that the expression of ERF107 is tightly regulated by multiple factors (Supplementary Figure 6). Moreover, our model suggests that the specific spatial expression of LBD13 in the pericycle cells as they transition to founder cells is key for proper lateral root initiation (Figure 5b). The model demonstrated that LBD13 in the founder cells induces PDF2 at the stage of nuclear migration and in dividing founder cells. Together, PDF2 and LBD13 activate ERF107, potentially inducing nuclear migration and founder cell divisions. To computationally predict how the dynamics among LBD13, PDF2 and ERF107 affect lateral root development upon nitrogen-limiting conditions, we adjusted LBD13 expression in our model to its expression observed in our time course data under nitrogen-limiting conditions, which is a 43% increase. After running our model under these new conditions, we observed an increased and extended expression domain of PDF2 and ERF107 (Figure 5c).

We explored the robustness of our model by changing the duration of LBD13 induction. Lower expression of LBD13 at a later time point resulted in a lower peak of ERF107 expression and no expression in the dividing founder cells (Figure 5d). In contrast, earlier expression of LBD13 at the same magnitude as in the original model, but with a shorter duration, resulted in a peak of expression in founder cells during nuclear migration (Figure 5e). In each of these cases of modulation of LBD13 expression magnitude and developmental stage, PDF2 expression duration was the same (from nuclear migration to the initiation of the auxin response), but its relative magnitude shifted during these time periods (Figure 5c–e). The results show that LBD13’s spatiotemporal expression matters for the proper transcriptional progression of PDF2 and ERF107 and thus also potentially for lateral root initiation.