Signalling pathway crosstalk stimulated by L-proline drives mouse embryonic stem cells to primitive-ectoderm-like cells

Amino acids are present in the high micromolar to millimolar range in mammalian reproductive fluid (Aguilar and Reyley, 2005; Cetin et al., 2005; Harris et al., 2005) and are necessary to support normal embryo development in vivo (Van Winkle, 2001; Van Winkle et al., 2006; Bazer et al., 2015). Consequently, supplementation of culture media with selected amino acids or certain groups of amino acids can be used to improve preimplantation development (Gardner and Lane, 1993; Lane and Gardner, 1997; Harris et al., 2005). For example, L-proline is a conditionally non-essential amino acid present in tubal fluid at ∼140 µM in mice, ∼150 µM in humans, ∼100 µM in rabbits, 50-300 µM in sheep and ∼200 µM in cows (Aguilar and Reyley, 2005; Cetin et al., 2005). In vitro, L-proline improves bovine oocyte maturation rates (Bahrami et al., 2019), promotes development to the blastocyst stage in the mouse system when added during fertilisation (Treleaven et al., 2021) and improves development when added to mouse embryo culture following fertilisation (Morris et al., 2020).

Pluripotent mouse embryonic stem cells (ESCs) serve as an in vitro model of mammalian embryo development. When L-proline is added, either in purified form or as part of HepG2 conditioned medium (MEDII), it stimulates cells to transition to a second distinct pluripotent population known as early primitive ectoderm-like cells (EPL cells; Rathjen et al., 1999; Washington et al., 2010) or proline-induced cells (PiCs; Casalino et al., 2011; Comes et al., 2013; D'Aniello et al., 2015, 2017; Patriarca et al., 2021; Cermola et al., 2021; Minchiotti et al., 2022). EPL cells/PiCs are a metastable primed pluripotent population which can revert to naïve ESCs upon removal of L-proline (Rathjen et al., 1999; Washington et al., 2010; Casalino et al., 2011).

The transition to EPL cells recapitulates many of the features of the conversion of inner cell mass (ICM) cells in the 4.5 days post coitum (dpc) mouse embryo to pluripotent primitive ectoderm at ∼5.5 dpc. The primitive ectoderm is now primed to gastrulate and form the three multipotent germ layers (Snow, 1977; Washington et al., 2010; Rivera-Pérez and Hadjantonakis, 2015). The similarities include the following: EPL cells are more prone to exit pluripotency than ESCs and represent a primed or partially primed pluripotent population (Smith, 2017; Hoogland and Marks, 2021; Wang et al., 2021; Glover et al., 2022); expression of the ICM marker Rex1 (also known as Zfp42) in EPL cells is reduced; and there is increased expression of the primitive ectoderm markers Fgf5 and Dnmt3b (Rathjen et al., 1999; Washington et al., 2010; Glover et al., 2022). Colonies undergo a change in morphology, from round and domed to flattened monolayers, and cell-cycle time is reduced from 11 h to 8 h (Stead et al., 2002; Washington et al., 2010; Glover et al., 2022). The continued presence of L-proline in culture primes EPL cells for differentiation to neural cells through a series of embryologically-relevant intermediate cell types (Rathjen et al., 2002; Shparberg et al., 2019a,b; Cermola et al., 2021).

In ESCs, L-proline is taken up via the Snat2 (Slc38a2) transporter (Tan et al., 2011). The mechanisms by which L-proline stimulates development and the transition along the pluripotency continuum include: acute activation of signalling pathways, epigenetic remodelling and regulation of intracellular metabolism (Washington et al., 2010; Casalino et al., 2011; Comes et al., 2013; D'Aniello et al., 2015; 2017; Tan et al., 2016). Collectively, these mechanisms modify a range of emergent properties that drive developmental progression (Washington et al., 2010) and are consistent with L-proline behaving as a growth factor (Morris et al., 2020).

In the mouse embryo, the mTorc1 pathway is required for L-proline-mediated improvement in preimplantation development. L-proline also activates the Erk1/2 and Akt pathways during this time (Morris et al., 2020). When added to ESCs, L-proline acutely activates the same signalling pathways (Lonic, 2006; Washington et al., 2010), as well as the p38 pathway (Tan et al., 2016). Selective inhibition of mTorc1 (with rapamycin) or Mek1/Erk1/2 (with U0126) or P38 (with SB203580 or PP2) prevents upregulation of the EPL cell marker Dnmt3b (Lonic, 2006; Washington et al., 2010). Conversely, inhibition of the Pi3k/Akt pathway with LY294002 blocks the morphology change and increase in proliferation but allows the associated gene-expression changes to occur (Lonic, 2006). Thus, several signalling pathways are involved in the transition from ESCs to EPL cells. Selective inhibition of these pathways blocks different aspects of the transition, and collectively this shows that L-proline modulates a complex signalling network. These experiments did not comprehensively measure changes in a range of emergent properties or expression of marker genes required to better understand this complex network.

To explore this, we employed inhibitors of Mek1 (Map2k1), Fgf receptor (Fgfr), Pi3k, mTorc1 and P70-S6 kinase (S6k) individually and in combinations. The results of these factorial experiments were analysed with supervised learning and unsupervised machine learning models including kmeans clustering, multiple linear regression (MLR), MLR with interaction terms and Bayesian regularised neural network with a Gaussian prior (BRANNGP; Woolf et al., 2005; Burden and Winkler, 2008; Winkler and Burden, 2012; Epa et al., 2013). MLR with interaction terms was used to calculate synergistic and antagonistic effects. This approach is most commonly used to determine drug interactions (Sorokin et al., 2018; Julkunen et al., 2020; Panina et al., 2020) and is becoming increasingly used in stem cell biology (Chang and Zandstra, 2004; Prudhomme et al., 2004; Audet, 2010; Jakobsen et al., 2014; Ireland et al., 2020).

ESCs were maintained in either 330 or 1000 U/ml LIF then directed to transition into EPL cells by addition of 400 μM L-proline for 6 days. In addition, cells were also allowed to undergo spontaneous differentiation without added LIF (Fig. 1A). All ESCs grown in 330 or 1000 U/ml LIF maintained their dome-shaped colonies and exhibited no differences in morphology score. Colony morphology changed significantly to flattened epithelial-like colonies in the presence of L-proline. Cells allowed to differentiate spontaneously underwent a more robust morphology change (Fig. 1A,B), consistent with these cells undergoing differentiation beyond the EPL cell stage (Tan et al., 2011; Minchiotti et al., 2022).

Cell number, apoptosis and proliferation were quantified at days 2, 4 and 6 of culture. Cell number was normalised to cells growing in 1000 U/ml LIF. Cells growing in 330 U/ml LIF+400 μM L-proline increased cell number by 1.7-fold at day 2 (Fig. 1C), consistent with previous results for ESCs grown in 1000 U/ml LIF+L-proline (Washington et al., 2010). No significant changes in proliferation or apoptosis were observed (Fig. 1D,E). Cells undergoing spontaneous differentiation without LIF or L-proline were dying by day 6, with cell number reduced by 85% (Fig. 1C). Apoptosis increased by 4.2-fold at day 4 and 2.3-fold at day 6 (Fig. 1D). Proliferation decreased by 50% at day 6 (Fig. 1E), suggesting deficiencies in medium formulation and therefore a reduced capacity to support growth of differentiating cells.

After 6 days of culture, gene expression was profiled, focusing on pluripotency genes Rex1 and Oct4 (Pou5f1), primitive ectoderm markers Dnmt3b, Fgf5, Lefty2 and Otx2, and the mesendoderm marker Mixl1. There were no differences in the expression of any of the genes between ESCs grown in 330 or 1000 U/ml LIF (Fig. 1Fi). Cells grown in 330 U/ml LIF+400 μM L-proline had comparable expression of Rex1 and Oct4, indicating that maintenance of pluripotency and the expression of the mesendoderm marker Mixl1 also did not change. The expression of all primitive ectoderm markers increased (Fig. 1Fii). Cells allowed to spontaneously differentiate had a gene expression profile consistent with rapid, unregulated differentiation – reduced expression of Rex1 and Oct4 and a transient wave of expression of Dnmt3b, Fgf5 and Otx2. The expression of Lefty2 remained strongly reduced throughout, whereas Mixl1 expression increased significantly (Fig. 1Fiii).

We examined the phosphorylation status of signalling pathway intermediates drawn from the Stat3, Fgf, Mek1/Erk1/2, Pi3k/Akt, mTor, p38 and Pkc pathways (Fig. 2A), each of which is known to play a role in pluripotency and cell state transition of ESCs (Kunath et al., 2007; Lanner and Rossant, 2010; Washington et al., 2010; Cherepkova et al., 2016; Tan et al., 2016).

Naïve, self-renewing ESCs were switched from 1000 U/ml LIF to EPL cell medium containing 330 U/ml LIF+400 μM L-proline and the phosphorylation status of the pathway intermediates was quantified by western blot over the short (0-12 h) and long term (1-6 days). Only phosphorylation of Fgfr increased modestly but significantly from 12 h onwards (Fig. 2B; Fig. S1A). The P38 and Pkc pathways, including the downstream Hsp27, previously shown to be altered by addition of L-proline (Tan et al., 2016), were not detected in ESCs±L-proline under our conditions (Fig. S1B). After 4 h of serum starvation [Dulbecco's Modified Eagle Medium (DMEM)/0.1% foetal bovine serum (FBS)/β-mercaptoethanol (β-Me) only], treatment with 1 mM L-proline acutely increased phosphorylation of Erk1/2 and Stat3^Y705 (Fig. 2C; Fig. S1C).

The effect of signalling pathway inhibitors, individually or in combination, on L-proline-mediated pathway activity was assessed using the following: (1) Mapk pathway using the Mek1/2 inhibitor U0126 (U) (Favata et al., 1998); (2) Fgfr pathway using receptor antagonist SU5402 (S) (Mohammadi et al., 1997); (3) Pi3k/Akt pathway using the Pi3k inhibitor LY294002 (L) (Vlahos et al., 1994); (4) mTorc1 pathway using the mTorc1 complex inhibitor rapamycin (R) (Sabers et al., 1995); (5) mTorc1 pathway at the downstream kinase, S6k, using the inhibitor PF-4708671 (P) (Pearce et al., 2010).

ESCs were cultured under standard conditions (1000 U/ml LIF). At 0 h, 400 μM L-proline was added, along with the first inhibitor (1i). Second doses (2i) of an inhibitor were added 2 h before sample collection. An initial dose of U0126 suppressed L-proline-mediated Erk1/2 (downstream of Mek1) phosphorylation only for 6 h, after which phosphorylation returned to the level seen with L-proline only (Fig. 2D,E; Fig. S1D,E). Furthermore, any second dose of U0126 failed to suppress Erk1/2 phosphorylation for longer than 6 h (Fig. 2E; Fig. S1D,E), suggesting that Erk1/2 phosphorylation was no longer controlled by Mek1.

None of the other four pathway inhibitors suppressed Erk1/2 phosphorylation on its own (Fig. 2E; Fig. S1D,E). However, extended suppression of Erk1/2 phosphorylation occurred up to 10 h when SU5401 or LY294002 or PF-4708671 were added after the initial addition of U0126 (Fig. 2E; Fig. S1D,E), suggesting that Erk1/2 phosphorylation was now controlled by crosstalk involving the Fgfr/Pi3k/S6k axis.

Rps6 phosphorylation showed similar interplay between pathways. Rps6 phosphorylation was suppressed for up to 10 h by rapamycin and up to 8 h by LY294002 (Fig. 2E; Fig. S1D,E), consistent with it modulating the mTorc1/S6k signalling axis (Fig. 2A). However, Rps6 phosphorylation was also temporarily suppressed by inhibition of Mek1, Fgfr and Pi3k, or combinations of Mek1 inhibition followed by Fgfr inhibition, or Mek1 inhibition followed by Pi3k inhibition (Fig. 2E; Fig. S1D,E). These results suggest complex, dynamic changes in signalling pathway activity over time in the presence of L-proline.

To further elucidate pathway interactions and their effect of emergent properties, a factorial experiment was designed to monitor the ESC-to-EPL transition over 6 days in the presence of all possible combinations of the five inhibitors (Fig. 3A). On days 2 and 4, cells were counted during replating, and samples were collected to quantify apoptosis and cell proliferation by flow cytometry. At day 6, in addition to measurements of apoptosis and cell proliferation, cells were imaged for colony morphology and qPCR was used to quantify changes in the expression of marker genes.

At day 2, the eight conditions containing both LY294002 and rapamycin had poor viability. Cell number had reduced by 90% (Fig. S3), proliferation reduced by nearly 40% (Fig. S4) and apoptosis had increased by more than 50% (Fig. S5). By day 4, very few live cells were present. These conditions were considered to be non-viable and were not considered in the day 4 and 6 measurements.

To get a broader understanding of relationships between emergent properties and gene expression, a correlation matrix was generated using all data from all viable inhibitor combinations (Fig. S6A). Expected correlations were observed, such as: (1) a positive correlation between cell number and proliferation (r≤0.57); (2) a negative correlation between cell number and apoptosis (r≥−0.43) across the 6 days of transition to EPL cells; and (3) the coupling of expression of pluripotency markers Oct4 and Rex1 (r=0.80) and EPL-cell markers Dnmt3b and Fgf5 (r=0.52). However, more nuanced correlations were observed, including positive correlations between: (1) proliferation and morphology (r=0.31 at day 4 and 0.41 at day 6); and (2) proliferation (at days 4 and 6) and expression of primed pluripotency genes: Dnmt3b (r≤0.49), Fgf5 (r≤0.32) and Mixl1 (r≤0.50). A high correlation was also observed between expression of the primitive ectoderm marker Lefty2 and the pluripotency markers Oct4 (r=0.55) and Rex1 (r=0.49).

These relationships were further explored using principal component analysis (PCA) with kmeans clustering. Two clusters were identified based on the silhouette score (Fig. S7A). The first cluster included the EPL cell control (Fig. S7B) and contained the highest proportion of conditions with LY294002 (Fig. S7C). This cluster most closely resembled EPL cells, including a higher morphology score, cell number and proliferative rate, and higher expression of primed pluripotency genes Dnmt3b, Fgf5 and Lefty2 (Fig. S7D). The second cluster appeared to be more ESC-like and included a high proportion of conditions with SU5402 and rapamycin (Fig. S7C).

To understand which signalling pathways or pathway combinations drive changes in gene expression and emergent properties, we employed supervised learning techniques to generate MLR and BRANNGP models (Figs 3Bi-Ei, 4Ai-Gi).

MLR models are a supervised learning approach used to understand linear relationships and to assign the relative importance of factors and their interactions. MLR models are easy to interpret but are less accurate at predicting the responses to multiple inhibitors when the underlying relationships are nonlinear. The coefficients of the MLR models denote the sign and magnitude of the contribution each inhibitor makes to the responses (Figs 3Bii-Eii, 4Aii-Gii). In addition to standard MLR, we also generated MLR models with two interaction terms to determine whether inhibitors were acting independently (additive), synergistically or antagonistically (Fig. 3Biii-Eiii). For these models, the coefficients for the inhibitors alone were compared with those for the interaction term. An additive effect occurs when the response to multiple inhibitors is the sum of the individual contributions of the inhibitors; i.e., no significant interaction effect. A synergistic effect occurs where the interaction terms enhance the response more than the sum of the individual contributions of the inhibitors. An antagonistic effect occurs where the interaction terms reduce the response below the sum of the individual contributions of the inhibitors. The most robust and predictive MLR models have higher adjusted R² and lower standard error (σ) values (Table S1). An F-test was used to determine whether MLR models with increasing numbers of parameters (e.g. a two-way versus a one-way MLR) resulted in a significantly improved fit (Table S2).

These machine learning models are used to capture nonlinearity and interactions in stimulus-response relationships. They predict the responses of cells to new combinations of factors more accurately than MLR models but are harder to interpret, so their most important use is in more accurately predicting the effect of different factors in stem cell behaviour. They are therefore complementary to linear MLR models. We used BRANNGP to assess nonlinearity and to create models without issues of overfitting and overtraining (Burden and Winkler, 2008; Winkler and Burden, 2012). This algorithm automatically optimises model complexity. Most datasets had substantially improved adjusted R² values when modelled using BRANNGP (Table S1), reinforcing the fact that nonlinearity and/or interaction were important. The average adjusted R² value for BRANNGPs was 0.64±0.07 compared with 0.48±0.07 for MLR models and 0.53±0.09 for MLR models with two-way interaction terms. Only apoptosis and proliferation models had adjusted R² values for BRANNGP models similar to those for the MLR models, suggesting that the relationship was approximately linear. The improved predictions of the neural network model for the other datasets imply that nonlinearity or interactions between factors are important for these biological endpoints. In nonlinear models, the importance of model features are local, rather than global, so they depend on where they are assessed on the response surface defined by the nonlinear model (Rengasamy, et al., 2022). Thus, the nature and magnitude of these factors cannot be as readily deconvoluted from the BRANNGP models.

The morphology dataset (Fig. S2) was used to train MLR, MLR with two-way interaction terms and BRANNGP models (Fig. 3Bi). The MLR model with no interaction terms showed that addition of SU5402 or U0126 prevented changes in colony morphology normally expected in the presence of L-proline (Fig. 3Bii), resulting in cells which retained a domed, ESC-like appearance. An F-test indicated that MLR with two-way interaction terms provided a better fit than the MLR (Table S2). This improved model showed that SU5402, U0126 and rapamycin all prevent morphology change (Fig. 3Biii). There was no significant interaction between SU5402 and U0126, indicating that this effect was largely additive. There was a significant interaction between SU5402 and rapamycin mediating morphology change. This interaction coefficient was reversed, indicating an antagonistic effect.

Modelling was generated for cell number and proliferation data from all inhibitor combinations. MLR produced robust fits for both cell number and proliferation with adjusted R² of 0.69 and 0.79 respectively (Fig. 3Ci,Di). All inhibitors significantly reduced cell number and proliferation, with rapamycin having the largest effect (Fig. 3Cii,Dii).

For both cell number and proliferation, the MLR with two-way interactions had an improved adjusted R² (0.79 and 0.83, respectively), and the F-test showed significant improvement (Table S2). The individual effects were retained for each inhibitor, but multiple interaction effects were noted: (1) U0126 and PF-4708671 were antagonistic in both cell number and proliferation models; (2) antagonism between rapamycin and PF-4708671 in the cell number model; (3) LY294002 and rapamycin were strongly synergistic for both cell number and proliferation models (Fig. 3C-Diii). An alternate model, which attempts to overcome non-normality of the proliferation input data by dividing the data into octiles, exhibited very similar results to the MLR model (Fig. S8).

Apoptosis data generated an adequate fit using the MLR (adjusted R² of 0.38), but this was significantly increased using the MLR with two-way interactions (adjusted R² of 0.58 and F-test P<0.05; Fig. 3Ei; Table S2). In the MLR model without interaction terms, SU5402 reduced apoptosis and LY294002 and rapamycin each increased apoptosis (Fig. 3Eii). In the MLR with two-way interactions, the individual effects for LY294002 and rapamycin were lost, and instead there was a strong synergistic effect. The day 2 parameter was also significant in the model (Fig. 3Eiii). In conjunction with the reduction in proliferation, this explains the early cellular lethality of the combination of LY294002 and rapamycin, which does not occur when the inhibitors are used individually.

The apoptosis model highlighted differential effects on apoptosis between SU5402 and PF-4708671 (Fig. 3Eiii). The increase in apoptosis in PF-4708671-treated cells combined with the decrease in proliferation (Fig. 3Diii) underpins the net decrease in cell number (Fig. 3Ciii). Contrarily, SU5402-treated cells had reduced apoptosis combined with decreased proliferation (Fig. 3Diii) resulting in a net decrease cell number (Fig. 3Ciii), indicating that proliferation was the main driver of decreased cell number. U0126 did not affect apoptosis in the presence of L-proline (Fig. 3Eii), indicating that the decrease in cell number elicited by U0126 is due entirely to a decrease in proliferation (Fig. 3Cii,Dii).

Computational modelling was also used to assess how inhibitors impacted gene expression at day 6 (Fig. 4; Fig. S9). The MLR model for Dnmt3b, Fgf5 and Lefty2 had a robust fit (adjusted R²>0.5, Fig. 4Ci-Ei). Models for Rex1, Oct4 and Otx2 expression had a modest fit (adjusted R² of 0.24, 0.21 and 0.31, respectively; Fig. 4Ai,Bi,Fi). None of these models was improved significantly by adding two-way interaction effects (F>0.05; Table S2; Fig. S10).

The MLR model showed that U0126 decreased expression of the pluripotency genes Rex1 and Oct4 (Fig. 4Aii,Bii). The decrease was modest when compared with the strong downregulation of these genes for cells undergoing spontaneous differentiation (Fig. 1Fiii).

Inconsistent changes in expression patterns occurred frequently for the EPL markers Dnmt3b, Fgf5, Lefty2 and Otx2. For example, SU5402 decreased expression of Dnmt3b and Fgf5 (Fig. 4Cii,Dii) but increased expression of Lefty2 and Otx2 (Fig. 4Eii,Fii). LY294002 increased Dnmt3b expression (Fig. 4Cii) but reduced that of Lefty2 (Fig. 4Eii). Rapamycin decreased Dnmt3b, Fgf5 and Otx2 expression (Fig. 4Cii,Dii,Fii) but did not alter that of Lefty2. PF-4708671 only reduced Lefty2 expression (Fig. 4Eii). These results highlight that expression of individual genes associated with the identity of EPL cells are associated with complex signalling pathway control. In all cases, two or three of the pathways we probed regulated expression of each gene. The MLR models for the mesendoderm gene Mixl1 were not sufficiently statistically significant to be able to make strong biological statements about this gene (adjusted R² of 0.18; Fig. 4Gi).

The complexity of how inhibitors drive expression of EPL marker genes was further exemplified by data showing that many inhibitor combinations suppressed some but not all primitive ectoderm genes (Fig. S9). We ran a functional assay to identify the pluripotency capacity of each inhibitor combination. After 6 days of culture in L-proline and the various inhibitors, cells were allowed to spontaneously differentiate as embryoid bodies (EBs). Samples were collected on days 2, 3 and 4 and qRT-PCR was used to quantify expression of the primitive streak marker brachyury (T). In the absence of inhibitors, more naïve cells – like ESCs – took 4 days to upregulate T expression, whereas more primed cells – like EPL cells – upregulated expression of T at day 2 (Fig. S11A). Across all inhibitor experiments, conditions that contained U0126 or LY294002 tended to upregulate T expression earlier, and conditions that contained rapamycin tended to upregulate T later (Fig. S11Bii). We also assessed the correlations between the slope of T upregulation and other genes. Significant positive correlations were found between T upregulation and Dnmt3b, Fgf5 and Mixl1, but not the pluripotency markers Rex1 and Oct4 or the more recently adopted primitive ectoderm markers Lefty2 and Otx2 (Fig. S6B).

The small molecule inhibitors used in this study helped elucidate the role of various signalling pathways that mediate self-renewal, cell state transitions, differentiation and other emergent properties such as colony morphology, cell number, proliferation and apoptosis during the transition from ESCs to EPL cells (Fig. 5A), namely the Mapk (using the Mek1 inhibitor U0126), Fgfr (using the antagonist SU5402), Pi3k (using the Pi3k inhibitor LY294002) and mTor (using the mTorc1 complex inhibitor rapamycin, or the S6k inhibitor PF-4708671) pathways. These signalling pathways are acutely activated by L-proline (Fig. 2B,C) or have been previously associated with L-proline-mediated formation of EPL cells (Lonic, 2006; Washington et al., 2010; Tan et al., 2016).

In the absence of inhibitors, L-proline increased pathway phosphorylation (Fig. 5B), including acute phosphorylation of Stat3^Y705 and Erk1/2 within 10 min (Fig. 2C; Fig. S1C). This suggests that L-proline can rapidly induce changes in pathways known to be important for maintenance or loss of pluripotency (Stavridis et al., 2010; Huang et al., 2014). Over 6 days, Fgfr phosphorylation increased but there was no change in the canonical intermediate Erk1/2 (Fig. 2B; Fig. S1A), suggesting that Fgfr is signalling through other intermediates such as Pkc, Pi3k, Src, Stat1, P38 and Jnk (Dailey et al., 2005).

When signalling pathway inhibitors were used in the presence of L-proline, signalling pathway crosstalk led to maintenance of Mapk signalling: Erk1/2, immediately downstream of U0126 target Mek1, had decreased phosphorylation in the presence of this inhibitor but only for 6 h. A second dose of U0126 did not extend this time (Fig. 2E; Fig. S1D,E). This effect is not unique to U0126. Erk1/2 phosphorylation is only transiently reduced when a variety of Mek1 inhibitors (PD98059, PD184352, PD0325901 and U0126) are added to the culture medium of ESCs (Chen et al., 2015). Together these results argue against the loss of U0126 activity, but rather that Erk1/2 phosphorylation is maintained by pathway crosstalk which bypasses Mek1. As reduced phosphorylation of Erk1/2 in the presence of U0126 could be extended to 10 h by adding the Fgfr inhibitor SU5402 or the inhibitor LY294002, it is plausible that the Fgfr-Pi3k-Akt axis sustains L-proline-mediated phosphorylation of Erk1/2 (Dailey et al., 2005). This complex network with multiple inputs speaks to the importance of cells maintaining Erk1/2 phosphorylation to avoid widespread apoptosis, as seen in Erk1^−/−/Erk2^−/− ESCs (Chen et al., 2015).

Our factorial study assessed how signalling pathways influence a variety of properties during the L-proline-mediated transition from ESCs to EPL cells (Fig. 5). No single inhibitor was sufficient to explain all the changes in gene expression and emergent properties during ESC transition to EPL cells. Rather, our modelling suggests that these signalling pathways have discrete roles within this transition, likely supported by signalling pathway crosstalk.

The combinational experiments revealed the following information for each pathway inhibitor: (1) When the Mapk/Erk1/2 pathway was inhibited by U0126, cells did not undergo the morphological change associated with the presence of L-proline (Fig. 3Aii), despite a decrease in Rex1 and Oct4 expression (Figs 4A-B, 5). This decrease was less than in cells undergoing spontaneous differentiation (Fig. 1Fiii), indicating disruption of the pluripotency gene regulatory network (Kim et al., 2008). U0126 did not alter expression of any of the primitive ectoderm markers (Fig. 4C-F). (2) When Fgfr was inhibited by SU5402, cells did not undergo a morphology change (Fig. 3Bii) and expression of the EPL-cell markers Dnmt3b and Fgf5, which is increased in the presence of L-proline alone, was blocked. In contrast, expression of EPL-cell markers Otx2 and Lefty2 increased in the presence of this inhibitor (Figs. 4C-F, 5). This suggests that Fgfr inhibition at least partially blocks the transition. (3) When the Pi3k/Akt pathway was inhibited with LY294002, the L-proline-mediated change in colony morphology still occurred, as did the increased expression of EPL-cell markers Dnmt3b and Fgf5. The L-proline-mediated increase in Otx2 expression also occurred, but L-proline-mediated increase in Lefty2 expression was suppressed. An early increase in Lefty2 expression is associated with the transition of ESCs to EPL cells, but Lefty2 expression is low when cells undergo lineage restriction (Harvey et al., 2010). These results suggest that increased Lefty2 expression is not obligatory for the transition. (4) When mTorc1 was inhibited by rapamycin, ESCs underwent L-proline-mediated change in morphology, but increased expression of Dnmt3b, Fgf5 and Otx2 was suppressed. This suppression is consistent with previously published data (Washington et al., 2010). However, in that study rapamycin blocked the morphology change, which we did not observe. Earlier protocols generated EPL cells using 1000 U/ml LIF and L-proline, which inconsistently upregulated expression of the primitive ectoderm marker Fgf5 (Washington et al., 2010). Here, we reduced LIF to 330 U/ml, which provided robust upregulation of Fgf5 expression (Harvey et al., 2010; Glover et al., 2022). These results highlight the sensitive balance between the cytokine LIF and the growth-factor-like properties of L-proline in promoting directed cell state transitions. (5) When the S6k branch of the mTorc1 pathway was inhibited with PF-4708671, like with rapamycin, it failed to prevent the L-proline-mediated change in colony morphology. However, unlike rapamycin it did not suppress the L-proline-mediated increase in the expression of the EPL-cell markers Dnmt3b, Fgf5 and Otx2 (Figs 4C-D,F, 5). This suggests that stimulation of expression of these markers by L-proline requires the 4ebp1 (Eif4ebp1) branch of the mTorc1 pathway.

All five inhibitors reduced cell number and reduced the rate of proliferation compared with L-proline (Fig. 3Cii,Ciii,Dii,Diii) but different effects were seen on apoptosis (Fig. 3Eii,Eiii). Cells treated with LY294002 or rapamycin each significantly increased apoptosis in the MLR model, but individual effects were lost in the MLR with interaction terms in favour of the strong synergistic effect. This strong synergistic effect is also seen in the proliferation and cell number models and underpins the lack of cell viability in conditions containing LY294002 and rapamycin. These results are consistent with Pi3k being a strong mediator of cell survival and progression (Chang et al., 2003; Takahashi et al., 2005; Tsurutani et al., 2005; Yu and Cui, 2016).

The synergistic effect of LY294002 and rapamycin on cell death was not seen in interactions with other the other mTOR-associated inhibitor, PF-4708671, which blocks the S6k branch of the mTorc1 pathway. Rather, the combination of LY294002 and PF-4708671 had only an additive effect and the combination of rapamycin and PF-4708671 had and antagonistic effect. Both results support the hypothesis that mTorc1 signalling via the 4ebp1 branch is anti-apoptotic (Nawroth et al., 2011; Pons et al., 2011; Yellen et al., 2011) and pro-proliferative (Dowling et al., 2010; Nawroth et al., 2011).

Inhibition of the Mek1/Erk1/2 pathway with U0126 did not affect apoptosis, and the inhibition of Fgfr with SU5402 decreased apoptosis. These conditions had a net reduction in cell number, driven by reduced proliferation. Fgfr signalling produces cell- and state-specific effects on apoptosis and proliferation (Dailey et al., 2005), and the reduced apoptosis observed with SU5402 provides further support for this.

Collectively, these results highlight biological system complexity and make it difficult, if not impossible, a priori to determine outcomes even when a single inhibitor is used. These pathways could be better elucidated by phosphoproteomics using mass spectrometry (Li et al., 2011) to identify additional signalling pathways that may be involved, including the recently identified L-proline target Tgf-β (D'Aniello et al., 2017). These pathways can be functionally addressed using automated high throughput factorial screening for key predictive properties such as morphology or proliferation, which were shown to have high positive correlations with EPL markers Dnmt3b and Fgf5 (Fig. S6A).

We selected four primitive ectoderm markers, Dnmt3b, Fgf5, Lefty2 and Otx2, to assess how cells transitioned to EPL cells. These genes had similar expression patterns during the transition to EPL cells (Fig. 1F) but this behaviour changed under inhibitor-treated conditions (Fig. 5B).

Under standard culture conditions, L-proline-treated cells had significantly increased expression of Dnmt3b and Otx2 between days 2 and 6 and Fgf5 between days 4 and 6, but Lefty2 expression was transiently increased at days 2 and 4 (Fig. 1Fii). Two inhibitors, SU5402 and LY294002, produced gene expression profiles that were less straightforward. Cells treated with SU5402 exhibited decreased Dnmt3b and Fgf5 expression and increased Lefty2 and Otx2 expression, whereas cells treated with LY294002 had increased Dnmt3b expression but decreased Lefty2 expression (Figs 4, 5).

These differences may reflect temporal expression patterns. Our control data showed transient increase in Lefty2 expression on days 2 and 4 (Fig. 1Fii). By day 6, when we measured gene expression in inhibitor-treated cells, Lefty2 expression had returned to baseline. This is in line with previous studies using EPL cells derived from embryoid bodies cultured in MEDII, which showed transient upregulation of Lefty2 from days 1 to 4 (Harvey et al., 2010). This also explains the negative correlations between Lefty2 and morphology changes (Fig. S6A). During the 6 days of culture to EPL cells, we did not note any significant changes in expression of Mixl1 (Figs 1Fii, 4G), indicating that cells remained within the pluripotency continuum and did not form mesendoderm.

To help assess this ambiguity between markers, we employed a functional assay which measured T expression as cells underwent spontaneous differentiation (Fig. S11A). EPL cells and some inhibitor-treated conditions (i.e. LY294002, low Lefty2, and SU5402, high Lefty2) tended to upregulate T fairly early (Fig. S11B), suggesting that the cells were further along the pluripotency continuum. ESCs and some inhibitor-treated conditions (i.e. PF-4708671 and rapamycin) took the longest to upregulate expression of T, suggesting that these were the most naïve cells. Dmnt3b and Fgf5 expression at day 6 correlated highly with robust upregulation of T during the functional assay. No correlations were seen with Lefty2, suggesting that it is a poorer indicator of the position of a cell along the pluripotency continuum.

Biological complexity was further highlighted when two or more inhibitors were used together. Two-way interaction effects were used to determine whether these pathways were independent (no interaction effects), antagonistic (where blocking two pathways simultaneously leads to a lower inhibition than the sum of the two inhibitors individually) or synergistic (where blocking two pathways simultaneously leads to higher inhibition than the sum of the two inhibitors individually). Only the emergent property MLR models were improved by including interaction effects (Table S2).

Antagonistic effects were seen for Mek1 and S6k, where the combination of inhibitors U0126 and PF-4708671 attenuated the inhibition of both cell number and proliferation compared with the use of each of the inhibitors alone (Figs 3Ciii-Diii, 5). The combination of mTor and S6k inhibitors (rapamycin and PF-4708671) attenuated the inhibition of cell number (Fig. 3Ciii) and the combination of inhibitors for Fgfr and mTor (SU5402 and rapamycin) promoted a change in colony morphology that was blocked by the individual inhibitors (Fig. 3Biii). These pathways likely coalesce on common downstream intermediates or transcription factors or suppress other pathways through crosstalk (Mendoza et al., 2011; Aksamitiene et al., 2012; Wang et al., 2013; Arkun, 2016).

Addition of both LY294002 and rapamycin resulted in strong synergistic effects that reduced cell numbers and proliferation and increased apoptosis (Figs 3Ciii-Eiii, 5), resulting in non-viable cells. Both pathways, when blocked individually, reduce proliferation and increase apoptosis (Fingar et al., 2002; Jirmanova et al., 2002; Murakami et al., 2004; Gross et al., 2005), and result in large defects in cell survival when blocked in combination in T cells, glioma cells and small cell lung cancer cells (Breslin et al., 2005; Takeuchi et al., 2005; Tsurutani et al., 2005).

We have shown that L-proline activates multiple signalling pathways including the Mapk, Fgfr, Akt and mTor pathways in driving the transition of ESCs to EPL cells (Fig. 5). L-proline uptake through the Snat2 transporter likely directly alters cell signalling to modulate gene expression and emergent properties. It is also possible that other mechanisms such as metabolic flux and epigenetic changes alter the cellular landscape to facilitate autocrine cell signalling, activate signalling through crosstalk or change chromatin accessibility. This has been seen previously with autocrine Fgf4 activation of Fgfr as cells undergo lineage commitment (Kunath et al., 2007).

The L-proline-mediated transition of ESCs to EPL cells exemplifies the progression of cells from a naïve to primed or partially primed state in the pluripotency continuum (D'Aniello et al., 2017; Morgani et al., 2017; Cermola et al., 2021), and recapitulates aspects of peri- and post-implantation embryogenesis. These results are consistent with other growth factor-like roles for L-proline, including facilitating preimplantation embryo development (Morris et al., 2020; Treleaven et al., 2021) and steering differentiation of pluripotent cells towards neuroectoderm (Rathjen et al., 1999; 2002; Pelton et al., 2002; Harvey et al., 2010; Washington et al., 2010; Shparberg et al., 2019a,b). The L-proline-mediated transition along the pluripotency continuum provides a useful model to study embryonic development in vitro.

All cell culture was performed at 37°C, 5% CO₂ in a humidified incubator. D3 ESCs (Doetschman et al., 1985) were maintained in ESC self-renewal medium containing DMEM (Sigma-Aldrich), 10% FBS (AusGeneX), 1000 U/ml LIF (Neuromics), 0.1 mM β-Me (Sigma-Aldrich) and Pen/Strep (50 U/ml penicillin and 50 μg/ml streptomycin; both Sigma-Aldrich). Cells were grown as a monolayer, and passaged using Trypsin-EDTA (Sigma-Aldrich), and replated at 2000-20,000 live cells/cm² (Glover et al., 2022).

ESCs were cultured to EPL cells by growing 20,000 cells/cm² in EPL cell medium (90% DMEM, 10% FBS, Pen/Strep, 0.1 mM β-Me, 330 U/ml LIF, 400 μM L-proline; Sigma-Aldrich) for 6 days, with passage every 2 days. As controls, ESCs were also cultured for 6 days with 330 U/ml LIF or allowed to differentiate spontaneously without LIF or L-proline (Glover et al., 2022).

The effect of signalling pathway inhibitors (alone and in combination) on the transition of ESCs to EPL cells was assessed (Fig. 2D,E; Fig. S1D,E). The inhibitors were as follows: Mek1 inhibitor U0126 (U; 5 μM, Selleck); Fgfr inhibitor SU5402 (S; 5 μM, MedChem Express); Pi3k inhibitor LY294002 (L; 5 μM, Selleck); mTorc1 inhibitor rapamycin (R; 10 nM, Selleck) and S6k inhibitor PF-4708671 (P; 10 μM, MedChem Express). All inhibitors were dissolved in DMSO, and a vehicle control containing the maximum concentration (0.22%) of DMSO was included.

At days 2, 4 and 6, ESCs cultured in 1000 U/ml LIF or 330 U/ml LIF, and ESCs treated with L-proline±inhibitor(s) were analysed for three emergent properties (cell number, apoptosis and proliferation) and/or phosphorylation of various signalling pathway intermediates. At day 6, cells were also assessed for colony morphology, changes in gene expression and differentiation as embryoid bodies, as described below. Data were collected over five independent experiments.

Cell counts were measured with a haemocytometer following the addition of 0.4% Trypan Blue solution (Glover et al., 2022) to a single-cell suspension obtained following trypsinisation.

Colony morphology was quantified by images from an Olympus IX-81 inverted microscope. Images were deidentified and colony morphology scored based on a predetermined scale: round, domed (ESC) colonies were scored as 0; flat, irregular, partially transitioned colonies were scored as 1; and fully transitioned colonies consistent with EPL cells were scored as 2 (Glover et al., 2022). Scoring was performed on all colonies (10-40 per image) over three representative images from each condition. The sum of the score was divided by the total number of colonies scored, and then averaged across the three images to produce a final score.

After 6 days of culture in adherent culture, cells were passaged and 1.5×10⁶ were transferred to suspension culture plates and allowed to spontaneously differentiate without LIF or L-proline as EBs. EBs were collected at days 2, 3 and 4 and analysed by qRT-PCR for expression of the primitive streak marker T (primer sequences are provided in Table S3).

Total RNA was extracted from cells using GeneElute Mammalian Total RNA MiniPrep Kit (Sigma-Aldrich), including on-column DNase treatment to remove any contaminating DNA. RNA was converted to cDNA using a High-Capacity cDNA Reverse Transcriptase Kit (Applied Biosystems). qPCR was run on 10 μl reaction volumes containing 3 μl of 0.5 ng/μl cDNA, 2 μl of 1 μM primer (equal mix of forward and reverse primers; Table S3) and 5 μl of 2× SYBR Green master mix (Sigma-Aldrich) in a 384-well plate using a Roche LightCycler 480 with the following parameters: 15 min at 95°C, followed by 40 cycles of 30 s at 95°C, 60 s at 60°C, 30 s at 72°C. Thermal melt curves were obtained following this by increasing from 60°C to 95°C at 2.5°C/s. Threshold (C_t) values were used to calculate relative expression to the reference gene Actb (encoding β-actin), employing REST v9 software. Results were normalised to untreated ESCs and transformed to log₂ fold changes. All samples were tested to ensure that the C_t values for Actb were consistent (20±1 SD).

Cell samples were washed in ice-cold PBS and lysed (1 μl lysis buffer per 4×10⁴ cells) in the presence of protease and phosphatase inhibitors (Table S4). For data in Fig. 2C, cells were serum starved in 90% DMEM, 0.1% FBS and 0.1 mM β-Me for 4 h before sample collection. Cell lysates were incubated on ice for 10 min and then centrifuged at 4°C at 12,000 rpm (20,000 g). The supernatant was loaded onto a 1.5 mm 12% polyacrylamide gel with a 4% stacking gel. Molecular weight markers (Bio-Rad Precision Plus Protein Standards) were also loaded. Electrophoresis was carried out in a Bio-Rad western blot chamber at 100 V for 2 h.

Following electrophoresis, proteins were transferred to a 0.45 μm nitrocellulose membrane (Bio-Rad) for 120 min at 100 V using a Bio-Rad transfer system. The membrane was blocked overnight in Odyssey Blocking buffer (LiCor) at 4°C, then washed 3×5 min in Tris buffered saline with Tween 20 (TBST) and then incubated with primary anti-phosphoprotein antibody overnight at 4°C with rocking. Anti-β-tubulin antibody was used to stain for the reference protein. The membranes were then washed 3×5 min in TBST before 2 h incubation at room temperature in the dark with fluorescently labelled secondary antibody. Primary and secondary antibodies were diluted in Odyssey Blocking buffer with 0.1% (v/v) Tween 20. For details of antibodies and dilutions, see Table S5.

Membranes were imaged using an Odyssey Infrared Imaging system (LiCor). Bands were identified based on the protein size, and the integrated intensity was calculated by Image Studio (LiCor). Selected bands for each antibody are shown in Fig. S1F. Background signal was accounted for by subtracting the integrated intensity of a size matched adjacent region (Fig. S1F). Data were normalised to β-tubulin to correct for differences in loading, and then to untreated ESCs.

Flow cytometry was performed on a FACS Calibur and the results quantified using FlowJo software. Apoptosis was assayed using detection of annexin V. Live cells were centrifuged at 1200 rpm (2400 g) for 2 min, washed in PBS and recentrifuged, and then resuspended in 100 μl annexin V binding buffer with either FITC-annexin V conjugated antibody (1:33 dilution in TBST) and Propidium Iodide staining solution (BD Pharmingen) or PE-annexin V conjugated antibody (1:33 dilution in TBST) and 7-AAD as per kit instructions (BD Pharmingen). Samples were analysed by flow cytometry within 30 min.

Proliferation was assayed using BrdU incorporation and processed using the FITC BrdU Flow Kit (BD Pharmingen). Briefly, BrdU was added to cells in culture at a final concentration of 10 μM and incubated for 1 h. Cells were passaged, washed in PBS, fixed in BD Cytofix/Cytoperm and stored at −80°C in BrdU freezing buffer until required. Thawed samples were then stained according to the manufacturer's instructions before flow cytometry.

We employed several methodologies to explore relationships within our data (correlation matrices and unsupervised learning) and to understand which signalling pathways contribute to changes in gene expression and emergent properties (supervised learning). Additional statistical analysis was performed, including one-way ANOVA with Tukey's or Dunnett's multiple comparisons test or a two-tailed t-test. Details for each test are included in the figure legends. All analyses were performed in R.

The correlation matrices were generated from paired data using Pearson correlations with significance levels based on rank correlation (P<0.05) (Hoyt et al., 2008). Matrices were calculated using the Hmisc R package, and parameters were ordered based on hierarchical clustering.

PCA was generated using the paired data matrix. Data were z-scored, and the silhouette score was calculated using the factoextra package. Kmeans clustering was performed using the stats package, using two centres, as indicated by the silhouette score, and 25 iterations. PCA was generated using the factoextra package, and samples were coloured by cluster (Fig. S7A). Boxplots were generated from the original paired data frame using the kmeans cluster assignments.

Gene expression and emergent properties (cell number, proliferation, apoptosis and morphology) were modelled with (1) standard MLR (Vittinghoff et al., 2012); (2) MLR with two-way interaction terms (Flanders et al., 1992); or (3) BRANNGP (Burden and Winkler, 2008; Winkler and Burden, 2012).

Before modelling, each inhibitor was encoded using a 1-hot descriptor (1 when present, 0 when absent). For modelling cell number, proliferation and apoptosis, 1-hot variables were also used to represent each experimental day (either 2, 4 or 6). Each condition with an average of three replicates was used as input for modelling, with replicates averaged before modelling. As conditions containing both LY294002 and rapamycin resulted in cells not being viable after day 2, these were excluded from modelling on days 4 and 6. Data were subject to Shapiro-Wilks test, and morphology, apoptosis and proliferation data were transformed to improve normality. To ensure linearity, the residuals for each model were also measured for normality using a Shapiro-Wilks test. Model fitting parameters, including adjusted R² and σ values can be found in Table S1. Adjusted R² was used for comparability across modelling styles. Models with a higher adjusted R² and lower σ values are considered to have better fit. F-tests were also calculated to compare MLR with MLR with interaction effects (Table S2), where a P-value<0.05 indicates that the more complex models significantly improve the explanatory power of the model.

As this data contained all permutations and no predictive capacity was required, all the data were used to train models. To assess the range of responses from splitting the data, we generated 50 random 80% training/20% test models and profiled the range of responses seen in Table S6.

Models were generated using both the MLR and BRANNGP (sparse three-layer feedforward Bayesian regularised neural network with a Gaussian prior; Burden and Winkler, 2008). These were implemented in the in house software package Biomodeller. The latter method automatically optimises the complexity of the model (number of weights) to maximise predictivity. Models were trained until they reached the maximum of the evidence for the model so no validation set was required to provide a stopping criterion, which is important given the small dataset sizes. These models employed two neurons in the hidden layer, linear transfer functions in the input and output layer neurons and sigmoidal transfer functions in the hidden layer neurons. Data applied to the input layer was column scaled. See Burden and Winkler (1999) for a detailed explanation of BRANNGP methodology.

The authors acknowledge support, training and help from Bosch Core Facilities at The University of Sydney and, in particular, Drs Donna Lai and Sheng Hua (Bosch Molecular Biology Core Facility), as well as Drs Angeles Sanchez-Perez and Shirley Nakhla (Bosch Live Cell Analysis Facility).

The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.201704.reviewer-comments.pdf.