Combined transcriptomic and phosphoproteomic analysis of BMP4 signaling in human embryonic stem cells

combination of transcriptomic and phosphoproteomic approaches. Results indicate that in a pluripotent state, hESCs maintain WNT signaling under negative regulation by expressing pathway inhibitors. Initial stages of differentiation are characterized by upregulation of WNT pathway ligands, TGF β pathway in- hibitors which have been shown in Xenopus to expand the BMP signaling range essential for embryonic patterning, and mesendodermal transcripts. Moreover, BMP4 enhances the phosphorylation of proteins associ- ated with migration and transcriptional regulation. Results further indicate the vital regulatory role of Activin A and BMP4 in crucial fate decisions in hESCs. was carried out to identify novel, differentially regulated transcripts as primary and secondary targets of the cytokine. Results indicate that in a pluripotent state, hESCs maintain WNT signaling under negative regulation by expressing pathway inhibitors. Initial stages of differentiation are characterized by upregulation of specific TGF β pathway inhibitors and mesendodermal transcripts. Interestingly, secondary events include activation of the WNT pathway and suppression of neuroectodermal markers. Furthermore, stably labelled hESCs were induced with BMP4 and mass spectrometry was utilised to identify phosphoproteomic changes. Results suggest that hESCs preserve the phosphorylation status of WNT pathway and pluripotency associated at the early stages of differentiation. BMP4 stimulation increases the phosphorylation of proteins involved in signaling, migration and Our dataset contains several interesting differentially regulated phosphosites. Our findings reflect the changes undergone by as exit and acquire a fate. novel in this study further the ways by which a single stimulus orchestrates differentiation.


Introduction
Due to their pluripotent profile, entailing limitless in vitro selfrenewal potential and the ability to adopt nearly any cell identity, human embryonic stem cells (hESCs) have countless applications in regenerative medicine (Dupont et al., 2019). However, harnessing their therapeutic potential requires an extensive understanding of the elaborate signaling pathways implicated in the maintenance of pluripotency and control of differentiation. To understand the complexity of this network, hESCs need to be propagated in a fully characterised and versatile culture system. In that direction, hESCs are commonly cultured on matrigel-coated surfaces in mTeSR1. Matrigel is the extracellular matrix produced by mouse teratocarcinoma cells and mTeSR1 is a chemically defined medium relying primarily on TGFβ1 and FGF2 for pluripotency preservation (Kleinman, 2001;Ludwig et al., 2006). TGFβ1 leads to the phosphorylation of SMAD2/3 proteins which promote the upregulation of NANOG by directly binding to its proximal promoter (Xu et al., 2008). NANOG forms an autoregulatory loop with OCT4 and SOX2 to safeguard pluripotency (Vallier et al., 2009). In parallel, FGF2 leads to the activation of PI3K-AKT and MEK-ERK pathways that have opposing effects on pluripotency (Armstrong et al., 2006;Kang et al., 2005). Elevated PI3K signaling suppresses ERK and inhibits the pool of GSK3β involved in the destabilization of c-MYC while promoting cell proliferation (Singh et al., 2012). Low PI3K signaling leads to inhibition of GSK3β by ERK and activation of the WNT/β-catenin cascade (Singh et al., 2012). In this context, TGFβ-phosphorylated SMAD2/3 proteins cooperate with WNT components to drive differentiation towards endoderm (Chen et al., 2012).
Previous developmental studies in mice have highlighted the significance of TGFβ superfamily members on early embryonic patterning. In the pre-implantation blastocyst, NODAL preserves epiblast pluripotency by activating the SMAD2/3 cascade (Mesnard et al., 2006). Interestingly, the same factor has a prominent role during gastrulation, acting as a morphogen in concert with BMP and WNT. This process starts with BMP ligands originating from the extra-embryonic ectoderm activating WNT in both the epiblast and hypoblast (Rivera-Pérez et al., 2005). This leads to the induction of NODAL which acts by maintaining BMP signaling at the extra-embryonic ectoderm (Brennan et al., 2001). Recent experiments using hESCs on micropatterned surfaces demonstrated that BMP4 induction can recapitulate human gastrulation in vitro by sequentially activating WNT and NODAL (Deglincerti et al., 2016;Warmflash et al., 2014). In standard hESC culture conditions, activation of the SMAD2/3 cascade which is pivotal during both pluripotency and differentiation can be accomplished by TGFβ1 or Activin A (Beattie et al., 2005). Here, the outcome of BMP4 induction relies on the activation status of ERK. In the absence of ERK signaling, BMP4 leads to NANOG downregulation via the SMAD1/5/9 cascade and differentiation towards extraembryonic lineages. In the presence of FGF2, ERK signaling prolongs NANOG expression and switches the fate of differentiation to mesendoderm (Yu et al., 2011).
In the present study, we utilise a modified Activin A-based mTeSR variant (ActA-mTeSR) for the propagation of hESCs and employ a combined transcriptomic and phosphoproteomic approach to study the onset of BMP4-induced mesendoderm differentiation. We provide insights into the initial transcriptional changes occurring using RNA-seq, and identify novel phosphopeptides and map their location by mass spectrometry (MS).

Dynamic BMP4-induced transcriptomic changes in hESCs post 1 h and 4 h of treatment
To elucidate the transcriptional changes induced by BMP4, H1 cells cultured in ActA-mTeSR were treated with the cytokine for 1 h or 4 h and RNA-seq was carried out. Principal component analysis (PCA) was employed to assess the transcriptional relatedness of the samples. Each experimental condition within the RNA-seq dataset grouped into a defined cluster within the PCA plot and biological replicates within each experimental group also clustered together (Fig. 1A). Results show that 61 and 367 genes were significantly regulated after 1 h and 4 h of BMP4 induction, respectively. Taken together, these suggest that two snapshots of the transcriptomic profile of differentiating H1 cells were documented, originating from primary and secondary responses to BMP4 signaling. Further analysis revealed that 58 genes were upregulated and 3 downregulated after 1 h of induction (Fig. 1B). After 4 h, the numbers increased to 303 and 64, respectively (Fig. 1C). While investigating the kinetics of all differentially regulated genes, four separate clusters emerged (Fig. 1D). To decipher the biological pathways (BPs) they regulate, gene ontology (GO) analysis was performed. Results indicate that in a pluripotent state, cells maintain WNT signaling under negative regulation by expressing pathway inhibitors TLE4 and SHISA2. TLE proteins are known to recruit negative chromatin modifiers on canonical WNT transcription sites (Brantjes et al., 2001). In mice, SHISA2 has been found to affect WNT receptor FZD3 trafficking to the plasma membrane by inhibiting its glycosylation (Onishi et al., 2017). Both genes are progressively downregulated upon BMP4 addition. To the best of our knowledge, this is the first study reporting these WNT inhibitors being regulated by BMP4. In parallel, WNT activators WNT9B and WNT3 are upregulated as differentiation progresses. This is in agreement with studies showing WNT to be detrimental for pluripotency but essential for driving, along with BMP4, primitive streak formation (Sumi et al., 2008). Furthermore, hESCs are known to have a differentiation preference towards the neuroectodermal fate (Muñoz-Sanjuán et al., 2002). This is further validated by our data showing that genes involved in brain development are initially present but gradually downregulated as the cells acquire a mesodermal phenotype characterised by expression of the cardiac markers ACTC1, ANKRD1 and HEG1. In agreement, WNTinduced cardiomyogenesis has been well characterised in the literature, with BMP4 and Activin A participating in the onset and its progression (Kim et al., 2015). Furthermore, immediately after BMP4 addition several cell cycle regulators are transiently upregulated. Amongst them were HES1, a suppressor of p27 Kip1 (Murata et al., 2005), and Cyclin E2, a regulator of late G1 and early S phase (Caldon et al., 2010). This correlates with recent research suggesting that mesoderm specification is influenced by cell cycle effectors (Yiangou et al., 2019).

Transcriptomic landscape of hESCs post 1 h BMP4 treatment
Here, we carried out GO analysis considering all upregulated genes after 1 h of BMP4 treatment. Among the results, we see an enrichment of NOTCH and canonical WNT signaling coupled with negative regulation of SMAD-dependent pathways and an overall priming for cardiac mesoderm specification ( Fig. 2A). Inhibitory SMAD binding and regulation of TGFβ pathway components terms are also prevalent among the molecular functions (MFs) enriched (Fig. 2B). While NOTCH activation is in agreement with research considering it indispensable for differentiation to any germ layer (Yu et al., 2008), upregulation of TGFβ pathway inhibitors BAMBI, SMAD6 and SMAD7, in this context, is an interesting and novel finding. Research in Drosophila embryos demonstrated that BMPs collaborate with Sog, a BMP inhibitor, and Tld, a protease that cleaves Sog, to modulate development of the dorso-ventral axis (Eldar et al., 2002). In Xenopus, development of the same axis is selfregulated by BMP levels. Low ligand concentration leads to the expression of BMP agonist ADMP while increased levels upregulate the antagonist Crossveinless-2 (Ambrosio et al., 2008). BMP4-induced upregulation of BAMBI, and SMAD6/7 has been correlated with a dynamic output to different ligand doses, as expected from a morphogen. Knock-down of BAMBI in Xenopus embryos reduced the expression range of developmental genes ventx1/2, causing them to saturate in lower ligand concentrations (Paulsen et al., 2011). Overall, this could be a preserved mechanism of inhibitors being stimulated by and synexpressed with ligands in order to fine-tune the differentiation outcome. Furthermore, two interconnected KEGG pathways were significantly enriched in response to BMP4 (Fig. 2C). The TGFβ signaling pathway is finely modulated by inhibitory SMADs and BAMBI. In addition, inhibitors of DNA-binding (ID) are also enriched. IDs heterodimerize with basic Helix-Loop-Helix (bHLH) transcription factors to mitigate DNA binding and have also been reported to promote cardiac mesoderm in hESCs (Cunningham et al., 2017). As expected, several stem cell differentiation cues were present post BMP4 induction. KLF4 is known to modulate the expression of NANOG (Chan et al., 2009) and HAND1 is involved in mesoderm and cardiac development (Richter et al., 2014).

Transcriptomic landscape of hESCs post 4 h BMP4 treatment
GO analysis of the 303 upregulated transcripts after 4 h of BMP4 induction was carried out. Results indicated a further cardiac specification along with negative SMAD regulation and negative oligodendrocyte differentiation (Fig. 3A). Here, cardiac differentiation relies primarily on MSX2. The upregulation of MSX2 is in agreement with research showing BMP4 and WNT to synergistically activate it for mesendoderm and cardiac differentiation to progress (Wu et al., 2015). Negative regulation of SMAD signaling is further assisted by the enriched transcripts NOGGIN and PMEPA1. Recent studies on hESCs cultured on micro-patterned surfaces suggest that BMP4-induced NOGGIN confines phosphorylated SMAD1 to the colony edges (Etoc et al., 2016;Sorre et al., 2014). Our results verified NOGGIN upregulation in a different culture setup and pinpointed the exact timepoint. While research around PMEPA1 in hESCs is limited, cancer studies have shown it promotes epithelial to mesenchymal transition by activating BMP signaling (Zhang et al., 2019). In agreement with this, Cadherinmediated cell-cell adhesion is significantly reduced at this timepoint (Fig. 3B). Furthermore, PMEPA1 is also known to sequester R-SMADs and promote PI3K/AKT signaling (Singha et al., 2014). Therefore, elevated PI3K activity along with R-SMAD binding upregulation (Fig. 3C) could be attributed to PMEPA1 diverting TGFβ signaling to the non-canonical route.
The KEGG pathways significantly enriched are the same between the two BMP4 induction points (Fig. 3D). However, DUSP9 and FGFR3 which are members of the FGF/ERK cascade are only enriched after 4 h of BMP4 induction. In mouse ES cells, DUSP9 has been found to reduce ERK phosphorylation in response to BMP4 (Li et al., 2012). In hESCs, the same process is effected by DUSP6 which leads to prolonged NANOG expression during BMP4 induction that drives mesendoderm differentiation (Yu et al., 2011). A recent study reported high DUSP9 expression levels and low ERK phosphorylation in naïve hESCs. This trend was reversed as the cells entered a primed state (Jiapaer et al., 2018). Therefore, DUSP9 might be exerting similar effects in mouse and hESC.

Phosphoproteomic changes induced by BMP4 in hESCs
Studying the phosphoproteome of hESCs is a relatively recent endeavour. BMP4 addition combined with FGF2 removal has been reported to trigger changes in approximately half the phosphoproteome of hESCs (Van Hoof et al., 2009). To investigate the changes induced by BMP4, H1 cells cultured in ActA-mTeSR were metabolically labelled and treated with the cytokine for 1 h. Samples were injected in the MS and analysed with MaxQuant (Fig. 4A).
BMP4 induction led to the upregulation of 39 (Table S1) and downregulation of 174 phosphosites, respectively (Table S2). Network enrichment analysis on the upregulated phosphosites identified the most significantly regulated BPs which are depicted in (Fig. 4B). To gain insight in the signalling cascades functioning in early differentiation steps, we screened upregulated phosphosites for putative responsible kinases. Our results suggest that CDK1, CDK5, MAPK1, MAPK3, GSK3a, TTK, PKBa and PKC-B are potentially triggered by BMP4 (Table S1). CDK1/2 were previously reported to play a pivotal role both in pluripotency and differentiation of hESCs (Van Hoof et al., 2009). Furthermore, PKC is involved in extraembryonic endoderm specification in hESCs (Feng et al., 2012). CDK5 is a kinase with crucial roles in neurons and is associated with motility and survival (Roufayel et al., 2019). Although its role in hESCs has not been thoroughly characterised, CDK5 is required in EMT in cancer cells (Liang et al., 2013). In the same context, BMP4 inhibition causes TTK downregulation (Yan et al., 2012) with the latter being also involved in EMT (Chen et al., 2018). hESCs exposed to BMP4 undergo ЕМΤ at the onset of mesodermal differentiation (Richter et al., 2014), therefore future research could dissect the possibility that BMP4 induces EMT in hESCs via CDK5 and TTK.
14% of the upregulated phosphosites have been previously reported to be induced by BMP4 in hESCs (Van Hoof et al., 2009) while 37% were not identified. Within the latter group, identified phosphosites with potentially important implications in hESCs are PAK4 Ser 99 , BMP2K Ser 391 , Ser 392 , Ser 394 and PDH1 Ser 232 . PDK1 is known to phosphorylate PAK4 at Ser 99 and regulate its activity in migrating cells (Bastea et al., 2013). Therefore, our finding implies a possible crosstalk between BMP4 and PI3K signaling. Endogenous BMP2 signaling in hESCs is implicated in extra-embryonic endoderm differentiation (Pera et al., 2004). Consequently, we could hypothesize that BMP4 regulates fate acquisition through regulation of several TGFβ superfamily members. PDH converts pyruvate to acetyl-coA to fuel the TCA cycle and ATP production. PDH phosphorylation at Ser 232 displays an inhibitory effect on the protein (Rardin et al., 2009). Stem cells primarily use pyruvate for lactate production (Song et al., 2019). The shift from lactate production to the TCA cycle occurs in gastrulation (Conaghan et al., 1993). This finding indicates that metabolic reprogramming to favour the TCA cycle probably occurs at a later stage of the differentiation process.
Network enrichment analysis on the downregulated phosphosites identified the most significantly regulated BPs which are depicted in (Fig. 4C). In this group, NEDD4L phosphorylation on Ser 448 was revealed. NEDD4L is a ubiquitin ligase that mediates SMAD2/3 proteasomal degradation with Ser 448 phosphorylation abolishing SMAD2/3 binding (Gao et al., 2009). This suggests that BMP4 negatively regulates TGFβ signaling via increased SMAD2/3 proteasomal degradation. Studies in human gastruloids concluded that cell fate acquisition after BMP4 induction relies on a balance of SMAD1/SMAD2 activation that is regulated by endogenously expressed TGFβ and BMP inhibitors (Warmflash et al., 2014). Our results indicate an additional level of complexity in the BMP4 morphogenic potential based on the cytokine's ability to affect the phosphorylation status of key signaling regulators, such as NEDD4L.
Next, we mined our dataset for phosphorylation sites on proteins encoded by genes which were previously described to strongly correlate with NANOG (pluripotency associated genes) (International Stem Cell Initiative et al., 2007) and identified sites on 5 proteins out of 20 ( Fig. 4D upper panel). DNMT3B Ser 136 , Thr 383 , Lin28A Ser 200 and SOX2 Ser 251 are phosphorylations that have been previously identified (Van Hoof et al., 2009) while to the best of our knowledge, Lin28B Thr 202 , Ser 203 and GRB7 Ser 361 are novel, in the context of hESCs. Interestingly, BMP4 downregulated Lin28A phosphorylation at Ser 200 , whereas the rest of the sites were not regulated after 1hr of ligand induction. The actual role of these modifications remains to be elucidated but since LIN28B and GRB7 are involved in metabolism and migration (Morris et al., 2013;Shen et al., 2004), they are expected to regulate these processes.
Furthermore, given the known importance of WNT in hESCs we attempted to identify phosphorylation sites on proteins involved in this pathway (Fig. 4D lower panel). We found known phosphosites on CTNNB (Ser 191 , Ser 552 , Thr 556 , Ser 675 ) and on proteins that modulate WNT either positively or negatively. BMP4 upregulated CTNNB phosphorylation at Ser 191 while the rest of the sites were not regulated by the ligand. Ser 191 is associated with nuclear translocation (Wu et al., 2008) and possibly degradation protection (Muñoz et al., 2007), thus indirectly implying WNT pathway activation. CTNNB phosphorylation at Ser 552 is linked to protein translocation from cell junctions towards the cytoplasm and the nucleus (Fang et al., 2007) as well as with increased CTNNB-dependent transcriptional activity (Taurin et al., 2006). Sequential amino-terminal phosphorylations that are known to target CTNNB for proteasomal degradation were not present in our dataset (Liu et al., 2002(Liu et al., , 1999. Therefore, available evidence does not suffice to conclude about CTNNB status. Nevertheless, findings suggest that as expected the WNT pathway is under regulation in hESCs.

Conclusion
In the present study we utilise a modified mTeSR medium that enables the Activin A -dependent culture of hESCs to study the early BMP4induced changes in their transcriptome and phosphoproteome. Cells cultured in that system were induced with BMP4 and RNA-seq was carried out to identify novel, differentially regulated transcripts as primary and secondary targets of the cytokine. Results indicate that in a pluripotent state, hESCs maintain WNT signaling under negative regulation by expressing pathway inhibitors. Initial stages of differentiation are characterized by upregulation of specific TGFβ pathway inhibitors and mesendodermal transcripts. Interestingly, secondary events include activation of the WNT pathway and suppression of neuroectodermal markers. Furthermore, stably labelled hESCs were induced with BMP4 and mass spectrometry was utilised to identify phosphoproteomic changes. Results suggest that hESCs preserve the phosphorylation status of WNT pathway components and pluripotency associated proteins at the early stages of differentiation. BMP4 stimulation increases the phosphorylation of proteins involved in signaling, migration and metabolism. Our dataset contains several interesting differentially regulated phosphosites. Our findings reflect the changes undergone by hESCs as they prepare to exit pluripotency and acquire a mesodermal fate. The novel targets described in this study further broaden the ways by which a single stimulus orchestrates differentiation.
Funding This work was supported by funding from the School of Biosciences, University of Birmingham.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data
Supplementary data to this article can be found online at https://doi. org/10.1016/j.scr.2020.102133.