Arabidopsis thaliana wild-type and mutant plants used in this study were in the Col-0 background, with the exception of the gai-1 mutant (NW63; Koornneef et al., 1985) that was in Ler background. The mutants sly1-11 (N868440), spl3-1 (FLAG_173C12; Wu and Poethig, 2006), spl9-4 (SAIL_150_B05; Wu et al., 2009) and spl15-1 (SALK_074426; Wu et al., 2009), correspond to T-DNA insertion lines while the xal2-2 mutant has a En1 transposon insertion (Garay-Arroyo et al., 2013). gai-1 and sly1-11 were provided from the Nottingham Arabidopsis Stock Centre.

Plants were grown on soil (Sunshine Mix 3 Sungro) under SD photoperiod (9 h light/15 h dark) at 22°C, and GA₃ treatment was performed as described in Pérez-Ruiz et al. (2015). The entire aerial tissue was collected for semi-quantitative RT-PCRs, while meristem enriched shoot apices were used for RT-qPCRs. All RT-qPCR assays were performed with 3–4 biological replicates (8 plants each) of 26, 36, or 45 days after sowing (das) plants, depending on the experiment.

RNA was extracted with the Direct-zol RNA MiniPrep kit (Zymo Research, United States) and cDNA was obtained with the SuperScript III RT (Invitrogen, Germany). qPCRs were performed in a StepOnePlus thermocycler (Applied Biosystems, United Kingdom) with primers enlisted in Supplementary Table 1. Relative accumulation was obtained as in Pérez-Ruiz et al. (2015) using PDF2, RNAH, and UPL7 as internal controls.

The flowering transition gene regulatory network (FT-GRN) was built based on an exhaustive search in the scientific literature. The genes were selected based on the most comprehensive reviews at the time (Srikanth and Schmid, 2011; Andrés and Coupland, 2012), and on the published dynamic models of flowering time, SAM transitions and flower development (Espinosa-Soto et al., 2004; Jaeger et al., 2013; Wang et al., 2014; Leal Valentim et al., 2015; Pérez-Ruiz et al., 2015; Dinh et al., 2017). Data supporting the FT-GRN is available in Supplementary Tables 2, 3.

Network topological properties were calculated with the igraph R package. For visualization yEd graph editor was used.

A Boolean model was used to formalize the FT-GRN interactions and study its dynamics. In short, the state of a gene i (x_i) is equal to 1 if the gene is expressed and its product is active; and x_i = 0 if it is not expressed or is inactive. The system is solved in discrete time steps. The temporal evolution of the state of gene i can be represented by:

Where the state of gene i at time t+1, is x_i(t+1), and it is described by a function F_i that depends on the state of its k regulators on a previous time step (t). Logic functions (Supplementary Table 4) were obtained from reported functional experimental data and also top-down experiments that include expression of reporter lines and RT-PCR in different mutant backgrounds, RNA-seq experiments, chromatin immunoprecipitation (ChIP)-PCR, ChIP-seq and ChIP-ChIP experiments documented in Supplementary Table 3.

For a Boolean system with n number of genes, the number of possible states is finite (Ω = 2ⁿ). The dynamics of a Boolean network can be well-described by a state transition table. In this table each row represents a state of the Boolean network a ∈ {0,1}ⁿ at time t and the corresponding state of the Boolean network f(a) at time t + 1. Therefore, the system will return to a previous visited state and then stay in a cycle of states. This cycle of the states is called an attractor. An attractor is represented by a set of states {a₀,a₁,…,a_p−1} where a_{i + 1} = f(a_i)(i0,…,p−2) and a₀ = f(a_p−1) hold. If p = 1, an attractor is called a fixed-point attractor. If p > 1, an attractor is called a cyclic attractor. In both cases, p is called the period of an attractor. The set of states that will evolve to the same attractor is called the basin of attraction (Akutsu et al., 2008).

The dynamic model was solved with both, a synchronous updating scheme, where all nodes are updated at each time step or with a random asynchronous updating scheme, where at each time step a single node is randomly selected to be updated. Simulations with both types of updating schemes were run exhaustively, meaning that all the possible initial states (2²³) were considered for the simulation. Input nodes were considered as constant. Both schemes gave the same 32 attractors in the FT-GRN model. The dynamic analysis of the model was performed with R (R Core Team, 2016) packages BoolNet (Mussel et al., 2010) and BoolNetPerturb (Martínez-Sánchez, 2020). Code for all simulations in this work is available at the repository https://github.com/CaroChavez/FT-GRN (Chavez, 2022).

The Boolean model was converted to a continuous system using ordinary differential equations (ODE) following the method developed by Mendoza and Xenarios (2006) and Sánchez-Corrales et al. (2010) to automatically translate the topology of a discrete regulatory network into a continuous standardized qualitative dynamical system and study its behavior. The Boolean functions were approximated to continuous differentiable sigmoidal functions that are bound between 0 and 1 (Sánchez-Corrales et al., 2010). For each node of the FT-GRN, we had a function w_i that is the continuous form of the truth table describing the state of the node x_i. Thus, the logical functions were replaced with a set of continuous functions that satisfy Zadeh’s rules of fuzzy propositional calculus (Zadeh, 1965). The continuous system was modeled as a series of coupled ODE, describing the activation of each node with an equation of the form:

The right-hand side of the equation has two parts: a sigmoid function to account for the rate of production and a linear term to account for the decay of x_i. The parameter h determines the strength of the interactions and controls the curve form: a step function, a logistic function, or a straight line. This sigmoid function was constructed to pass through the points (0, 0), (0.5, 0.5), and (1, 1) for any positive value of h (Sánchez-Corrales et al., 2010). The parameter γ_i is the decay rate of node x_i, and w_i is the continuous form of the Boolean function for node x_i, transformed using fuzzy logic following the rules:

The explicit form of the FT-GRN ODE system including w_i can be found in Chavez (2022). The parameters analysis of the FT-GRN continuous model showed that the model gives similar attractors to the Boolean version for ranges of 10–50 and 0.5–1 for h and γ, respectively. Parameters values selected for the simulation were: h = 20, γ_i = 1. We set the parameter γ_i for all nodes equal to 1, so that the expression level of the node at its stationary state is merely determined by the degree of truth of the fuzzy proposition w_i. By giving h > 0, and γ_i ≥ 1, the equations keep the variables x_i in the closed interval [0, 1]. This means that the methodology is designed to describe normalized values of activation, not absolute values of concentration (Mendoza and Xenarios, 2006; Sánchez-Corrales et al., 2010). The FT-GRN continuous model gives gene steady states close to either 0 or 1; the steady states of the system are then comparable to the FT-GRN Boolean model phenotypes. The numerical analysis of the ODE system was conducted using R with packages deSolve (Soetaert and Petzoldr, 2010) and BoolNetPerturb (Martínez-Sánchez, 2020).

Previous work showed that XAL2 participates in floral transition, specifically, it was established that CO positively regulates XAL2 while SOC1 represses it. Also, XAL2 upregulates SOC1, LFY, AP1, and TFL1 (Pérez-Ruiz et al., 2015; Figure 1A). All these regulatory interactions were documented under a LD photoperiod, and they were observed when the entire seedling was collected. But the strongest XAL2 mutants phenotypes were observed under SD photoperiod (Pérez-Ruiz et al., 2015). Since flowering transition relies on developmental and physiological signals under SD conditions (Wang et al., 2009; Quiroz et al., 2021), we searched for the role of XAL2 in flowering regulation in response to aging or GA₃ addition.

XAL2 mRNA levels increase in aerial tissues during plant development, peaking around 42 days after sowing (das), and GA₃ treatment accelerated XAL2 accumulation when they growth in short-day (Figure 1B). Besides, XAL2 is induced by GA₃ in shoot apices as soon as 26 das, compared to non-treated plants, and it is downregulated in a DELLA gain of function mutant; the gibberellic acid insensitive 1 (gai-1; Koornneef et al., 1985), and the sleepy1-11 (sly1-11) mutant affected in GA signal transduction (Eriksson et al., 2006), indicating that XAL2 is positively regulated by GA (Figure 1C). To investigate the role of XAL2 in the regulation of genes involved in flowering transition in response to GA₃, we treated wild-type and xal2 plants with this hormone and collected shoot apices at 26 and 36 das to evaluate early and late developmental events. Other than XAL2, only SOC1 was induced earlier (at 26 das) by GA₃ (Figure 1D and Supplementary Figure 1). Interestingly, from eight genes analyzed at both times (Figure 1D and Supplementary Figure 1), XAL2 mediates the upregulation of SOC1, LFY, and AP1 in response to GA₃ (Figure 1D). While GA₃ upregulates SPL9, SPL15, FD, and FRUITFULL (FUL) genes independently of XAL2 (Supplementary Figure 1). We also found that XAL2 transcript accumulation was reduced in the spl9-4 and spl15-1 mutants (Wu et al., 2009) under SD conditions, but not in the spl3-1 background (Figure 1E), suggesting that SPL9 and SPL15 mediate XAL2 expression in response to aging.

Thus, XAL2 is induced early in development probably by SPL9 and SPL15, and in response to GA signaling; but most importantly, XAL2 partly mediates the response of SOC1, LFY, and AP1 to GA₃. These results further document the various roles of XAL2 during plant development. Though, the question remains: if XAL2 induces those genes similarly in SD and LD, why the xal2 mutant has a more severe flowering transition phenotype under SD? Therefore, we were interested in analyzing its role in the context of the GRN.

For a GRN system-level understanding of XAL2 function in the flowering transition, we integrated the novel results (Figure 1) and previous experimental data into a flowering transition-GRN (FT-GRN) model (Figure 1A and Supplementary Table 2). Based on XAL2 regulatory interactions we incorporated endogenous aging and GA signals, essential for the flowering response under SD conditions, and the vernalization response which previous models did not consider. Therefore, we aim to resolve a new FT-GRN that could respond to different environmental and internal signals.

The FT-GRN presented here (Figure 2) was derived from a hand-curated literature search reported on the molecular regulators of flowering transition in A. thaliana. It integrates regulatory interactions from 146 publications and 45 genes (Supplementary Tables 2, 3), including the novel XAL2 data shown here (Figure 1). The network nodes represent genes and gene products such as microRNAs or proteins (mainly transcription factors). It also incorporates GA hormone signaling (GA), and processes such as aging (AGE) and vernalization (VER). Genes were selected according to their relevance in the plant’s flowering response to external and internal cues, and functionally redundant genes were collapsed into a single node of the network (Supplementary Table 2).

The FT-GRN has five input nodes: CO, VER, FCA, AGE, and PNY. CO was selected as the input for the simulation of a LD photoperiod because flowering response relies on this protein, although there are many circadian cycle proteins upstream of CO (Suárez-López et al., 2001; Valverde et al., 2004). Vernalization (VER) corresponds to a prolonged period of cold that resembles winter and it is important for FLC epigenetic silencing (Sheldon et al., 2002; Whittaker and Dean, 2017), but also in the upregulation of AGAMOUS-LIKE 24 (AGL24) independently of FLC inhibition by mechanisms still unknown (Michaels et al., 2003). FCA is a RNA binding protein that acts within a larger protein complex in FLC epigenetic silencing in response to aging, independently of photoperiod and vernalization, thus it has been considered as an autonomous pathway of flowering transition (Macknight et al., 1997; Wu et al., 2020). FCA is also involved in the post-transcriptional regulation of miR172 and SVP depending on ambient temperature (Lee et al., 2007; Jung et al., 2012). Details about FCA transcriptional regulation are unknown, but it regulates itself at the post-transcriptional level to achieve correct levels of its active form (Quesada et al., 2003). The AGE node represents a growth period after some molecules are depleted, accumulated, or become active/inactive. For example, GA levels (Eriksson et al., 2006), FD (Abe et al., 2005; Klepikova et al., 2015) and XAL2 mRNAs (Figure 1) increase as the plant grows. On the contrary, miR156 levels decrease during development (Wang et al., 2009; Jung et al., 2011), like FLC (Searle et al., 2006; Klepikova et al., 2015). PNF and PNY (both included in the PNY node) were included because they are essential in flower meristem formation (Smith et al., 2004; Kanrar et al., 2008), although they have not been considered classical flowering genes. PNY regulation is unknown, and thus we kept it as an input in the FT-GRN. Therefore, the model inputs simulate the response to multiple signals from both the environment and endogenous cues.

The FT-GRN also incorporates the FT and TSF genes (both in the FT node), which are induced by LD photoperiod in the leaves, and then, the protein is transported by the phloem to the SAM (Wigge et al., 2005; Corbesier et al., 2007). FT/TSF and TFL1 compete for the binding of FD and FDP, resulting in the promotion or repression of flowering transition, respectively (Abe et al., 2005; Zhu et al., 2020). The GA node represents high levels of this hormone and the concomitant DELLA proteins degradation (Davière and Achard, 2013). We kept separated SPL9/15 (in the SPL9 node) and SPL3/4/5 (in the SPL3 node) because they form different clades with unique functions in plant development (Xu et al., 2016). Similarly, AP2 and AP2-like genes were separated into two nodes because AP2 is related to floral development and AP2-like to flowering repression (Mathieu et al., 2009; Yant et al., 2010), although they also shared common functions. The FT-GRN includes MADS-box transcription factors such as SVP and FLC that function as flowering repressors, and others that promote flowering: SOC1, FUL, AGL24, and XAL2. Finally, this FT-GRN incorporates the flower meristem identity genes: LFY and AP1 (together with CAL regulation); and the inflorescence meristem identity gene TFL1 (Ratcliffe et al., 1999).

The first result from this network is that most of the nodes are highly connected among themselves (Figure 2). The FT-GRN has 23 nodes and 121 edges (Figure 2 and Supplementary Table 3) and a directed net density of 0.24 that measures the proportion of edges present within all possibilities in the directed network. In fact, the whole network shows an average of 10.5 regulatory interactions per node (average degree). The analysis of the FT-GRN topology shows that the FT-GRN has a high global transitivity of 0.52, which measures the probability that the adjacent nodes of a node are connected. The network reciprocity is 0.33; this number represents the proportion of reciprocated ties in a directed network; so, there is a relatively high proportion of direct feedback loops in the FT-GRN. Therefore, the FT-GRN is highly connected (density and average degree), topologically redundant (transitivity), and has a high number of feedback loops (reciprocity).

As expected, despite the high connectivity, some nodes are more connected than others (Supplementary Figure 2). The nodes with the highest degree (total number of edges) are SOC1 and AP1. Also, FUL, LFY, AP1, and SOC1 showed the highest indegree values, which measures the number of regulators of each gene (Supplementary Figure 2). It is interesting that FUL has an indegree closer to SOC1, which is considered an integrator of the flowering pathways, and surprisingly, the number of FT regulators is like the number of regulators for other nodes in the FT-GRN. Additionally, the number of genes regulated by a particular node in the network, indicated by the outdegree, showed that SOC1, SVP, GA, and AP1 have the highest values (Supplementary Figure 2), accordingly with their importance as regulators in the flowering transition network.

Contrary to what would be expected if independent or parallel genetic pathways converged to a few integrators to induce flowering transition, we found that all nodes are highly interconnected among themselves. Therefore, the so-called integrators were not identified in the network based on purely topological information. Even if the topology of a network constrains its dynamics, the topological properties of the network such as the node’s degree are not enough to understand the function of specific components whiting the whole system, the specific dynamic details captured in the node’s logic rules which are different for example in FT and SPL9, can explain why even if they have the same degree, they have different roles in the dynamics of the system. One way to explore the node’s function in the context of the system-level dynamic properties of the model is by perturbing individual components (see section “Differentiation and Plasticity Are Emergent Properties of the FT-GRN Signal Sensing and Integration”).

We hypothesized that a multi-stable FT-GRN module underlies the meristem identity transitions in the development from the vegetative to the reproductive state (Figure 3). In this study, we have integrated many of the molecular components that have been characterized for flowering transition in response to different environmental and internal signals, including genes like the SPL transcription factor family that have been shown to participate in developmental transitions and their negative regulation by miR156 (Huijser and Schmid, 2011). During aging, miR156 is downregulated with the concomitant increase of SPLs and miR172 expression (Wu and Poethig, 2006; Wang et al., 2009; Wu et al., 2009). The latter represses AP2 and AP2-like TF’s family that act as repressors of flowering transition (Mathieu et al., 2009; Yant et al., 2010). The miR156-SPLs-miR172 module is also important for vegetative transition from juvenile to adult plants (Wu et al., 2009).

We defined the upregulation of SOC1 (Samach et al., 2000; Immink et al., 2012) and TFL1 expression (Liu et al., 2013; Baumann et al., 2015) as markers of VM to the IM transition, while the FM identity acquisition was defined by the expression of LFY and AP1 with the concomitant repression of TFL1 (Shannon and Meeks-Wagner, 1993; Liljegren et al., 1999; Ratcliffe et al., 1999).

The FT-GRN dynamic behavior was formalized with a Boolean model (Supplementary Table 4) that captures the regulatory logic between the genes necessary and sufficient for SAM phase transitions and flowering. First, we simulated the dynamics of a wild-type FT-GRN by an exhaustive dynamic exploration, initializing with all 2²³ possible states. Despite the immense size of the state space of the FT-GRN model, both, synchronous and asynchronous updating schemes gave the same 32 fixed-point attractors, and zero cyclic attractors (see section “Materials and Methods” and Supplementary Figure 3). This result indicates that the non-linear interactions integrated into the model, based on available experimental data, effectively canalize a large and complex dynamic system into a relatively small set of attractors, demonstrating the system’s robustness to initial conditions.

We recovered attractors with similar expression patterns as those observed during SAM phases of development (Supplementary Table 5). Thus, the attractors were classified into the different meristem types: JVM (10 attractors), AVM (2 attractors), IM (4 attractors), and FM (16 attractors). There is more than one attractor per meristem type depending on the initialized input states, thus, we use the term meristem phenotype to distinguish them from the individual attractors obtained. For instance, the model recovers two attractors with the AVM phenotype that differ in the state (on or off) of the CO node, which simulates a LD or SD photoperiod, respectively. The consensus states for the attractor’s phenotypes are shown in Figure 4.

Based on the analyses of the attractors recovered by the model, the conditions driving the JVM to the AVM state are the plant aging, or the acceleration of vegetative development when plants are vernalized. The model predicts that vernalization has an impact on the juvenile-adult phase change because it downregulates FLC, and FLC induces AP2 and AP2-like genes and represses SPL9 and MIR172 (Deng et al., 2011; Mateos et al., 2015). As expected, downregulation of the flowering repressors, FLC and SVP, is necessary for the AVM to become an IM. Likewise, the formation of the FM in the model requires LFY and AP1 positive feedback and the presence of PNY to block the expression of TFL1, driving the system to the floral state. Therefore, the FT-GRN proposed here constitutes a multi-stable, dynamical system’s mechanism that sustain previous observations underlying SAM phase transitions.

Some expression patterns obtained from the model attractors differed from the observed profiles in planta (Figure 4 and Supplementary Table 5). Such discrepancies can be interpreted in terms of the discrete nature of the Boolean model that is limited to qualitative but not quantitative changes in expression or activity. For example, FD is expressed in the JVM, and its expression rises in the AVM and IM (Abe et al., 2005; Wigge et al., 2005); the model does not account for the low expression levels of FD at early stages, but it recovers that FD is upregulated at the SAM in later stages of development. Besides, FD mRNA has been detected in the FM (Abe et al., 2005); but it is negatively regulated by AP1 (Kaufmann et al., 2010b; Pajoro et al., 2014b). Therefore, the model predicts that this gene is turned off in the FM attractors.

Likewise, there is evidence suggesting that SOC1 and AP1 downregulate XAL2 (Pérez-Ruiz et al., 2015). In the model, these regulations result in the complete repression of XAL2. Also, RNA in situ hybridization of miR156 shows it at the shoot apex during the plant life-cycle, although its levels decline as the plant ages (Wang et al., 2008, 2009). In the Boolean model miR156 is present in the JVM and absent in the AVM and IM. Basal GA levels are present at the SAM throughout the life span (Eriksson et al., 2006), but they increase in the mature stages of development. In the Boolean model, GA is absent in the vegetative attractors and present in the IM (Figure 4). This can be interpreted as increased GA levels during the flowering transition, rather than their absence in the JVM.

Overall, the model recovered the great majority of documented gene expression patterns during A. thaliana developmental stages.

Zhang et al. (2021) recently published a study in which they used single-cell RNA sequencing (scRNA-seq) to distinguish seven cell types and subgroups from 23 cell clusters based on genes specific expressions that distinguish each cell type in the shoot and leaves (Zhang et al., 2021). This information can be visualized in two-dimensional maps using the uniform manifold approximation and projection (UMAP), and we downloaded those images from their platform http://wanglab.sippe.ac.cn/shootatlas/. Sixteen genes of our network are expressed in different cell types, while FT and AP1 are absent, confirming that the shoots and leaves are at the vegetative stage. Additionally, there was no expression of miR156 and miR172 in the single-cell data, probably due to their instability.

As seen in Supplementary Figure 4, even when some genes are expressed in the same cell type, they have apparent differences in the clusters in which they are identified. For example, SOC1 is densely expressed in cells that give rise to mesophyll cells (MC), and a group of undefined cells (UC), while AGL24 is expressed only in one cluster of the MC. Conversely, AGL24 has a broader distribution in cells that will give rise to vascular cells (VC), proliferating cells (PC) and is highly expressed in shoot meristematic cells (SMC) and epidermal cells (EC) in contrast to SOC1. Similarly, SVP is expressed in all cell types, whereas FLC is restricted to five of them (Supplementary Table 6).

SPL3 and SPL9 are expressed differentially. While SPL3 mRNA accumulates in EC, VC, and companion cells (CC), SPL9 has higher expression in the MC and different clusters of the EC-type than SPL3 (Supplementary Figure 4), supporting the idea that they have at least partial independent roles during the vegetative phase of development.

Another example that caught our attention was the expression of FD. It is known that TFL1 binds FD in the absence of FT (Zhu et al., 2020). Both genes expression overlaps in the SMC, similarly to the inflorescence meristem. However, FD is also expressed in EC, PC, CC, and a different cluster of VC (Supplementary Figure 4). This opens the possibility that FD may act independently of TFL1 or binds other proteins. The fact that TFL1 was detected in the VC made us suspect that it could be mistaken with other members of the FT family; thus, we analyzed the expression in a root single-cell assay http://wanglab.sippe.ac.cn/rootatlas/ (Zhang et al., 2019). TFL1 is also expressed in root VC (Supplementary Table 6), and FT was not found, suggesting that TFL1 mRNA might be present in the vasculature of vegetative organs. Moreover, scRNA-seq from the root, showed that AGL24 is almost undetectable with scarce presence in lateral roots and root hair cells, while SOC1 is highly expressed in different cell types. Further experimental functional data are now required of each SAM stage to be able to recover additional interactions and cell-specific attractors.

Robust GRNs should buffer noisy signals and stochastic biochemical processes while preserving cell phenotypes (Kitano, 2004). The FT-GRN model’s functional robustness was assessed by the invariance of attractors in response to perturbations in the logical rules. First, we tested the robustness of attractors to perturbations in the Boolean functions; afterward, we compared the state transition robustness of the model to random models with the same topological properties; and finally, we approximated the Boolean model to a continuous model to test if the attractors were independent of the modeling framework and the kinetic parameters of the functions.

For each node, copies of the model were constructed with one-bit random systematic perturbations to the truth tables of the Boolean functions. The attractors recovered by each perturbed model were compared with the attractors of the unperturbed model (Figure 5A). The perturbed networks either recovered identical attractors to the original model, a subset of the original attractors, or both original and novel attractors. Only when the input nodes were perturbed, fewer attractors were recovered compared to the original FT-GRN (Figure 5A). Furthermore, the mean basin size of the new attractors in the perturbed models were relatively small compared to the original ones (Figure 5B).

As expected, the perturbations that recovered identical attractors, were the ones that implied nodes with the highest indegree (Figure 5C) and higher redundancy in their regulation. Finally, we tested if some attractors were susceptible to be lost by random perturbations to the Boolean functions and we found that the 32 attractors were maintained with the same frequency when the network was randomly perturbed (Supplementary Figure 5). Thus, the functional robustness analysis revealed that the FT-GRN attractors are robust to perturbations in the Boolean functions (Figure 5).

A biological network is assumed to be robust to small perturbations (Kitano, 2004). We expected attractors in real networks to be less susceptible to noise than attractors in randomly generated networks, since biological processes are relatively stable and emerge even in the face of ever-changing conditions prone to stochasticity. The functional robustness was compared between the FT-GRN model and randomly generated models with the same topological properties. As expected, the FT-GRN model was considerably more robust than random models (Supplementary Figure 6). Therefore, robustness of the FT-GRN is not related only to its topology, but to the dynamics of its gene regulatory interactions.

The Boolean model was then transformed to a continuous system of ODE to test whether the mathematical formalism used was determinant of the attractors obtained. This was not the case, as the continuous FT-GRN model recovered the same steady states as the Boolean model (Supplementary Figure 3). The parameters analysis of the FT-GRN continuous model showed the attractors were stable to a range of values in the kinetic parameters. This result indicates that the attractors recovered in the original FT-GRN Boolean model are not an artifact due to the discrete nature of the Boolean model. Instead, such attractors emerge from the overall network topology and logical structure, independently of the kinetic parameters. Besides, it supports the thesis that Boolean models can recover the complexity of biological regulatory mechanisms and thus, they are powerful tools to unveil how genotypes map into phenotypes.

To summarize the results obtained from the robustness analyses, the attractors recovered by the FT-GRN model are stable to discrete perturbations of the Boolean functions and represent the expression patterns of the SAM and FM in the development of A. thaliana. Additionally, the robustness of the FT-GRN is not a causality related to the topological properties of the network, but a consequence of the biological regulatory mechanisms involved. Furthermore, the specific meristem configurations are emergent properties of the logic behind the regulatory interactions of the FT-GRN, and they are independent of kinetic parameters.

For a further validation of the FT-GRN model, we tested its ability to recover the observed mutant plant phenotypes compared to those observed in wild-type (WT) plants. Single loss-of-function (LOF) mutants and constitutive expression (CE) lines were simulated by maintaining each node state at 0 or 1, respectively. We assessed if novel phenotypes were recovered compared to the WT plants simulation (Figure 6). For this purpose, the cumulative basin of attraction size for a particular phenotype, corresponding to the sum of the state space that converge into attractors with the same meristem type in relation to all the state space, was compared between the mutant and the WT. Finally, the in silico mutant predictions were compared to the experimental data available for each mutant phenotype to assess if such data was recovered by the model (Figure 6 and Supplementary Table 7).

Most of the LOF simulations showed similar attractor’s phenotypes as the WT (Figure 6A), consistent with most of the single LOF mutant plants having mild phenotypes (Supplementary Table 7). However, a few LOF mutants completely lost some of the WT attractors, also changing the basin size of the remaining attractors in the state space. This was the case for low levels of GA, soc1 and fca mutant simulations, in which the IM was lost, and the basin size of AVM attractors was increased compared to the WT. The same happened to the JVM basin of the fca mutant. In agreement with this, the ga1-3 mutant is late flowering in a LD and in a SD, flowering fails completely (Wilson et al., 1992). soc1 and fca plants are very late flowering (Koornneef et al., 1991; Onouchi et al., 2000), and the latter produce more vegetative juvenile and adult leaves than WT plants (Telfer et al., 1997).

Of interest, the spl9/spl15 simulation showed that plants remain at the JVM state and transit to the IM without developing an AVM, in agreement, the spl9 spl15 double mutant has a prolonged juvenile phase (Schwarz et al., 2008; Xu et al., 2016). The PNY LOF simulation does not show the FM attractors. Correspondingly, the pny pnf double mutant does not develop flowers when the flowering time genes are upregulated (Smith et al., 2004). Conversely, nullifying TFL1 in the model predicts only FM attractors, losing the vegetative and IM states, likewise the tfl1 plants flower very early, and their SAM is rapidly transformed into a terminal flower (Shannon and Meeks-Wagner, 1991).

By observing some puzzling results when we compared the mutant simulations with the documented mutant phenotypes, the model revealed a connection between the vegetative phase change and its effect in the flowering transition. In particular, the simulation of FLC or miR156 LOF mutants recovered a smaller JVM basin compared to the one observed in the WT, but they also showed an increased AVM basin compared to the WT. These basin changes may be interpreted biologically as a faster vegetative phase change. It has been reported that 35S::MIM156 lines that show reduced function of miR156, only produce rosette leaves with adult morphology (Franco-Zorrilla et al., 2007; Xu et al., 2016), while flc mutation promotes adult vegetative traits both dependently and independently of the flowering time (Willmann and Poethig, 2011). Nonetheless, LOF of FLC or plants with reduced miR156 activity, also have an early flowering phenotype, that was not recovered by the simulated LOF mutants due to the presence of the flowering repressor SVP, which is independent of miR156 and FLC levels, preventing early flowering transition in the model. These results suggest unknown negative regulators of SVP, which would be necessary for the early flowering phenotype in the model when FLC or MIR156 are mutated. The double LOF mutants svp flc or miR156 svp simulations had a larger IM basin size compared to the WT (Supplementary Figure 7). svp flc and svp miR156 had an IM basin size of 25 and 18.8% of the state space, respectively, while the IM basin size was 12.5% in the WT, this predicts early flowering phenotypes and shows that SVP has a redundant role as a flowering repressor.

Furthermore, the svp simulation resulted in reduced AVM and JVM basins, a larger IM basin, and a cyclic attractor with an AVM/IM phenotype (Figure 6A). Accordingly, the svp mutant plants are early flowering, and they produce fewer rosette and cauline leaves than wild-type plants (Hartmann et al., 2000). The predicted AVM/IM phenotype in svp, suggest that this gene is relevant for maintaining the vegetative stage. The role of SVP during the vegetative phase beyond its role as a flowering repressor needs further research.

lfy mutant plants have partial transformations of flowers into inflorescence shoots, and secondary flowers generally arise from the axils of the outer floral organs; these plants also have additional secondary inflorescence shoots subtended by cauline leaves (Weigel et al., 1992). In agreement, it is remarkable that the lfy simulation recovered novel cyclic attractors with a mixed FM/IM phenotype at the expense of a slight reduction in the basin of FM attractors under SD conditions. Besides, the AP1 LOF simulation does not present the FM attractors while it recovered chimeric inflorescence-floral meristem (IFM) attractors as well as additional JVM/FM/IFM cyclic attractors. The novel attractors of the AP1 LOF simulation are consistent with the observed phenotypes in real ap1-1 mutant flowers, including the formation of bract-like first whorled organs and floral buds in their axils (Irish and Sussex, 1990).

In contrast to LOF mild phenotypes, simulations of constitutive expression (CE) lines predicted more drastic alterations, including loss of some types of meristems and novel meristem phenotypes (Figure 6B). CE of the flowering repressors FLC, SVP, AP2, and AP2-like genes, resulted in the disappearance of the IM attractors and a larger basin of the vegetative meristem ones. The predictions agree with phenotypes observed in the overexpressor lines of these genes, which are late flowering and, in some cases, remain vegetative throughout their complete life cycle (Supplementary Table 7; Michaels and Amasino, 1999; Chen, 2004; Masiero, 2004; Lee et al., 2007; Mathieu et al., 2009; Yant et al., 2010). In the FCA CE simulation, the JVM basin was reduced compared to the WT, and the AVM state disappeared, while the IM basin of attraction increased. The FCA overexpressor line is early flowering, although it is not as drastic as the simulation, probably because there are different ratios in the abundance of FCA alternative transcripts, adding another level of regulation that is not considered in the model (Macknight et al., 1997). Like the FCA CE simulation, the AVM attractors in the miR156 CE simulation disappeared, but in contrast to CE of FCA, a larger JVM basin was obtained in the miR156 CE line than in the WT. This result agrees with the prolonged juvenile vegetative phase phenotype observed in miR156 overexpressor plants (Wu and Poethig, 2006; Schwarz et al., 2008).

Conversely, the JVM attractors disappeared in the SPL3 and SPL9 CE simulations, and instead, a larger AVM basin was predicted compared to the WT. Indeed, SPL3, SPL4, SPL5, SPL9, and SPL15 overexpressor plants transit from their juvenile phase to further stages faster, and adult traits appear earlier than in WT plants (Cardon et al., 1997; Wang et al., 2008; Yamaguchi et al., 2009). In SOC1 CE simulation, the vegetative state of the SAM is lost, which correlates with the observed effect of SOC1 overexpression, causing very early flowering in both SD and LD (Borner et al., 2000). In addition, the model predicts that the PNY CE forms only FM attractors; however, this overexpression line has not yet been described in the literature.

In some CE simulations, new phenotypes were predicted besides the ones observed in WT plants. Increased GA response simulation had larger AVM and IM basins and a smaller JVM basin compared to the WT, but it also recovered additional cyclic JVM/AVM attractors. Mutants in SPINDLY (SPY), a GA signaling repressor, resemble wild-type plants that have been repeatedly sprayed with GA₃; spy plants and GA₃-treated WT plants are early flowering and they also produce adult vegetative leaves earlier (Jacobsen and Olszewski, 1993; Telfer et al., 1997). Similarly, simulation of FT CE also predicted a smaller JVM basin and recovered additional cyclic attractors with a JVM/AVM phenotype. Indeed 35S::FT plants are early flowering (Kardailsky et al., 1999; Kobayashi et al., 1999; Hanzawa et al., 2005), but there is no information related to the morphological traits like the presence of abaxial trichomes distinctive of mature leaves versus juvenile leaves in these transgenic lines. Finally, other CE mutant simulations had novel phenotypes at the expense of losing others. For example, AP1 and LFY CE simulations, predict that the SAM converts into the FM or a chimeric IFM meristem. Actually, 35S::AP1 plants are early flowering, and the SAM becomes a floral meristem that produces a terminal flower (Mandel and Yanofsky, 1995), besides, 35S::LFY plants are also early flowering and they develop solitary flowers instead of secondary shoots (Weigel and Nilsson, 1995). On the contrary, the FM disappeared in TFL1 and XAL2 CE simulation, and a chimeric IFM meristem is predicted, as has been reported in the literature for these mutant lines (Ratcliffe et al., 1998; Pérez-Ruiz et al., 2015).

Although most of the LOF and CE model predictions were congruent with mutant plants described in the literature (Figure 6), some discrepancies were found. The model did not recover the chimeric IFM observed in the overexpression lines of AGL24/SVP and the agl20-101D (Yu et al., 2004; Liu et al., 2007), which are similar to the 35S:XAL2 phenotype (Pérez-Ruiz et al., 2015). Our previous discrete model, had different conditional rules for loss of function mutants, the WT expression, and gene overexpression lines (Pérez-Ruiz et al., 2015). In order to translate the discrete multi-valued dynamics into a Boolean model, we followed the methodology in Chaos et al. (2006). Therefore, the conditional induction of TFL1 in those overexpressors was lost.

We also found that the flowering time phenotypes of co, ft, fd, agl24, ful, and xal2 mutants (Supplementary Table 7) were not reflected in changes of the size of the VM or IM basins (Figure 6). This can be explained because the model is not quantitative and as such, it is not meant to give quantitative information of flowering time as other models have done (Jaeger et al., 2013; Wang et al., 2014; Leal Valentim et al., 2015; van Dijk and Molenaar, 2017; Haspolat et al., 2019).

A third group of in silico mutants were partially recovered because they lost the IM state which correlates with their late-flowering phenotypes, but they still recovered the FM states, which should be lost. This group includes the simulation of GA insensitive mutants, SOC1 and FCA loss-of-function, and constitutive expression of AP2, AP2L, SVP, and FLC. The few discrepancies between the model’s predictions and experimental mutant data yielded novel hypotheses that should be tested experimentally.

We further investigated the role of XAL2 in the context of the FT-GRN. We explored XAL2 mutant phenotypes simulating SD conditions by keeping CO input turned off in the model and we also analyzed the GA response in this condition (Figure 7). Simulation of wild-type plants treated with GA shrank the JVM basin size and enlarged the AVM and IM state basins compared to non-treated plants in SD, this agrees with the promotion of flowering by GA₃ treatment in SD. However, simulation of xal2 mutant treated with GA₃ showed a smaller basin size of the IM state than the WT treated with GA.

In SD, A. thaliana flowering depends on endogenous cues like hormonal and aging signals (Quiroz et al., 2021). The miR156 – SPLs – miR172 – AP2-like module is important for aging developmental changes, and SPL9 and SPL15 positively regulate XAL2 (Figure 1E). Simulation of spl9/spl15 in SD lacked the AVM stage, congruent with spl9 spl15 plants having a longer juvenile phase (Schwarz et al., 2008; Xu et al., 2016). However, simulating spl9/spl15 treated with GA in SD recovers the AVM and behaves similarly to GA-treated WT plants, suggesting that GA signaling is at least partially independent of SPL9/SPL15. Interestingly, when the double mutant spl9 xal2 was simulated in SD, the IM state was lost. Moreover, when spl9 xal2 double mutant was treated with exogenous GA in silico, it could not recover the IM basin size and instead had a larger vegetative basin size, showing that the introduction of xal2 null-function mutant compromised the GA response. The model predictions support the notion that XAL2 is part of the flowering response to GA in SD, but generation of spl9 xal2 and spl15 xal2 double mutants and empirical confirmation are required.

By considering the attractors recovered by the FT-GRN as specific meristem types, we were interested in understanding how the SAM stages are acquired, and how the sequential pattern of such transitions emerges. Furthermore, we wanted to understand how such phase transition pattern is altered under different environmental and developmental cues. First, we searched for the nodes that could reshape the attractors landscape in a continuous model; then, we explored how to induce differentiation by transitory perturbations of the nodes in the Boolean model and finally, we found that the temporal sequence of differentiation observed in A. thaliana is a response to input signals propagation throughout the FT-GRN.

The attractors landscape (Waddington, 1957) corresponds to the space of dimension “n” (n = the number of nodes in the network), formed by the attractor states at the bottom of the “landscape” valleys, together with all the other states of the system. Such n-dimensional space emerges from the network structure, topology and dynamics, and it constrains the transitions among attractors and also underlies generic time-ordered patterns of temporal transitions among cell types that are observed in actual organismal morphogenesis (Villarreal et al., 2012; Davila-Velderrain et al., 2015b). By recovering and analyzing such attractors landscape for the network proposed in this study, we uncovered which nodes of the FT-GRN may be key for phase transitions. Specifically, the attractor landscape of the FT-GRN continuous model was explored by modifying the decay rates of each node (Davila-Velderrain et al., 2015b) and identifying those that were key for the transition from one SAM phase to another one using previously proposed approaches (Alvarez-Buylla et al., 2008; Huang et al., 2009; Villarreal et al., 2012; Davila-Velderrain et al., 2015b).

In silico, it was possible to transit between all meristem types when forcefully maintaining specific nodes shutdown (Figure 8). We found that the vegetative phase transition is controlled by: MIR156, FLC, SVP, SPL9, VER, and AGE nodes in the model, whereas FCA, SOC1, GA, and SVP controlled flowering transition. Shutting-down TFL1 results in FM fate while silencing AP1 or PNY mediated FM transition to other meristem types. For example, we found the appearance of an additional chimeric inflorescence-floral meristem (IFM) phenotype can be reached by the system when AP1 is silenced if it starts from the FM state. Not all nodes have the same propensity to modulate differentiation; altering one of ten nodes: MIR156, FLC, SPL9, SVP, SOC1, GA, FCA, PNY, TFL1, and AP1, was sufficient to cause differentiation from one phenotype to another one, whereas the other thirteen did not cause meristem transitions when they were conditionally downregulated in silico.

Surprisingly, FT or LFY in silico perturbation were not enough to drove the system to differentiation. Furthermore, even FT and SPL9 have the same in and out-degree (Supplementary Figure 2), in silico perturbation of SPL9, but not FT, was enough to change the phenotype in the simulation. We searched for components of the FT-GRN driving differentiation when transiently perturbed; the results showed that only input’s perturbations, but not others, were sufficient to drive differentiation (Supplementary Figure 8). This result further supports that the FT-GRN is robust, and more important, that differentiation driven by the FT-GRN occurs in response to signal sensing and integration. The developmental temporal sequence JVM-AVM-IM-FM observed in A. thaliana is driven by changes in the inputs AGE, VER, FCA, and PNY, but not by transient perturbations in the CO node. This may seem to be contradictory since A. thaliana reproductive development is induced by LD photoperiod, however, considering that it is a facultative flowering plant; it can bloom under SD conditions. The results of the FT-GRN model highlight that the developmental sequence JVM-AVM-IM-FM is robust and independent of photoperiod conditions even if reproductive development is faster under LD than SD. Therefore, developmental transitions are driven by specific signals such as aging (AGE and FCA), vernalization (VER), and spatial cues at the SAM (PNY).

Reshaping the attractors landscape associated with a GRN dynamical model by modulation of gene’s expression helps to understand the role of particular genes in development in the context of the complex GRN involved (Davila-Velderrain et al., 2015b). This approach was applied to the continuous FT-GRN model to assess how the network transduces input signals into developmental transitions (Figure 9). This analysis is like the differentiation trajectories found in the Boolean model (Supplementary Figure 8), but it is possible to visualize the change of each gene over time following the development in response to the input signals.

The inputs that drive the developmental trajectories from JVM to AVM (AGE and VER), AVM to IM (FCA), and IM to FM (PNY) were “turned on” in silico. This allowed us to observe how the information was propagated through the network as the relative gene expression levels changed until reaching a new steady state (Figure 9). Interestingly, similar expression trajectories occurred in both SD and LD development simulations for most nodes, except for FT which, as expected, was upregulated in LD, but not in SD during flowering transition. Interestingly, XAL2 was differentially expressed in response to aging between SD and LD during the vegetative phase transition. Considering that CO or other LD signal regulates XAL2 independently of SPL9 (or other component of the aging module), the model predicts that XAL2 should be expressed in the JVM earlier in LD than in SD growing plants (Figures 9B,C). Another difference between SD and LD was the dynamics of TFL1 expression to allow the FM fate acquisition. Under SD and LD, the model predicts that TFL1 is shut down but not in a linear way; the curve of TFL1 expression during the floral fate acquisition (IM to FM) in LD is smoother than under SD, where a small bump is observed before the subsequent decline (Figures 9B,C). This may imply that an additional dynamical barrier needs to be surpassed under SD conditions to attain FM identity compared to LD. Indeed, it is more likely to observe an incomplete commitment to the FM state under SD than under LD in A. thaliana and other plants (Okamuro et al., 1996; Tooke et al., 2005).

The model showed that PNY is an essential component of the network, allowing the system to converge to the floral meristem (Figures 8, 9), even so, there is not much information about its regulation. SOC1 may regulates PNY, but the evidence was inconclusive (Immink et al., 2012; Andres et al., 2015). The results obtained in this work led us to hypothesize that XAL2 or SOC1 may regulate this gene. Therefore, we tested this hypothesis and found that PNY relative mRNA accumulation in soc1 and particularly in xal2 mutants showed lower levels than WT in apices of SD growing plants (Supplementary Figure 9), suggesting that MADS-domain TF regulate this important homeobox. Hence, this FT-GRN is helpful to identify missing interactions when there is insufficient information about some gene’s regulation.

Previous work showed that XAL2 participates in flowering transition, both as a regulator of flowering time and as a regulator of meristem identity genes, particularly as a positive direct transcriptional activator of TFL1 (Pérez-Ruiz et al., 2015). Here we focused on the role of XAL2 as a flowering promoter under short days because the strongest phenotypes of xal2 mutants are observed under this growth condition (Pérez-Ruiz et al., 2015), while a previous model had only considered its role for long days. Therefore, we explored XAL2 regulatory interactions under SD and found that XAL2 is positively regulated by SPL9 and SPL15 (Figure 1), which are part of the aging module (Quiroz et al., 2021). Results showed that XAL2 is necessary for the GA regulatory response of SOC1, LFY, and AP1 under SD conditions (Figure 1). The integration of the FT-GRN model allowed us to simulate different “growth” conditions in silico for both wild-type and mutants (Figure 7). The FT-GRN model shows that XAL2 is important for GA flowering promotion response particularly under SD, in agreement with the experimental data shown here. The model also predicts a XAL2 developmental effect in the vegetative phase change since this transition is compromised in the simulated spl9 xal2 double mutant that can be validated experimentally in the future (Figure 7). Therefore, the prediction suggests that SPL9/SPL15 and XAL2 participate in several phase change transitions of the SAM. In this regard, it is known that SPL9 participates during the vegetative phase transition in A. thaliana (Xu et al., 2016).

However, new findings showed that SPL9 is not only controlled by miR156 and GAs, but other endogenous signals also participate. Brassinosteroids (BR), promotes vegetative phase transition through regulating the activity of SPL9 at transcriptional and post-transcriptional levels. Moreover, the BRASSINAZOLE-RESISTANT 1 (BZR1), a key transcription factor of BR signaling pathway, interacts with SPL9 to cooperatively regulate the expression of downstream genes (Wang et al., 2021). It is interesting, that BZR1 and BRI1-EMS-SUPPRESSOR 1 (BES1), also upregulate FT and SOC1 and repress DELLAs and FLC (Bao et al., 2019), strongly suggesting BR involvement in SAM’s phase transitions. Further studies can be done on the role of this hormone in XAL2 regulation and its role in the context of our network.

To our knowledge, this work includes, for the first time, the miR156 – SPLs – miR172 – AP2-like aging module as part of a flowering dynamic model, and the inclusion of this module in the FT-GRN was surprisingly informative. With the FT-GRN we could differentiate the juvenile and the adult vegetative states of the SAM and postulate a system-level developmental mechanism to underlie both the vegetative phase change and the transition to a reproductive state. In this sense, the FT-GRN Boolean model showed that regulatory interactions underlying both juvenile to adult vegetative phases transitions, are tightly interlinked to regulatory modules involved in the transition to inflorescence and flowering developmental stages. During the vegetative phase of development, flowering-promoting genes are prevented from being expressed until they are required. The model also suggests that flowering promoting signals accelerate the vegetative phase transition.

The simulations of mutants (Figures 6, 7) often showed that flowering inductive genes also affect the vegetative phase change, while only at maturity they trigger the reproductive phase (Telfer et al., 1997; Willmann and Poethig, 2011). For example, the model suggests that high GA levels or constitutive expression of FT may promote the vegetative phase change. There is no information on the juvenile to the adult vegetative transition of extremely fast flowering plants, like FT, SOC1, XAL2, or AGL24 overexpressor lines. However, high levels of gibberellins promote abaxial trichome production in A. thaliana, a classic trait that characterizes the adult vegetative phase (Telfer et al., 1997) and, it has been discovered that FT-signaling regulates the transcription of SWEET10, a bidirectional sucrose transporter in the leaves (Andrés et al., 2020), that could accelerate vegetative development by mobilizing carbohydrate resources. Hence, our model revealed that flowering signals and genes known to function in the reproductive transition could also participate in the vegetative phase change, and conversely, the vegetative phase change is necessary for reproductive competence.

We also found that PNY is necessary for flower meristem formation in the context of the FT-GRN network (discussed below). PNY regulation is almost unknown, although it had been described as a possible target of SOC1 (Immink et al., 2012; Andres et al., 2015). The results obtained in this work led us to hypothesize that XAL2 or SOC1 may regulate this gene. Indeed, we found that PNY mRNA accumulation is lower in xal2 and soc1 mutants than in WT apices of SD growing plants (Supplementary Figure 9). Therefore, the model has proven to be a helpful framework to integrate available experimental data, detect evidence holes and to postulate novel hypotheses or predictions that could not be proposed without such a systems biology framework.

Previously, an exploration of plasticity by stochastic simulated perturbations of a floral organ specification (FOS) Boolean GRN model (Alvarez-Buylla et al., 2008) suggested that the inflorescence meristem hardly transited into other floral organ basins when low noise levels were used. When larger magnitudes of noise were simulated, the system directly jumped to the attractors corresponding to gene expression configurations of stamens or carpels (Alvarez-Buylla et al., 2008). We postulated that several components and interactions linking SAM states and floral organ attractors could be missing in the FOS model, and this could hamper a more detailed and joint exploration of the attractors landscape. In A. thaliana, floral meristems develop at the flanks of the SAM and only some mutations enable the SAM to differentiate into a terminal flower (Shannon and Meeks-Wagner, 1991; Mandel and Yanofsky, 1995; Kardailsky et al., 1999; Kobayashi et al., 1999). On the other hand, the model of Dinh et al. (2017) stressed the importance of spatial cues at the IM to establish the identity of lateral organs through explicit references to this aspect, by the inclusion of an “apex” node and an auxin-maxima node. Certainly, hormones such as auxin (Li et al., 2013), cytokinins (Denay et al., 2017), and BR (Hepworth and Pautot, 2015) have been implicated in FM initiation, formation and separation from the IM. Also, other signals like reactive oxygen species have a role in flowering transition and meristem maturation (Huang et al., 2021). Because FM forms at the flanks of the IM in a specific spatial distribution, coupling GRN dynamics with different physicochemical fields including those signaling cues, together with mechano-elastic forces emergent from cellular dynamics and meristem structure could explain FM initiation and separation from the IM.

In our model the PNY node drives the system from the SAM vegetative and inflorescence attractors to the FM attractors (Supplementary Figure 8 and Figure 9). Furthermore, explorations of the attractors landscape (Figure 8) suggest that development of floral meristems requires additional interactions involving PNY. The PNY/PNF-SHOOT MERISTEMLESS (STM) complexes (Cole et al., 2006; Rutjens et al., 2009) regulate maintenance of the SAM and they specify FM emergence in reproductive development by maintaining a boundary between the central and peripheral zones at the SAM and restricting organogenesis to the periphery of the apical meristem (Ung et al., 2011; Khan et al., 2015). The pny pnf plants initiate compact shoots that fail to form flowers even though a subset of inflorescence meristem identity genes are expressed in the SAM (Kanrar et al., 2008; Lal et al., 2011; Andres et al., 2015). PNY/PNF-STM are required for induction of the floral identity genes AP1 and LFY by the FT-FD-SPL3/4/5 or by SOC1-AGL24 complexes during inflorescence development (Kanrar et al., 2008; Lal et al., 2011; Smith et al., 2011; Andres et al., 2015). Thus, it is probable that the PNY node, which represents both PNY and PNF genes in the model, connects the FT-GRN to other regulatory modules that include spatial cues and cellular dynamics at the inflorescence apex. Conversely, there might be feedback mechanisms between the flowering genes and PNY as suggested by our experimental data.

Systems biology approaches are required to provide mechanistic understanding of development (Albert and Thakar, 2014; Davila-Velderrain et al., 2015a). In this work, we have used a GRN Boolean model to address the system-level mechanisms underlying A. thaliana reproductive development. Previous models captured important aspects of flowering time under LD conditions, and under the assumption of flowering pathways and integrator genes that transfer information to AP1 as a quantitative output (Jaeger et al., 2013; Wang et al., 2014; Leal Valentim et al., 2015; van Dijk and Molenaar, 2017; Haspolat et al., 2019). The latter models relied on continuous variables optimized to calculate flowering time, which was modeled as a variable dependent on AP1 relative expression levels. Therefore, those models have helped us understand important system-level aspects of flowering time such as: the commitment and irreversibility of development, integration of noisy signals, the importance of cooperativity and feedback loops, and prediction of missing interactions (Jaeger et al., 2013; Leal Valentim et al., 2015).

Here, we propose a different modeling framework to previous modeling approaches, that although being qualitative, it integrates experimental data on a larger number of regulatory components and enables formal analyses of the role of environmental cues in flower transition. Even if it does not capture quantitative differences in flowering time, it constitutes a dynamical, multi-stable system that integrates previously identified genetic flowering transition pathways to explain SAM phase transitions in both vegetative and reproductive development. Furthermore, the robustness analyses performed showed that the FT-GRN is a robust biological module, and the recovered attractors persist under different initial conditions, even when kinetic functions were altered to produce topological changes.

The model dynamically canalizes external and internal signals into the JVM – AVM – IM – FM developmental sequence as it is qualitatively observed in plants (Figure 9). This suggests that the dynamics of the underlying GRN establish developmental checkpoints, including the commitment to the reproductive phase. We showed that the same regulatory network explain how flowering transition is induced in response to internal cues modulated by environmental signals. The comparison between SD and LD developmental trajectories suggests that endogenous signals are necessary to overcome flowering repressive signals (Figure 9), even when long-day photoperiod accelerates this transition (Quiroz et al., 2021).

Development emerges from the complex interaction between growth and differentiation, and it is modulated by environmental signals. Hence, development is the outcome of coupled GRN dynamics, spatial cues and mechano-elastic fields or constrains without the need of a central “choreographer” (Barrio et al., 2010). Here, we focused on multi-stability as an emergent property of the FT-GRN that explains the attainment of different types of meristems at the apex, and differentiation in response to endogenous and environmental signals during plant development from vegetative to reproductive states. The FT-GRN constitutes an important step toward the systems biology understanding of A. thaliana reproductive development, and it can be coupled with cellular dynamics of a growing domain and with explicit spatial components in the future to further understand additional morphogenetic patterns and behaviors.

Although the present FT-GRN model recapitulates phase transitions of the SAM, it has limitations. It does not include some signals known to participate in the reproductive phase transition, such as changes in temperature (Susila et al., 2018) trehalose-6-phosphate signaling (Wahl et al., 2013), other hormones besides GA, or redox status (Huang et al., 2021).

Also, we did not include chromatin remodeling factors, implicated in SAM’s development (You et al., 2017), which may add another level of regulatory complexity to the FT-GRN network. In a study that analyzed genome-wide histone modifications by the distribution of H3K27me3 and H3K4me3 marks along with expression changes in different developmental series including flowering transition in A. thaliana, it is reported that during early flower morphogenesis, changes in H3K4me3 prevail over changes in H3K27me3 and quantitatively correlate with expression changes, while H3K27me3 changes occur later (Engelhorn et al., 2017). Also, they found that genes with higher levels of H3K4me3 and lower levels of H3K27me3 are more prone to coincide with expression changes. Two MADS-box genes SEPALLATA 1 and 2 (SEP1/2) have that characteristic. Furthermore, it has been proved that AP1 and SEP3 bind and help to open the chromatin before expression changes are observed at the apical meristem during reproductive transition (Engelhorn et al., 2017). Further functional analyses of these genes and their epigenetic regulation may be key to incorporate an upper-level of regulation in the FT-GRN.

Moreover, temporal patterns of the epigenetic marks are likely to vary among cell-types at the SAM. The expression patterns observed in different cell types by scRNA-seq (Zhang et al., 2021) can be related to many functions of each gene besides SAM’s phase transition. In fact, all genes detected in the aerial part of the plant are found in the root as well (Supplementary Table 6). Such complex expression patterns emphasize the need of regulatory multi-stable network models to integrate available functional data, postulate a network-level mechanistic model and validate the recovered attractors with single-cell top-down data. Nonetheless, for now we may say that many of the genes of the network expressed in both the RAM and SAM showed expression at the stem cells and/or the proliferation regions of the meristems, consistent with their role in generic meristem maintenance and cellular differentiation decisions at the root and the shoot.

Single-cell RNA sequencing is a powerful technique to be considered. The information given by this technique shows that FCA and PNY (to which there were no information at the vegetative stage) are expressed in a diverse set of cell types. Additional functional experiments are required to understand their regulation at the vegetative stage. It also helps us to postulate new hypothesis. SOC1 and AGL24 proteins heterodimerize allowing SOC1 to enter the nucleus (Lee et al., 2008). However, differential expression patterns given by scRNA-seq does not coincide completely with this idea, neither at the aerial part or the root cell types. Although further research is required, we may speculate that, similarly to FLC and SVP (Mateos et al., 2015), SOC1 and AGL24 may have joint and independent functions, at least for those vegetative stages. It is known that some of the proteins are able to migrate and probably the same happens between cell types. Thus, the expression results need to be taken with caution. Nevertheless, it would be important to compare the data from single cell transcriptomic analyses from inflorescence meristem with the data of the vegetative phase to find out if those expression patterns in each cell-type persist or if they change along development.

Additional new findings may be relevant to complement the network to understand SAM’s phase transitions. For example, it has been discovered that FLC/SVP bind and represses TARGET OF FLC AND SVP 1 (TFS1), while SOC1/SPL9 upregulates it in the presence of different chromatin remodeling factors (Richter et al., 2019; Madrid et al., 2021). Additionally, it has been reported that GA signaling participates in flowering transition under lower mild temperatures, independently of FT regulation (Galvão et al., 2015), and also by repression of flowering repressors via GAI ASSOCIATED FACTOR 1 (GAF1; Fukazawa et al., 2021). The present model can be modified as new interactions are updated. Finally, it would be interesting to include genes implicated in stem cell maintenance and the genes that regulate them during meristem termination, for example, WUSCHEL (WUS), CLAVATA (CLV), and AGAMOUS (AG; Landau et al., 2015).