Gata6, Nanog and Erk signaling control cell fate in the inner cell mass through a tristable regulatory network

In the mouse, two differentiation processes take place before the implantation of the egg in the uterus. The first one gives rise to the inner cell mass (ICM) and the trophoblast (TE). The second one is the differentiation of the ICM into primitive endoderm (PrE) and epiblast (Epi). Two antagonistic transcription factors control the differentiation of the ICM into Epi and PrE: Nanog is required for the differentiation into Epi cells (Mitsui et al., 2003; Silva et al., 2009; Messerschmidt and Kemler, 2010; Frankenberg et al., 2011), whereas Gata6 is necessary to produce the PrE epithelium in vitro and in vivo (Morrisey et al., 1998; Koutsourakis et al., 1999; Capo-Chichi et al., 2005; Cai et al., 2008; Morris et al., 2010). The zygotic expression of these genes starts around the 2/4-cell stage (Guo et al., 2010; Miyanari and Torres-Padilla, 2012), and from the 8-cell [embryonic day (E) 2.5] to the 32-cell (E3.0), stage, Gata6 and Nanog proteins accumulate in almost all the cells (Dietrich and Hiiragi, 2007; Plusa et al., 2008). From E3.0-E3.25, their expression becomes mutually exclusive asynchronously within the ICM cells. Hence, at E3.75, the ICM contains two distinct cell populations that have a salt-and-pepper pattern: Gata6-expressing PrE progenitors and Nanog-expressing Epi progenitors (Rossant et al., 2003; Chazaud et al., 2006; Kurimoto et al., 2006; Plusa et al., 2008; Guo et al., 2010). These two populations are then sorted, so that the PrE forms a layer of cells separating the Epi from the blastocoel (Rula et al., 2007; Plusa et al., 2008; Meilhac et al., 2009). After specification, PrE progenitors activate several tissue-specific genes, such as Pdgfra, Sox17, Gata4, Dab2 and Lrp2, which are required for their maturation (Stephenson et al., 2012; Artus and Chazaud, 2014).

Experimental findings indicate that Nanog and Gata6 inhibit each other's expression. First, the invalidation of Nanog induces the expression of Gata6 in the whole ICM (Frankenberg et al., 2011), while forced expression of Gata6 in ES cells downregulates Nanog and pluripotency markers (Fujikura et al., 2002; Shimosato et al., 2007). Moreover, Nanog can bind to Gata6 promoter and directly decreases its activity in vitro (Singh et al., 2007).

Besides the Nanog and Gata6 network of interactions, the Fgf/RTK signaling pathway also plays a crucial role in the balance between Epi and PrE cell fate specification. Embryos mutant for Grb2 – an adaptor of the Erk signaling pathway – do not produce any PrE cells, whereas all ICM cells express Nanog (Chazaud et al., 2006). Likewise, culturing wild-type embryos with a Mek inhibitor abolishes the expression of Gata6 and induces Nanog expression (Nichols et al., 2009; Yamanaka et al., 2010). Conversely, if these embryos are cultured with recombinant Fgf4, they present a larger proportion of cells differentiating into PrE (Yamanaka et al., 2010). Interestingly, there is a window of plasticity between E2.5 and E4.0 where ICM cells can change their identity through the influence of their Fgf/RTK environment (Yamanaka et al., 2010; Grabarek et al., 2012; Arias et al., 2013). Inhibiting the Erk signaling pathway also prevents ES cell differentiation into PrE and maintains them in a pluripotent state (Cheng et al., 1998; Burdon et al., 1999; Hamazaki et al., 2006; Ying et al., 2008). Experiments modulating the Fgfr/Erk pathway in Nanog mutants revealed that, in a first phase around E2.5, Gata6 expression is induced by the Erk pathway. Afterwards, Erk signaling progressively becomes dispensable for the maintenance of Gata6 expression in the absence of Nanog, but remains necessary to counteract the Nanog-induced Gata6 repression (Frankenberg et al., 2011). Thus, in this second phase, the Fgfr/Erk signaling pathway indirectly activates Gata6 through Nanog downregulation.

The analysis of Fgf4 mutants shows that this ligand is required for PrE differentiation. Although Fgf4 is not required to induce Gata6 expression at E2.5, it is necessary for its maintenance after the 32-cell stage to drive the cells towards a PrE fate (Feldman et al., 1995; Arman et al., 1998; Kang et al., 2013; Krawchuk et al., 2013; Ohnishi et al., 2014). Thus, another Fgf or RTK ligand must be present around E2.5 to induce Gata6 expression.

The examination of Nanog mutant embryos uncovered a non-cell-autonomous role for Fgf4 in the maturation of the PrE (Messerschmidt and Kemler, 2010; Frankenberg et al., 2011). Indeed, secretion of Fgf4 from Epi cells, stimulated by Nanog, induces the expression of PrE maturation markers, such as Sox17, Gata4 or Pdgfra, that are downstream of PrE specification. PrE maturation is also disturbed in Oct4^−/− embryos, whereas Epi versus PrE specification seems to occur correctly (Frum et al., 2013; Le Bin et al., 2014). Previous studies on the phenotype of Gata6^−/− embryos have shown that the PrE epithelium is not produced at E4.5 (Cai et al., 2008). However, analyses at earlier stages discriminating between a failure to specify or a failure to differentiate cell lineages have not been described.

In the present study, we further investigated the interactions between Fgfr/Erk signaling and the transcription factors Nanog and Gata6. We first analyzed the early phenotype of Gata6^−/− embryos, demonstrating that this factor is required for PrE specification and for the inhibition of an Epi fate. Then, to investigate in more detail the mechanism of cell fate specification, we built a mathematical model describing the gene regulatory network responsible for ICM differentiation, including the effects of Fgf/RTK signaling. This model, based on previously reported in vivo and in vitro experimental results on the interplay between Nanog, Gata6 and RTK signaling in the early mouse embryo, is the first one to propose a self-organized mechanism for the PrE versus Epi fate choice. We use the model to make predictions on Gata6^−/− and Nanog^−/− mutants and verify them experimentally. This interdisciplinary approach also allowed us to decipher the imbalanced and precocious Epi specification observed in Gata6^+/− embryos.

We analyzed Gata6^−/− embryos during Epi and PrE lineage specification from E3.25 to E4.5. A failure to produce the PrE epithelium reported by Cai et al. (2008) could be due either to a lack of PrE specification, shown by a conversion of all ICM cells into Epi, or to a failure of PrE maturation and differentiation. In the latter case, PrE precursors would be specified, but fail to express markers of PrE maturation or die, reducing the number of ICM cells. Litters from Gata6^+/− intercrosses produce Gata6^−/− embryos at Mendelian ratios at E3.75. Although wild-type embryos produce around 45% of PrE and 53% of Epi, all ICM cells express Nanog in the Gata6^−/− embryos (n=6) (Fig. 1A, see also Fig. 5D for quantification). Moreover, no PrE marker, such as Sox17 and Gata4 (n=8), are expressed in these mutant embryos. As ICM cell numbers are alike in wild-type and mutant embryos (see Fig. 4C), this experiment shows that Gata6^−/− mutants cannot specify PrE cells and that all ICM cells adopt an Epi fate.

It was previously shown that the RTK pathway is required to induce Gata6 expression and thus specify PrE precursors in wild-type embryos (Nichols et al., 2009; Yamanaka et al., 2010; Frankenberg et al., 2011). Alternatively, the RTK pathway could also act in parallel to specify PrE independently of Gata6. To test this hypothesis, we administered recombinant Fgf4 during embryo cultures at different stages of development. With wild-type embryos, Fgf4 treatment is able to induce Sox17 and PrE markers in all ICM cells as previously shown (Yamanaka et al., 2010). Conversely, Fgf4 is not able to rescue Sox17 or Gata4 expression at any time-window tested in the mutant embryos (n=14) (Fig. 1B,C; supplementary material Fig. S1A), meaning that Fgf4 alone cannot induce PrE specification. As we had previously shown that PrE maturation requires Fgf4 (Frankenberg et al., 2011), these analyses demonstrate that the expression of Sox17 and Gata4 requires an activation of both the Gata6 and the Fgf4 pathway.

In parallel to PrE markers, we analyzed Nanog expression after Fgf4 treatment. Nanog expression is inhibited by early Fgf4 treatments (from the 8-cell stage) in wild-type and Gata6 mutant embryos (n=5; supplementary material Fig. S1B). This shows that the RTK pathway can inhibit Nanog expression independently of Gata6 at this stage. By contrast, when we applied Fgf4 from E3.25 to E3.75, early during the salt-and-pepper set-up, Nanog expression was maintained in all Gata6^−/− ICM cells (n=4; Fig. 1C). Therefore, at this stage Nanog expression is insensitive to repression by Fgf4 in the Gata6^−/− embryos. Thus, as observed for Gata6 in Nanog mutants (Frankenberg et al., 2011), we can consider two phases for Nanog expression: phase 1, when Nanog is sensitive to the Fgf/RTK pathway inhibition; and phase 2, when Nanog expression can be maintained despite the activation of the Fgf/RTK pathway.

Finally, when cultures are prolonged until E4.5, Nanog expression is lost in the Gata6^−/− embryos, whether they have been treated with Fgf4 (n=7) or with the vehicle only (n=7; Fig. 1B; supplementary material Fig. S1A,C). In wild-type embryos, Nanog expression is also absent in the Epi cells at E4.5 (supplementary material Fig. S1C) (Chambers et al., 2003). This shows that the mutation of Gata6 does not influence the downregulation of Nanog at E4.5.

All mutants analyses and epistatic studies recently carried out indicate that Gata6, Nanog and Fgfr/RTK are sufficient and required to control ICM cell specification (Chazaud et al., 2006; Nichols et al., 2009; Yamanaka et al., 2010; Frankenberg et al., 2011; Kang et al., 2013; Krawchuk et al., 2013). Although these results bear significance to ICM cell plasticity, biological tools are not sufficient to understand the inter- and intracellular molecular dynamics that allow the Epi and PrE specification within the ICM. Mathematical modeling provides a useful approach for investigating in further detail the complexity of this developmental process. Based primarily on in vivo and in vitro results, we developed a model for cell fate specification of ICM into Epi or PrE cells, schematized in Fig. 2. At the core of the network, Gata6 and Nanog inhibit each other and activate their own expression in vitro (Molkentin et al., 2000; Boyer et al., 2005; Singh et al., 2007; Kim et al., 2008; Verzi et al., 2010). Additional support for these regulations comes from our observation that, by ectopically expressing Gata6 in F9 cells in the presence of RTK inhibitors, Gata6 inhibits Nanog expression independently of the Fgf/RTK pathway (supplementary material Fig. S2A). The Nanog and Gata6 auto-activation loops could be direct, as both proteins can bind their own regulatory sequences (Loh et al., 2006; Verzi et al., 2010), or indirect, acting through the networks of pluripotency or of PrE differentiation, respectively (Boyer et al., 2005; Loh et al., 2006; Artus and Chazaud, 2014). Recent publications have shown that Nanog controls its expression through an autorepression in ES cells (Fidalgo et al., 2012; Navarro et al., 2012). This autorepression does not seem to occur in the embryo, at least during the preimplantation stages, as both Nanog mRNA and protein can be detected simultaneously in Epi cells (supplementary material Fig. S2B).

Besides the interactions between Gata6 and Nanog, the model incorporates the role of the Fgfr/Erk signaling pathway. This pathway, activated through the binding of Fgf4 to the receptor Fgfr2, enhances Gata6 synthesis while repressing Nanog expression (Hamazaki et al., 2006; Santostefano et al., 2012; Kang et al., 2013; Krawchuk et al., 2013; Ohnishi et al., 2014). Finally, the model takes into account the downregulation of Fgfr2 by Nanog in ES cells, and its upregulation downstream of Gata6, as suggested by ChIP experiments (Niakan et al., 2010; Ma et al., 2011). This assumption is also reinforced by single-cell qRT-PCR analyses where the Fgfr2/Fgf4 level is correlated with the Gata6/Nanog expression profile (Kurimoto et al., 2006; Guo et al., 2010; Ohnishi et al., 2014).

The differentiation status of a single cell is defined in the model by the values of four intracellular variables whose temporal dynamics are described by a set of four ordinary differential equations. The first three variables represent the level of a protein: Gata6 (G), Nanog (N) and Fgfr2 (FR). Indeed, we chose to consider the concentrations of the proteins only, and not of their mRNAs, in order to keep a low number of variables and as proteins are the final species regulating gene expression. The fourth variable (ERK) represents the level of activity of the Fgfr/Erk signaling pathway. Its value depends on the extracellular concentration of Fgf4 (Fp). The regulatory network schematized in Fig. 2 is mathematically described by the equations given in the supplementary materials and methods. Gata6 synthesis involves two contributions: it is activated, respectively, by ERK and Gata6, and inhibited by Nanog. In a symmetrical manner, Nanog synthesis involves two contributions: it is inhibited by ERK and activated by Nanog, respectively, and inhibited by Gata6. Fgfr2 expression is repressed by Nanog and induced by Gata6. The Erk pathway is activated by Fgf4-bound receptors; this activation is reversible. The equations for G, N and FR include a linear term of decay. All regulations are represented by phenomenological Hill expressions.

The steady states of the modeled system are the solutions of the equations when the variables do not vary with time. Using an appropriate set of parameter values and a proper concentration of extracellular Fgf4, the model for a single cell accounts for the existence of three stable steady states: an ICM-like state where both Nanog and Gata6 are expressed; an Epi-like state where Nanog is expressed but Gata6 is not; and a PrE-like state where Gata6 is expressed but Nanog is not (supplementary material Fig. S3B,C). At this stage, parameter values were chosen phenomenologically to account for the existence of an Epi-like state at low Fgf4, a PrE-like state at high Fgf4 and the co-existence of both states and a third ICM-like state at intermediate concentrations of Fgf4. Analyzing the possible steady states of the model as a function of parameters such as the constants characterizing auto-activation or mutual inhibition of Gata6 and Nanog shows, however, that such tristability is not a punctual phenomenon but occurs in a sizeable parameter range outside which the system displays either mono- or bistability (supplementary material Fig. S3A).

In this work, we use this cellular model at the level of a population, modeled as a network of 25 cells arranged on a square 5×5 grid, which roughly corresponds to the ICM size at E3.75. Boundary conditions are periodic to account for the spherical structure of the blastocyst. In this population, cells interact through extracellular Fgf4, which now becomes a variable of the system (as opposed to the single cell model where it is a parameter). In each cell, an additional equation (supplementary material Eqn S5) thus describes the synthesis of Fgf4, which is activated by Nanog (Frankenberg et al., 2011). We assume that Fgf4 synthesis is immediately followed by its secretion into the extracellular medium. Thus, each cell within a population is characterized by a set of five variables (see equations in the supplementary materials and methods). The first four describe an intracellular protein level (G, N, FR) or pathway activity (ERK), whereas the fifth one corresponds to the amount of Fgf4 proteins secreted by the cell (Fs).

The concentration of Fgf4 perceived by cell i (Fp_i) corresponds to the averaged level of Fgf4 produced by the cell itself and by its four closest neighbors. This is a simple way of simulating Fgf4 diffusion, which must be taken into account as Fgf4 concentration cannot be assumed to be homogeneous given the high degree of cellular compaction in the developing embryo. Importantly, the model also includes some noise in this diffusional process in the form of a random deviation (γ_i) around the average local Fgf4 concentration (Fig. 2; supplementary material Eqn S6). This extrinsic noise is the only source of stochasticity in our deterministic simulations. Once attributed randomly for each cell at the beginning of the simulations, all γ_i values remain fixed in the course of time. The level of ERK within a cell depends on the concentration of Fgf4 it perceives (Fp_i) and on the concentration of receptor Fgfr2 (FR) at its surface. As to the values of the parameters in the population model, they were phenomenologically adjusted in order to account for the available experimental data on the time evolutions of the mRNAs of Nanog, Gata6 and Fgf4 during early embryogenesis (Yamanaka et al., 2010) and on the proportions of ICM cells differentiating into PrE or Epi from wild-type and Nanog^−/− mutant embryos in a variety of conditions (Frankenberg et al., 2011; this paper).

The dynamics of Nanog (N), Gata6 (G) (Fig. 3A) and Fgf4 (Fs) (supplementary material Fig. S4) obtained with the model correspond to the expression pattern of the respective proteins during early embryogenesis, both for a future Epi cell and for a future PrE cell (Chazaud et al., 2006; Kurimoto et al., 2006; Plusa et al., 2008; Guo et al., 2010; Frankenberg et al., 2011). When multiple stable steady states co-exist, the outcome of the system depends on the initial conditions. In the embryo, Nanog and Gata6 proteins start to be detected at the 8-cell stage and, at the 32-cell stage, they are co-expressed in almost all blastomeres (Plusa et al., 2008; Lavial et al., 2012). In the model, which involves arbitrary time units, the beginning of the simulations corresponds to the stage at which Nanog and Gata6 are null (Fig. 3A). By contrast, the initial level of ERK is elevated (supplementary material Fig. S4) because of the high initial concentration of Fgf4 in the extracellular space, in line with experimental observations (Guo et al., 2010; Krawchuk et al., 2013). The results of the model do not change if other RTK ligands contribute to the initial activation of Erk, as suggested by recent experimental results (Tang et al., 2011; Kang et al., 2013; Krawchuk et al., 2013).

With these initial conditions, all cells first reach the ICM-like state, which is reflected by an increase of Gata6 and Nanog levels (Fig. 3A). Simultaneously, the level of RTK ligands in the extracellular space decreases, according to observations on Fgf4 mRNAs levels (Guo et al., 2010). As a consequence, ERK diminishes, which forces a subset of cells to leave the ICM-like state and to reach the Epi-like state, where Gata6 expression is arrested and Nanog levels are upregulated. Hence, these cells synthesize and secrete Fgf4 at a higher rate than before. The subsequent increase in the local concentration of Fgf4 pushes their neighboring cells towards the PrE-like state, where they stop expressing Nanog (Fig. 3A; supplementary material Fig. S3). In summary, the model reproduces the emergence of ICM cells, which co-express Nanog and Gata6, as well as their specification into Nanog-expressing, Fgf4-secreting, Epi-like cells and Gata6-expressing, PrE-like cells. Furthermore, these two populations reach a random salt-and-pepper distribution at the end of the simulations, which is consistent with experimental data (Chazaud et al., 2006; Plusa et al., 2008; Frankenberg et al., 2011). Importantly, the fate of a cell is determined at the level of the population and is not imposed by its specific value of γ_i. The role of this parameter is to introduce some spatial heterogeneity in the extracellular Fgf4 concentration, allowing some cells to perceive a different concentration of Fgf4 and thus quit the ICM-like state. Interestingly, the salt-and-pepper distribution can already emerge when a single cell in the population has a different value of γ_i (e.g. when γ_i=−0.1 for the cell at the center of the grid and γ_i=0.1 elsewhere). The heterogeneity then propagates over the entire field through the interactions between neighboring cells.

With the set of parameters used in Fig. 3A, 54.6±5.4% of the cells differentiate into Nanog-expressing Epi cells (Fig. 4E), which fits well with our observations (Fig. 4B). The proportions obtained with the model can be modified with small changes in the rate of Fgf4 degradation (kdf), and the model can thus reproduce the proportions obtained with mice from various genetic backgrounds, from 40% to 55% of Epi progenitors (Batlle-Morera et al., 2008; Frankenberg et al., 2011). Moreover, a subset of ICM cells undergoes apoptosis around E3.75. This mechanism was proposed to eliminate the cells that still express both Nanog and Gata6 after E3.75 (Plusa et al., 2008; Meilhac et al., 2009; Frankenberg et al., 2011). Consistent with this possibility, 13.0±6.6% of the modeled cells continue to express Nanog and Gata6 at the end of the simulations and thus remain in an ICM-like state, even though all their neighbors have already differentiated (Fig. 3A).

The model predicts that individual cell specification is asynchronous and heterogeneous in space between ICM cells. This phenomenon has been partially observed by immunofluorescence in independent studies (Gerbe et al., 2008; Plusa et al., 2008; Lavial et al., 2012). However, these global analyses did not quantify Nanog and Gata6 cell-to-cell variation at different stages. We thus performed a detailed analysis of in vivo Nanog and Gata6 protein levels at different time points and compared them with the model (Fig. 3B-E). At E2.75 (12- to 15-cell stage), most of the cells co-express Nanog and Gata6 at relatively low and comparable levels (Fig. 3B). The results are essentially the same at E3.0 (Fig. 3C). At E3.25, the large majority of cells have increased their levels of Nanog and Gata6 expression (Fig. 3D). Interestingly, some cells express both proteins, although at different levels, whereas others seem to have started their specification by decreasing one marker much more than the other (Fig. 3D). Such an asynchrony fits with the model predictions, as the temporal dynamics of Gata6 and Nanog are highly different among the simulated cells (Fig. 3A). The progressive decrease in either Nanog or Gata6 is also in agreement with single cell RNA analyses (Guo et al., 2010; Ohnishi et al., 2014). At this stage, a bias can be observed with more cells expressing Nanog high/Gata6 low (Fig. 3D, blue and red clusters, representing 48.5% of cells) compared with Gata6 high/Nanog low (Fig. 3D, yellow cluster, representing 10.7% of cells), showing that Epi cells specify first. At E3.5, many cells are already specified in PrE or Epi, displaying low levels of either Nanog or Gata6, respectively (Fig. 3E). However, more Epi cells are specified (Fig. 3E, blue cluster), as visualized by Gata6 levels below background detection, compared with future PrE cells (Fig. 3E, yellow and purple clusters). Interestingly, this asynchrony is predicted by the model where, on average, Epi cells are specified earlier than PrE progenitors (18 versus 28 arbitrary units of time to specify Epi and PrE cells in the simulation shown in Fig. 3A).

Another method to validate the model is to analyze ICM cell plasticity depending on Fgf/RTK modulation. It is worth mentioning that, in the model, once a cell has chosen its identity (Epi or PrE), it will not change its physiological condition. As both states co-exist over a large range of extracellular Fgf4 concentrations (supplementary material Fig. S3B,C), the level of Fgf4 encountered in the simulations (of untreated wild-type embryos) is never high enough to push a cell out of the Epi-like state, nor sufficiently low to push a cell out of the PrE-like state. Hence, once the Epi and PrE progenitors are specified (i.e. the cells have reached one of these two steady states), their fate cannot be changed by modifying the concentration of Fgf4 in the range explored during the developmental process. By contrast, treatments with either Fgf4 or Fgfr/Mek inhibitors suffice to push a wild-type cell out of the Epi-like or PrE-like state, respectively, in both the model (see bifurcation diagrams shown in supplementary material Fig. S3) and the experimental data (Nichols et al., 2009; Yamanaka et al., 2010).

Interestingly, the quantity of undecided cells dramatically increases if the model is simulated with a globally or even locally homogeneous (and high) distribution of Fgf4 in the intercellular space (all γ_i=0) (data not shown). Thus, the model predicts that heterogeneities in the extracellular distribution of Fgf4 are essential for the specification process to work. This was recently confirmed by experiments of Fgf4 administration on Fgf4^−/− embryos (Kang et al., 2013; Krawchuk et al., 2013). These mutants do not produce any PrE progenitor, and the uniform administration of exogenous Fgf4 most often fails to rescue their phenotype, suggesting that local heterogeneities in Fgf4 concentration or availability are required for the emergence of the salt-and-pepper distribution, which is consistent with the results of the model.

In Gata6^−/− mutants, all the ICM cells differentiate into Epi (Fig. 1A). Accordingly, if the model is simulated with a null rate of Gata6 synthesis, the level of Nanog increases in all the cells (data not shown). Furthermore, the model reproduces the dynamics in two phases obtained when Gata6^−/− ICM are treated with exogenous Fgf4 (supplementary material Fig. S1B, Fig. 1C). If Fgf4 is added from the beginning of the simulation (t=0), Nanog levels do not increase in any cell (supplementary material Fig. S5A); by contrast, if it is added when Nanog levels have already reached a maximum, Nanog expression is maintained (supplementary material Fig. S5B). Thus, in the model, Nanog self-activation is strong enough – when its level of expression is sufficiently high – to counteract the direct inhibitory effect of the Fgf/RTK signaling pathway.

This different sensitivity to exogenous Fgf4 at the level of expression of Nanog determines the two phases. The model predicts that the transition from one phase to the other does not occur at the same time for all cells. Indeed, if Fgf4 is administered when Nanog is already expressed, but not yet at its maximal level, only a subset of cells maintains Nanog expression, showing that they are already in phase 2 (supplementary material Fig. S5C). The proportion of cells in phase 2 increases with the time of Fgf4 administration (supplementary material Fig. S5D). An earlier specification of Epi cells in Gata6^−/− mutants compared with wild type can be observed in both the model (supplementary material Fig. S5C) and experimental data (Fig. 5D), meaning that Gata6 expression delays Epi specification.

In the course of the experiments, we noticed that Gata6^+/− embryos have a large deficit in the specification of PrE cells at E3.75 (n=11; Fig. 4A,B). This defect in PrE specification is counterbalanced by the acquisition of an Epi fate (Fig. 4A-C). Therefore, removing one allele reveals a dose-sensitive defect.

As Gata6 expression is biallelic (Miyanari and Torres-Padilla, 2012), one would expect that all ICM cells of Gata6^+/− embryos behave equally, i.e. either they are able to specify PrE and so behave as in wild-type embryos (with no effect of removing one allele) or the low dose of Gata6 does not allow PrE specification, as in Gata6^−/− embryos. It is thus puzzling to observe that a few cells manage to specify into PrE. A logical explanation is provided by the mathematical analysis. Adapting the model for Gata6^+/− cells consists of decreasing the rate of Gata6 synthesis. Because of the resulting decrease in the concentration of Fgfr2 (see Fig. 2), this change induces an increase in the extracellular concentration of Fgf4 needed to push a cell on the PrE-like state. Hence, a Gata6^+/− PrE cell needs to be surrounded by more Fgf4-producing Epi cells than a wild-type PrE cell. Consequently, the proportion of PrE cells obtained with the model is lower in Gata6^+/− than in wild-type populations (Fig. 4D). The best fit with the experimental results is obtained when the rate of Gata6 synthesis is decreased by 15%: while the model is calibrated to obtain 46.4±4.7% of PrE cells in the wild-type populations, this proportion decreases to 30.1±5.4% if the rate of Gata6 synthesis is reduced by 15% (Fig. 4E), in good agreement with experimental data (Fig. 4B). To confirm this assumption, using immunofluorescence quantification we investigated whether, in vivo, Gata6 levels are indeed decreased in the Gata6^+/− embryos. As shown in Fig. 5D and supplementary material Fig. S6, Gata6 levels are lowered by 19.3% compared with wild type at E3.25 and by 32.0% at E3.75. Thus, the model rightly suggests that the lack of one Gata6 allele in the heterozygous cells is partially compensated and does not correspond to 50% loss of Gata6 activity. Altogether, the slight alteration in Gata6 level leads to a PrE/Epi imbalance due to changes in Gata6, Nanog and Fgf/RTK dynamics.

When the model is simulated with a 15% decrease of Gata6 synthesis rate, the time required for a cell to differentiate into Epi appears to be reduced (compare black arrows in Fig. 3A and in Fig. 4D). Thus, the model predicts that the specification of a Gata6^+/− cell into Epi requires less time than for a wild-type cell (13 units of time instead of 18 for the wild-type cell). We addressed this question experimentally by counting cells expressing Fgf4 – considered as an early Epi marker). The proportions of cells expressing Fgf4 are significantly higher in Gata6^+/− embryos compared with wild type (Fig. 5A,B). The time-plot of Fgf4 expression also demonstrates a precocious induction of Epi cells in the heterozygous embryos (Fig. 5C). Quantification of Nanog and Gata6 proteins shows that this precocious specification is due to an imbalance between Nanog and Gata6 levels with lower Gata6 expression (Fig. 5D; Fig. S6). These experimental data validate the model and its underlying assumptions. Thus, the relative levels of Nanog and Gata6 not only modulate Epi/PrE ratios but also control the timing of specification. Surprisingly, we did not detect any difference in the number of PrE cells at E4.5 (n=5, see also supplementary material Fig. S7C) in vivo, meaning that the embryo is able to compensate for low levels of Gata6 at this stage.

We then assessed the sensitivity of Gata6^+/− Epi cells to the administration of exogenous Fgf4, as we did in the Gata6^−/− embryos. When treated with recombinant Fgf4 from the 8-cell stage, Gata6^+/− embryos behave like wild-type embryos, with an absence of Nanog expression and the whole ICM expressing Sox17 (supplementary material Fig. S7A). Surprisingly, whereas they possess a functional Gata6 allele, only a few ICM cells expressed PrE markers in Gata6^+/− embryos treated with Fgf4 from E3.25 (Fig. 6A-C; supplementary material Fig. S7B-D). Moreover, Nanog remained expressed in many cells at E3.75. In fact, the proportions of Epi and PrE cells were similar in treated and untreated Gata6^+/− embryos (Fig. 6B). Thus, the Epi cells are already insensitive to Fgf4 at E3.25 in Gata6^+/− embryos, as observed in Gata6^−/− Epi cells, but not in the wild type. We checked whether Gata6^+/− PrE cells are also insensitive to Fgfr/Erk signaling at E3.25. When embryos are treated with Fgfr and Mek inhibitors from this stage, Gata6 expression is inhibited, demonstrating that PrE precursors cells in Gata6^+/− embryos are still sensitive to Fgfr/Erk signaling at E3.25 (n=8; Fig. 6D) and even at E4.0 (data not shown). In summary, Epi and PrE lineages become unresponsive to RTK signaling through different and independent mechanisms. Moreover, the Epi cells of Gata6^+/− embryos become unresponsive to Fgf/Erk signaling more precociously than in wild type.

The fact that Gata6^+/− embryos exhibit a more precocious loss of sensitivity to Fgf4 compared with wild type could be the consequence of the earlier specification. As insensitivity to Fgf4 cannot be induced after a long exposure (from E1.5) to high levels of Nanog and/or RTK inhibition (Yamanaka et al., 2010), the loss of sensitivity of Epi progenitors to important variations of this signaling pathway is not due to an early and long exposure to Nanog expression/RTK inhibition, but rather is induced by other mechanisms – probably involved in the maturation of the epiblast – which come into play at E4.0 and E3.25 in wild-type and Gata6^+/− embryos, respectively.

Consistently, the loss of sensitivity of wild-type and Gata6^+/− Epi cells to Fgf4 cannot be reproduced by the model. This confirms that, in these two genotypes, the insensitivity of Epi cells to Fgf4 is most probably induced by a phenomenon not included in the gene regulatory network accounting for specification. Thus, the model suggests that other mechanisms and factors act after the specification program to consolidate Epi and PrE identities.

In this study, we investigated the role of Gata6 during preimplantation development. We first demonstrated the requirement of Gata6 for PrE specification. Indeed, as in Grb2 mutants (Chazaud et al., 2006), ICM cells of Gata6^−/− embryos can adopt only an Epi fate. Moreover, we showed that the Fgf pathway, while being required to induce Sox17 and Gata4 (Frankenberg et al., 2011), cannot rescue the Gata6^−/− PrE specification defect. This result not only demonstrates that Gata6 cannot be bypassed by Fgf4 administration but also shows that Sox17 and Gata4 expression require the activity of both Gata6 and the Fgf/RTK pathway.

In a second step, we showed that the previously identified gene regulatory network involved in the Epi/PrE specification (Stephenson et al., 2012; Artus and Chazaud, 2014) can account for experimental observations through a mechanism involving tristability. The model reveals how the salt-and-pepper pattern could be first triggered by a decrease – and not an increase – in RTK pathway activity from an initially elevated level. This event would induce the Epi specification of a subset of cells by favoring Nanog expression. The proposed mechanism is corroborated by protein quantifications showing that Epi cells specify earlier than PrE cells. As a consequence of Nanog upregulation, these cells produce more Fgf4. The resulting increase in extracellular Fgf4 induces the transition of the remaining cells towards the PrE state by activating their Fgfr/Erk pathway. This scenario does not exclude the possibility that some unknown factors could be involved in either helping the decrease of RTK activity or directly promoting Nanog expression.

The model shows that Fgf4 must be heterogeneously distributed throughout the ICM to implement the salt-and-pepper pattern. This result explains the difficulty of rescuing Fgf4 mutants with exogenous Fgf4 (Kang et al., 2013; Krawchuk et al., 2013). Additionally, simulations with the model indicate that once the cells are specified they do not change identity, unless they are confronted by an artificial activation or block of RTK signaling. This means that after specification the cell identity does not fluctuate, although it remains sensitive to high (non physiological) variations in RTK signaling.

In the model, the mechanism for the transition from the undifferentiated progenitor to two differentiated cell types involves the presence of three co-existing steady states. Such tristability, observed in other models in the context of cell fate specification (Huang et al., 2007; Tian et al., 2013), is generated here by multiple positive-feedback loops, which arise from reciprocal inhibition and self-activation of Gata6 and Nanog, mutual activation of Gata6 and Fgfr2/Erk, and mutual inhibition of Nanog and Fgfr2/Erk. The choice between the Epi and the PrE fate is predominantly determined by the status of surrounding cells, through the secretion of signaling molecules such as Fgf4, which controls RTK activity. This process is entirely self-organized: starting from a situation corresponding to the 2- to 4-cell stage, the differentiation into Epi and PrE cells results spontaneously from the changes in RTK signaling associated with the cell fate specification process.

The present model is the first one that describes the specification of ICM cells into Epi and PrE cells in vivo; other models proposed for cell fate specification in ES cells are based on noise-induced transitions (Kalmar et al., 2009) or oscillations (Glauche et al., 2010). Other groups developed bistable models for the specification of embryonic stem cells into PrE progenitors in vitro (Chickarmane and Peterson, 2008; Chickarmane et al., 2012). However, the latter models cannot be used to describe the emergence of PrE progenitors in vivo. Indeed, Oct4, which plays a key role in these models, is not involved in the core regulatory network as Epi and PrE can be specified in Oct4 mutants (Frum et al., 2013; Le Bin et al., 2014). In addition, these models do not focus on the emergence of common ICM progenitors and cannot account for the self-organized specification of these progenitors into a mixed population of Epi and PrE cells.

Removing one allele of Gata6 reveals a finely tuned dose for the balance between PrE and Epi specification. The present model suggests that this mutation results in a moderate reduction of Gata6 activity when compared with the wild type. Quantification of Gata6 protein indeed demonstrates a low decrease, with 19.3% at E3.25 and 32.0% at E3.75. Various factors could explain why the activity of Gata6 is not reduced by 50% in Gata6^+/− embryos: an increase in allelic transcription; a regulation of Gata6 mRNA targeted degradation (Elatmani et al., 2011); post-translational modifications, such as phosphorylation (Adachi et al., 2008); or modulation of protein degradation through the binding of Bmi1 (Lavial et al., 2012). Altogether, the analysis of Gata6 heterozygous embryos coupled to modeling reveals that a slight diminution of one factor can greatly perturb the system, demonstrating that precise levels of Gata6 and Nanog are crucial for the balance between PrE and Epi fates, and the timing of their specification. Remarkably, Nanog^+/− embryos do not have any phenotype in the specification of Epi and PrE cells (Frankenberg et al., 2011; Miyanari and Torres-Padilla, 2012). This might be due to the fact that, as Nanog expression is monoallelic during the specification stages (from the 8-cell to early blastocyst stage) (Miyanari and Torres-Padilla, 2012; Deng et al., 2014), its level is similar in Nanog^+/− and wild-type embryos. However, the issue of Nanog monoallelic expression is still being debated (Faddah et al., 2013; Filipczyk et al., 2013).

In the whole study, the model was used as a tool to help decipher the complex behavior arising from developmental transitions associated with tristability in the dynamics of the regulatory network. We identified novel important steps for Epi/PrE specification: (1) Nanog and Gata6 are co-expressed and both increase in all blastomeres, whereas an RTK ligand is present in the extracellular medium; (2) the decrease of RTK signaling in a subset of cells induces an Epi specification; and (3) Fgf4, produced unevenly by Epi cells, induces a PrE specification in neighboring cells. The Epi/PrE specification occurs when Gata6 and Nanog fall below critical levels and is asynchronous throughout the ICM due to heterogeneity in Fgf4 distribution in the extracellular medium. Once specification is completed, additional factors, such as Gata4, Sox17 and Pdgfra, for the PrE (Stephenson et al., 2012; Artus and Chazaud, 2014) push the cells through maturation by reinforcing cell identity or by eliciting further differentiation steps.

Our study sheds new light on the debate about the mechanism of the salt-and-pepper induction: some authors favor a stochastic induction (Ohnishi et al., 2014), while others promote a position-biased system (Morris et al., 2010,, 2013), involving a differential heritage between inner cells produced by the two asymmetric cell division cycles (8- to 16- and 16- to 32-cell stage). Our model and experimental data support the alternative mechanism proposed by Krupa et al. (2014), following the observation that it is not a difference between asymmetric divisions 1 and 2 that drives an Epi or a PrE fate, respectively, but the number of inner cells present. In this scenario, expression of Fgf4 by inner cells that accumulate through cell division would provide the level of Fgf4 required to induce a PrE identity. This interpretation holds with our view that Epi cells are specified first by promoting Nanog, and thus Fgf4 expression. Interestingly, Fgfr2 mRNA seems to be depleted in inside cells compared with outside cells at the 16-cell stage (Morris et al., 2013). This could provide an explanation for a decrease in RTK activity, which is the mechanism proposed by the model to increase Nanog levels. Alternatively, the decrease in RTK activity could result from the cell-to-cell variation of ERK.

The mechanism proposed here for the PrE versus Epi lineage specification and for the origin of the salt-and-pepper pattern is of a deterministic rather than stochastic nature, even if it requires a source of heterogeneity, which is introduced in the model via a cell-to-cell variability in the concentration of extracellular Fgf4 that they perceive. This variability is measured by parameter γ_i, which is fixed in a random manner for each cell at the beginning of the simulations and does not change thereafter. This extrinsic, random source of cellular heterogeneity is required to induce some asynchrony in the rates at which cells evolve towards the Epi state in response to the decrease of RTK ligand that they perceive. As a result, only a few cells specify into Epi and therefore secrete Fgf4, pushing the remaining ones towards a PrE state. The salt-and-pepper expression pattern is thus a natural consequence of this mechanism.

Because of the relatively large distance between the steady-state branches in the bifurcation diagram shown in supplementary material Fig. S3B,C, it seems unlikely that the sole intrinsic molecular noise associated with fluctuations in Gata6 or Nanog is able to induce the salt-and-pepper pattern. However, the possibility cannot be excluded that, for other parameter values, fluctuations due to molecular noise could trigger the transition from the ICM to either Epi or PrE if the distance between the steady-state branches is reduced. We are currently investigating this possibility, keeping in mind that excessive fluctuations might hinder the initial evolution towards the ICM state. Interestingly, the model can generate the same results, including the salt-and-pepper distribution, when decreasing γ_i and introducing intercellular variability in the initial levels of Nanog and Gata6, which represents another form of extrinsic noise. This suggests that the randomness causing the necessary heterogeneity between cells could additionally rely on different levels of expression of the components of the control network, arising from cell-to-cell variability in gene expression in these cells, as observed by Ohnishi et al. (2014).

Altogether, our mathematical and in vivo data unravel novel aspects of the mechanism governing preimplantation development in terms of the coexistence between three stable steady states corresponding to ICM, Epi and PrE, respectively. Our results reveal that it is the induction of Epi cells first, possibly through a decrease of RTK signaling, that is responsible for the initiation of the salt-and-pepper pattern of PrE and Epi cells, and that an alteration of the relative levels of Nanog and Gata6 in the cell fate regulatory network can instigate an earlier cell specification.

Experiments were performed in accordance with French and EU guidelines for the care and use of laboratory animals.

Gata6^+/− mice were obtained by mating Gata6^tm2.2Sad males (Sodhi et al., 2006) with Tg(Pgk1-cre)1Lni females (Lallemand et al., 1998). Litters with homozygous and heterozygous embryos were obtained by crossing Gata6^+/− mice through natural mating (see supplementary Materials and methods for genotyping and staging). All the phenotypes observed were fully penetrant in Gata6^−/− and Gata6^+/− embryos, and are observed in more than three littermates per experiment. Embryo cultures were carried out as previously published (Frankenberg et al., 2011).

Fluorescent in situ hybridization and immunostaining were performed as previously described (Chazaud et al., 2006; Gasnier et al., 2013) (see supplementary materials and methods for the list of the antibodies used). Cell counting was semi-automated with the Imaris software (Bitplane).

All embryos used for the quantification analysis were scanned by a Leica SP5 confocal with the same pinhole, laser intensity and z-section (according to Dietrisch et al., 2007). Quantification analysis was carried out with Imaris (Bitplane) coupled with a Matlab-based graphical interface (MathWorks). The mean of fluorescence intensity of each cell was divided by the gain of the photomultiplicator used for the detection and normalized by the quantification of DAPI fluorescence. Background level was defined as the average of the mean fluorescence intensities of randomly chosen cytoplasmic spots divided by the average of DAPI fluorescence.

The mathematical modeling is described in details in the supplementary materials and methods section.

Statistical tests were obtained with Prism (Graphpad Software). The normal distribution of values was verified with a Shapiro–Wilk normality test and results were analyzed by the Mann–Whitney test or by Student's t-test.

Clustering was carried out with R and XLSTAT software according to a k-mean clustering using the Hartigan and Wong algorithm. The number of groups (4 or 5) was attributed before the clustering in order to distinguish different cell identities. Each group is identified by a different color and with their corresponding percentages (group 1, black; group 2, yellow; group 3, blue; group 4, purple; group 5, red).

We thank M. Cohen-Tannoudji for providing PGK-Cre mice; J. Artus, M. Cohen-Tannoudji and V. Azuara for critical reading of the manuscript; and P. Pouchin and C. Belville for technical advice.