Mechanisms underlying WNT-mediated priming of human embryonic stem cells

Understanding the transition from pluripotency towards differentiation of human embryonic stem cells (hESCs) is a crucial step in elucidating specific aspects of human development. The activin/nodal signaling pathway, mediated by transcription factors SMAD2 and SMAD3 (SMAD2/3), plays a dual role in this process. On one hand, nuclear SMAD2/3 promote the expression of pluripotency factors OCT4 and NANOG and are essential for hESC pluripotency maintenance (Beattie et al., 2005; James et al., 2005; Vallier et al., 2005; Xiao et al., 2006); on the other hand, in combination with WNT/β-catenin signaling, SMAD2/3 can induce hESC differentiation into primitive streak derivatives (Loh et al., 2014; Singh et al., 2012; Sumi et al., 2008), including human organizer tissue (Martyn et al., 2018).

Classical work in vertebrate embryos, originally performed in Xenopus embryos and subsequently in fish and mouse, emphasized the importance of activin/nodal signaling pathway when acting as a morphogen to induce and pattern mesendodermal (ME) lineages, including the organizer/primitive streak (Dunn et al., 2004; Green and Smith, 1990; Wilson and Melton, 1994). The idea that activin acts as a morphogen has consequently been adopted by in vitro differentiation studies of hESCs: low levels are required to maintain pluripotency; moderate to high levels have been assumed to push the embryonic cells toward ME differentiation. Previous work has demonstrated that hESCs can maintain pluripotency over a wide range of activin concentrations and that WNT signaling is required for activin-mediated dose-dependent differentiation (Funa et al., 2015; Martyn et al., 2018; Singh et al., 2012; Yoney et al., 2018). The mechanism by which WNT alters the cellular response to different levels of activin/nodal signaling and, more generally, how these two signaling pathways synergize to induce ME differentiation is not fully understood.

Several mechanisms have been proposed to explain how WNT and activin signaling interact to drive ME differentiation in ESCs. For example, it was proposed that the WNT effector, β-catenin, physically interacts with SMAD2/3 to direct their binding and activation of target genes (Funa et al., 2015; Wang et al., 2017); alternatively, it was suggested that β-catenin may not directly interact with SMAD2/3, but that it activates transcription synergistically by promoting different transcriptional activation steps (Estarás et al., 2015); β-catenin binding may also help SMAD2 overcome a repressor complex, TAZ/YAP/TEADs, at ME-specific genes and thus enable their activation (Beyer et al., 2013; Estarás et al., 2017).

A common theme of these mechanistic models is that they are all based on simultaneous WNT and activin stimulation and thus involve co-binding of β-catenin and SMAD2/3 to the target ME genes. In contrast, our previous study revealed that the coordination of WNT and activin/nodal involves a signaling cascade of WNT acting upstream of activin: hESCs ‘primed’ with WNT for 24 h with subsequent exposure to activin can differentiate into ME cells, although each signaling pathway alone is unable to induce this differentiation (Yoney et al., 2018). Our finding that these two signaling pathways can be temporally separated suggests a new mechanism that mediates the synergy between WNT and TGFβ signaling in the context of ME differentiation; it also provides a new way to decouple and manipulate the differentiation process.

Using the WNT to activin stepwise stimulation protocol, we found that WNT priming does not affect the binding or chromatin opening by SMAD2/3, as suggested by previous models. A brief overlap between nuclear β-catenin and SMAD2/3 is also not sufficient for hESCs to exit pluripotency, and once hESCs have been primed, nuclear β-catenin is dispensable for ME differentiation. Instead, we find that differentiation crucially depends on the induction of the transcription factor EOMES during WNT priming, a key regulator of ME differentiation across vertebrates (Arnold et al., 2008; Brown et al., 2011; Bruce et al., 2003; Ryan et al., 1996). Exogenously expressed EOMES can completely replace WNT signaling and drive ME differentiation together with SMAD2/3. Therefore, this study reveals that the primary function of WNT signaling during ME differentiation is not to globally interact with SMAD2/3 but rather to induce EOMES expression.

We have previously shown that WNT and activin signaling can be presented in temporal succession to induce mesendoderm from pluripotent hESCs (Yoney et al., 2018). 24 h of WNT treatment followed by 24 h of activin stimulation leads to significant reduction in the expression of the pluripotency marker SOX2, and to strong induction of the primitive streak and ME marker genes TBXT/BRA and EOMES, and the anterior ME marker GSC (Fig. 1A-C and Fig. S1A,B). In a typical experiment, we observed a high number of cells with ME marker expression. From three independent replicates of the WNT/activin treatments we observed 94%, 99% and 88% of cells with at least one of the ME makers (BRA, EOMES or GSC) expressed at a level greater than three standard deviations above the average background level in pluripotency (–/activin). Although EOMES and GSC expression are highly correlated in ME cells, both markers show low correlation with BRA (Fig. S1C). These differences in correlation likely reflect the differences in the dynamics of gene expression, as BRA is induced in the primitive streak/mesendoderm but is turned off as cells further differentiate towards anterior mesendoderm and definitive endoderm (DE) (Gu et al., 2004; McLean et al., 2007; Teo et al., 2011).

The results following WNT/ACT contrasts with 24 h of WNT stimulation where the expression of OCT4 and SOX2 is maintained at pluripotency levels, GSC remains low, and EOMES and BRA are only mildly activated (Fig. 1B and Fig. S1A,B). NANOG expression is reduced, consistent with the removal of activin during WNT priming (Vallier et al., 2009). When cells were primed with WNT and then cultured in E7 base medium (Chen et al., 2011) without WNT or activin for an additional 24 h, EOMES and BRA transcript levels start to drop towards their basal levels in E7 (Fig. 1B). These observations suggest that 24 h of WNT signaling does not cause a stable fate change in hESCs. When activin was applied prior to WNT, ME gene expression remains low, strongly suggesting that ME differentiation does not simply require both signals in close succession but also in a particular temporal order (Fig. S1A,B).

As shown in Fig. S2, similar results were obtained with the small molecule CHIR99021 (CHIR), which activates the canonical WNT signaling pathway in a cell-autonomous manner by blocking β-catenin degradation (Kunisada et al., 2012; Naujok et al., 2014). Taken together, these results demonstrate that WNT priming through β-catenin is required to switch the output of activin signaling from pluripotency maintenance to ME specification.

The experiments presented above confirm our original observation that hESCs have a differential transcriptional response to activin alone, versus activin stimulation following WNT priming. To begin to decipher the mechanisms underlying this difference, we performed bulk RNA-sequencing (RNA-seq) under a set of WNT and activin conditions (Fig. 1A). We identified 1212 genes that showed a significant change in expression in at least one of the conditions relative to 24 h in E7 medium (n=2 biological replicates, adjusted P-value <0.01). K-means clustering of these genes revealed four major clusters. Cluster 1 includes genes for which the expressions show comparable increases in −/activin and WNT/activin, but less in WNT and WNT/−, indicating that these genes are induced by activin (n=349). Cluster 2 represents genes induced by WNT (n=206). Cluster 3 includes genes that show higher induction in WNT/activin than WNT or activin alone (n=222). Cluster 4 represents genes that are repressed by the WNT and/or activin condition (n=435) (Fig. 1C and Table S1).

Among all differentially expressed genes, 189 had significant changes (increased or decreased) after 24 h of WNT treatment (Table S1). Genes that were induced by WNT, including EOMES, BRA and NODAL, partially reversed to their basal level after a subsequent 24 h in E7 medium, which further supports the notion that WNT signaling primes cells rather than causing a stable fate change. The expression of pluripotency makers NANOG, OCT4, SOX2 and KLF4, either remained constant or was upregulated in response to WNT or activin (Fig. S1D), confirming the lack of differentiation with either signal alone. Consistent with broad transcriptional changes associated with differentiation, there were more extensive changes when WNT was followed by activin (WNT/activin), including 661 upregulated genes and 383 downregulated genes. In our WNT/activin condition, cells express markers of anterior mesendoderm and early DE, including EOMES, GSC, SOX17, FOXA2 and GATA4 (Fig. 1C and Table S1). We did not observe expression of key markers of cardiac, paraxial or lateral plate mesoderm, including MESP1/2, TBX6, MSGN1/2, HAND1 or ISL1. These results are consistent with previously published in vivo and in vitro studies, which demonstrated that cells differentiate towards the DE fate via a ME intermediate (Gu et al., 2004; McLean et al., 2007; Teo et al., 2011).

In our previous study we showed that the application of activin to hESCs induces transient translocation of SMAD2 and SMAD4 (SMAD2/4) into the nucleus, which elicits a transcriptional change in ∼3500 genes (Yoney et al., 2018). Many of the differentially expressed genes are activated rapidly following SMAD2/4 translocation and are repressed after its nuclear exit, suggesting that they are direct targets of SMAD2/4. Among these genes, 113 of them show significantly higher mRNA levels in WNT/activin compared to WNT- or activin-only conditions, indicating that these SMAD2 target genes respond to the cooperative signaling of WNT and activin (‘sensitive to WNT priming’, Table S2). This group of genes includes EOMES, BRA, GSC, MIXL1, GATA6 and NODAL, which agrees with previous reports that β-catenin and SMAD2/3 co-regulate expression of genes that are expressed within the primitive streak and ME derivatives (Estarás et al., 2015; Funa et al., 2015). We also collected SMAD2/4 target genes identified in our previous study that are not sensitive to WNT priming (Table S2).

In summary, our genome-wide transcriptional analyses support the idea that 24 h of WNT stimulation primes pluripotent hESCs for differentiation and that, together, WNT and activin drive cells towards a committed ME fate. We were able to unbiasedly identify a group of 113 activin/SMAD2 target genes for which the transcriptional responses were affected by WNT priming. In our subsequent analyses, we will compare WNT priming-dependent versus -independent genes with respect to additional features, including transcription factor (TF)-binding sites and chromatin accessibility, to further dissect the mechanism of WNT priming.

Motivated by the observations of ‘Wnt memory’ in other species ranging from fly to frog (Alexandre et al., 2013; Blythe et al., 2010), we further explored our original hypothesis that hESCs possess the ability to record WNT signals. To demonstrate extended WNT memory, we stimulated cells with WNT for 24 h, followed by one additional day in E7 medium before switching to activin for 24 h. The small molecule inhibitor SB431542 (SB) was also included to block the ability of cells to respond to endogenous activin/nodal signals and upregulate ME markers. Consistent with previous observations (Yoney et al., 2018), adding SB during Wnt priming had little effect on the subsequent activin response (Fig. 1D,E). Despite the extended interval between WNT and activin stimulation, Activin could still induce ME markers and downregulate SOX2 at the protein level, albeit to a lesser extent (Fig. 1D,E and Fig. S3). We interpret these results to mean that: (1) the WNT effect decays over time and the cells are not stably committed to the ‘WNT-primed’ state; and (2) the WNT priming effect is still available for at least 1 day after removal of the signal.

We next investigated the mechanism by which WNT priming modifies the cellular response to activin. We had previously shown that SMAD2 concentration and nuclear residency dynamics are not affected by WNT priming, suggesting that WNT affects SMAD2 activity in a subsequent step (Yoney et al., 2018). As SMADs have weak DNA specificity and often rely on co-factors to be directed to target genomic regions in a cell type-specific manner (Massagué, 2012), we first considered the possibility that WNT priming alters the target sites of SMAD2/4. Accordingly, we analyzed the SMAD2/3 ChIP-seq data of Kim et al. (2011) (ChIP in hESCs in conditioned media, which is similar to our –/activin condition) and of Tsankov et al. (2015) (ChIP in mesendoderm cells generated with WNT+activin treatment for 12h, which is similar to our WNT/activin condition). The overall SMAD2/3 ChIP profiles are similar in these two conditions. Several examples are shown in Fig. 2A, where the SMAD2 ChIP-seq signals are largely overlapping, despite the mRNA levels of these genes being a few fold or even orders of magnitude higher in the WNT priming followed by the activin condition (WNT/activin) than in the activin alone condition (–/activin) (Fig. S4A). Heatmaps confirm that SMAD2 binding sites detected in ESCs tend to be occupied by SMAD2 in ME cells, and vice versa (Fig. 2B). Tsankov et al. has one measurement of SMAD2 ChIP-seq in ESCs, which shows low signal-to-noise (Fig. S4B). Nevertheless, this dataset also shows weak peaks at annotated SMAD2-binding sites, both in ESCs and ME cells (Fig. S4B,C). We also quantified the area underneath the SMAD2/3 ChIP-seq peaks from these datasets. They are highly correlated, both at the genome-wide scale (R ∼0.8; Fig. 2C) or in the vicinity of genes that are sensitive to WNT priming (R=0.86; Fig. 2D; Table S2), whereas the average correlation between the three replica of the ME ChIP data is only ∼0.71. We therefore conclude that SMAD2 binding measured in these two datasets/conditions is highly similar.

As the published SMAD2/3 ChIP data were performed in cell lines and under conditions that are not identical to ours, we selected a few SMAD2-binding regions near genes that respond to WNT priming (GSC, CER1 and NODAL) and carried out ChIP-qPCR measurements at two different time points during activin treatment with or without WNT priming (Fig. 2E). For all genes tested, we observed an increase in SMAD2 binding at 2 h relative to 24 h after activin presentation, consistent with the adaptive change of the SMAD2 nuclear concentration over time (Yoney et al., 2018). However, the level of binding is comparable in the presence or absence of WNT priming. This analysis indicates that higher SMAD2 nuclear concentration promotes SMAD2 binding at target genes. However, SMAD2 binding per se is not correlated with transcriptional output and corresponding changes in cell fate.

Given that SMAD2 can bind to the same genes with or without WNT priming but activates them to a different extent, it is possible that SMAD2 may have different activation potentials with or without WNT priming. Because TF activation is often accompanied by local chromatin opening, we carried out ATAC-seq to probe the changes in chromatin accessibility. The ATAC-seq experiments were performed under the same conditions already described for our bulk RNA-seq analysis. The ATAC-seq profiles in all conditions are very similar, showing that there is no global rearrangement of open chromatin regions in the transition from pluripotency to ME differentiation (example traces are provided in Fig. 3A and global analysis in Fig. S5A,B). Furthermore, a vast majority of genes contain ATAC-seq peaks near their TSS and within transcript regions, regardless of their transcriptional status (see GEO deposit with accession number GSE176222).

Despite the global similarities, there are quantitative changes in the ATAC-seq signals over specific TF-binding sites. Two examples of ATAC-seq signal change near MIXL1, LEFTY1 and LEFTY2 are shown in Fig. 3A (peaks marked with asterisks). These dynamic ATAC-seq peaks overlap with the binding sites of several TFs, including β-catenin, SMAD2 and EOMES (Fig. S5C). For β-catenin, the ATAC-seq data clearly show that, globally, its binding sites are accessible in E7 medium before WNT stimulation (Fig. 3B). After 24 h of WNT treatment, β-catenin associates with these pre-opened chromatin regions and leads to mild but statistically significant enhancement in ATAC-seq signals (Fig. 3C). The level of enhancement correlates with the strength of β-catenin binding determined from the ChIP-seq peak area (Estarás et al., 2015). Interestingly, such enhancement is reduced when cells are moved back to E7 medium after WNT priming (WNT/–), indicating that β-catenin at least partially dissociates from chromatin after elimination of WNT signals. In addition, the ATAC-seq peaks drop to the same level when cells are primed with WNT and then stimulated with activin (WNT/activin) as in the WNT/– condition, suggesting that the presence of nuclear SMAD2 does not retain β-catenin on chromatin (Fig. 3C).

We performed the same analysis using SMAD2-binding sites (Kim et al., 2011; Tsankov et al., 2015). The ATAC-seq signals over the SMAD2 binding sites are essentially identical in E7 base medium and WNT conditions, indicating that WNT signaling does not alter the chromatin environment near these sites. In contrast, with both activin only (–/activin) and WNT priming followed by activin treatment (WNT/activin ), ATAC-seq signals are significantly enhanced over SMAD2-binding sites, and the level of enhancement correlates with SMAD2-binding strength (Fig. 3D,E and Fig. S5D). Therefore, such enhancement likely results from SMAD2 binding and subsequent opening of the local chromatin. Importantly, we observed similar enhancement in –/activin and WNT/activin samples over the SMAD2-binding sites, which supports our conclusion that SMAD2/3 binding is similar in these two conditions. Moreover, these data indicate that the activity of chromatin opening by SMAD2 is not dependent on WNT priming.

We have shown that WNT priming affects neither SMAD2 binding nor its ability to open chromatin. It has been proposed that the main TF that responds to WNT signaling, β-catenin, functions as a co-activator of SMAD2, either through direct physical interaction (Funa et al., 2015) or by cooperatively recruiting RNA pol II (Estarás et al., 2015). β-Catenin indeed tends to bind WNT-primed genes together with SMAD2 (Fig. S6A,B). It is thus possible that during activin treatment after WNT priming, the remaining β-catenin nuclear fraction synergizes with newly translocated SMAD2 to drive the expression of key ME genes (Fig. 4A).

Two pieces of evidence argue against this scenario as being the only requirement for ME gene induction. First, β-catenin rapidly accumulates in the nucleus after Wnt presentation, reaching the peak concentration at ∼6h in hESCs (Massey et al., 2019). If the overlap between nuclear β-catenin and SMAD2/4 is the key mechanism of WNT priming, a shorter period of WNT treatment should be sufficient for ME differentiation (Fig. 4B). However, 6 h or less of WNT priming followed by WNT washout and 24 h of activin in the presence of an in inhibitor of endogenous WNT secretion (IWP2), is not able to activate ME genes (Fig. 4C,D). The fact that WNT priming requires longer time strongly indicates that a product resulting from β-catenin activation, instead of β-catenin itself, is essential to turn on the ME differentiation program with SMAD2

Second, we tested whether the overlap between β-catenin and SMAD2 is necessary for ME differentiation. We repeated our WNT/activin protocol adding the small molecule inhibitor endo-IWR1 at varying times from 0 h to 12 h after the addition of activin (Fig. 5A). Endo-IWR1 promotes the degradation of β-catenin through stabilization of the axin destruction complex (Chen et al., 2009) and therefore should reduce its overlap with SMAD2 (Fig. 5A). Application of this drug together with WNT completely abolishes ME and DE differentiation, and maintains pluripotency gene expression (Fig. 5B-D). When endo-IWR1 was added immediately after WNT treatment (0 h) together with activin, the induction of some ME genes was reduced, indicating that a short overlap between β-catenin and activin signaling indeed is important for ME differentiation. However, when the drug was added a few hours (4-12 h) after the start of the activin treatment, it no longer blocked ME induction or the downregulation of SOX2 and NANOG (Fig. 5B,D). Importantly, marker genes for DE, including FOXA2, SOX17 and GATA4, were also induced with the later application of endo-IWR1, indicating that hESCs can differentiate into more mature endoderm after the elimination of β-catenin (Fig. 5C). We note that the expression of all genes tested is stable between 24 and 48 h of activin treatment, except for BRA, which is downregulated at 48 versus 24 h (Fig. 5B). This is consistent with previous findings that BRA is turned off as cells differentiate from a ME intermediate to more mature DE (D'Amour et al., 2005; McLean et al., 2007). Taken together with experiments presented in Fig. 4, we conclude that the presence of β-catenin in the nucleus is not sufficient to carry out the ME to DE differentiation program and, after 24 h WNT priming, continuous overlap between β-catenin and SMAD2 is also not necessary for ME and early DE differentiation.

Based on the results presented above, we hypothesized that WNT priming generates another TF that functions as an essential SMAD2 co-regulator (Fig. 6A). Given our data in Figs 4 and 5, we suspect that this TF should be activated by the end of the 24 h of WNT treatment and should bind in proximity to WNT-primed genes in a WNT/activin specific manner. Accordingly, we probed the binding signature of this TF by searching for enhanced ATAC-seq peaks that require the combined action of WNT and activin, e.g. the ATAC-seq peak near the MIXL1 gene in Fig. 3A. We identified 3983 peaks that are significantly enhanced in the Wnt priming followed by activin condition (WNT/activin) in comparison with both activin alone (–/activin) and Wnt alone (WNT/–) (Fig. 6B). The probability for enhanced peaks to fall within 20 kbp of WNT-primed genes is much higher than for genes that are not primed by WNT, indicating that these ATAC-seq peaks correlate with transcriptional activation and are therefore likely to represent enhancers (Fig. 6C).

We analyzed motifs that are enriched in the WNT/activin-enhanced ATAC peaks, both genome-wide and the subset that is proximal to WNT-primed genes. In both cases, we found the EOMES motif to be the most significantly enriched motif (the EOMES motif was detected in ∼80% of the sequences) (Fig. 6D; full list is shown in Table S3). Indeed, ∼80% of WNT/ACT-enhanced ATAC-seq peaks overlap with EOMES ChIP peaks (Tsankov et al., 2015, mesendoderm condition), whereas the overlap with β-catenin and SMAD2 peaks is much less (Fig. S6C). Consistent with our expectation, EOMES is activated by WNT alone and becomes highly expressed with WNT priming followed by activin (Fig. 1B and Table S1). EOMES tends to bind near SMAD2 in the vicinity of Wnt-primed genes (Fig. S6D,E), and is an essential TF for differentiation into ME cells and DE (Li et al., 2019; Teo et al., 2011; Tosic et al., 2019). The TF with the second most enriched motif, BRA (Fig. 6D), is not essential for DE differentiation in vitro (Li et al., 2019; Tosic et al., 2019). Although BRA-deficient mice die by embryonic day 10, they do not exhibit substantial defects in the anterior mesendoderm and definitive endoderm but fail to form midline and posterior mesoderm (Wilkinson et al., 1990). The other significantly enriched TF-family motif belongs to the GATA factors. GATA4/6 are key endoderm TFs, but their expression is not upregulated after 24 h WNT (Table S1). Therefore, we conclude that EOMES is the most likely effector of WNT priming.

If EOMES is the key factor mediating WNT priming, the prediction is that we can achieve ME differentiation in the absence of WNT by artificially expressing EOMES. To test this hypothesis, we generated an hESC line with doxycycline (DOX)-inducible EOMES (TT-EOMES,Fig. S7A) using the piggyBac transposon system (Lacoste et al., 2009). When these cells were treated with DOX, morphological changes were induced that were indicative of ME differentiation (Fig. S7B). According to previous studies, EOMES expression in the absence of activin drives cells towards cardiac mesoderm fate by activating MESP1, and addition of activin inhibits cardiac differentiation (van den Ameele et al., 2012). We observed the same trend of MESP1 expression in our dox-treated TT-EOMES lines in the presence or absence of activin (Fig. S7C). To prevent cells from differentiating towards this alternative path, we added activin when inducing EOMES with DOX in the protocols below.

We next tested whether a pulse of EOMES expression could replace WNT priming to drive ME differentiation. A titration experiment showed that exogenous EOMES is induced in a DOX concentration- and time-dependent manner (Fig. S8A). In particular, a 6 h treatment with 0.1 µg/ml DOX induces EOMES protein expression to a level comparable with the endogenous levels in our WNT/activin protocol (Fig. S8B,C). We induced TT-EOMES for 6 h with 0.05, 0.1 and 0.2 µg/ml dox in the presence of activin, washed out DOX and measured gene expression at different time points after DOX removal (Fig. 7A). TT-EOMES expression level decreases rapidly following the removal of DOX (Fig. 7B). Such transient TT-EOMES expression induced by 0.1 and 0.2 µg/ml, but not 0.05 µg/ml DOX, was sufficient to drive expression of ME and DE genes, including the endogenous EOMES. Except for BRA, the levels induced in the TT-EOMES line matched those obtained with WNT/activin treatment of the unmodified parental cell line (Fig. 7C). We noted that ME genes, such as EOMES and GSC, are activated earlier than bona fide DE genes, such as SOX17, indicating that the cells driven by TT-EOMES go through a differentiation process similar to that obtained with other in vitro DE protocols (McLean et al., 2007).

Consistent with a previous study (Pfeiffer et al., 2018), we also found that EOMES can induce the expression of WNT3 (Fig. S8D). To eliminate the possibility that TT-EOMES drives hESC differentiation through the production of endogenous WNT, we repeated our experiments in the presence of endo-IWR1. When applied to the unmodified hESCs in the WNT/ACT protocol, endo-IWR1 completely blocks the activation of ME and DE genes (Fig. 5B,C). In contrast, endo-IWR1 has little effect on the ME and DE gene expression driven by TT-EOMES (Fig. 7D). These results strongly suggest that EOMES is an effector of WNT priming, and exogenously expressed EOMES can drive ME differentiation in the absence of WNT signaling.

EOMES is a TF expressed during early development that is crucial for endoderm and cardiac mesoderm formation in mouse and human, as EOMES knockout completely abolishes these lineages (Arnold et al., 2008; Costello et al., 2011; Kartikasari et al., 2013; Li et al., 2019; Pfeiffer et al., 2018; Teo et al., 2011). It is also thought that EOMES is at the top of the ME gene regulatory network hierarchy, because it is one of the first genes that is induced during ME differentiation (Teo et al., 2011) and exogenous expression of EOMES in the presence of activin promotes the expression of endodermal genes (van den Ameele et al., 2012). However, because EOMES engages in a complex gene network with WNT and activin, e.g. EOMES can be activated with WNT signaling, and WNT3 is also a direct target of EOMES (Pfeiffer et al., 2018), it was not clear whether EOMES could completely bypass WNT and directly induce ME differentiation with activin. In fact, it was proposed that the activation of WNT signaling downstream of EOMES is crucial for ME and cardiac differentiation (Pfeiffer et al., 2018). Here, we show that induction of EOMES at near physiological levels in the presence of activin is sufficient for ME and DE induction without the involvement of β-catenin signaling. Although there is previous literature showing that EOMES is a co-factor for SMAD2 in ME differentiation, we specifically found EOMES and SMAD2 co-binding enriched in genes sensitive to WNT priming (Fig. S6E), which establishes a critical link between WNT and activin signaling in ME differentiation and provides a potential mechanism for WNT priming and memory (Yoney et al., 2018).

EOMES binds to its own promoter along with SMAD2; together EOMES and SMAD2 drive ME differentiation (Kartikasari et al., 2013; Teo et al., 2011). One question therefore raised by our study is why this positive-feedback loop is not stably activated with WNT priming alone, i.e. why does EOMES expression decrease when cells are returned to neutral media after WNT priming? In our protocol, EOMES is activated by 24 h WNT (mRNA level is increased by ∼200-fold in comparison with E7), but activin can induce its expression further (∼1500-fold in comparison with E7) (Fig. 1B). It is possible that EOMES is already associated with its own promoter in WNT but cannot drive itself to a high level without SMAD2. Alternatively, EOMES concentration may need to reach a critical level to jumpstart the feedback, and this level is not reached during the WNT priming phase. Further studies are needed to distinguish among these possibilities. Further investigation is also required to determine the mechanism by which EOMES and activin and/or SMAD2 lock in the ME fate. One possibility is that EOMES suppresses pluripotency factors or alternative differentiation pathways (Teo et al., 2011; Tosic et al., 2019). Our proposed mechanism for WNT priming is summarized in Fig. 8.

Chromatin modifications have long been implicated in cellular memory as the term epigenetics implies. In Xenopus, a delay between WNT signaling and transcriptional onset was linked to an epigenetic mechanism in which β-catenin recruits a histone methyltransferase to target genes and is required for their activation at the mid-blastula transition (Blythe et al., 2010). The competence of the Xenopus blastula to respond to WNT signals was also shown to have a chromatin component (Esmaeili et al., 2020). Closer to our study, Tosic et al. showed that EOMES and BRA are required for ME gene expression in mouse, and that their binding is associated with open chromatin. We did not find evidence for drastic changes in chromatin opening with WNT priming. β-Catenin-binding sites are already highly accessible in hESCs, which is likely because these sites are occupied by TCF/LEF proteins (Estarás et al., 2015), and β-catenin binding only mildly enhances the local ATAC-seq signals. In addition, the presence of β-catenin does not alter the accessibility near the SMAD2/3 sites, consistent with the idea that SMAD2/3 binding is not affected by WNT priming. Likewise, SMAD2/3 in the activin-alone condition does not enhance chromatin accessibility over the β-catenin-binding sites. Taken together, these findings suggest that, in contrast to what would be predicted from earlier models on ME differentiation, the chromatin opening and the direct binding of β-catenin and SMAD2/3 are independent events along the path to ME differentiation. Additionally, by superimposing our ATAC-seq analysis with our analysis of prior published datasets, we demonstrate a consistency among multiple similar protocols for mesendoderm induction.

Gene regulation by TGFβ signaling in different cell types often involves alternative SMAD2/3 binding partners. For example, SMAD2 binding in hESCs was proposed to depend on OCT4, NANOG and FOXH1 (Attisano et al., 2001; Brown et al., 2011; Chen et al., 1996; Mullen et al., 2011; Teo et al., 2011), and on EOMES and GATA6 in definitive endoderm (Li et al., 2019). However, we find that SMAD2 binding does not seem to differ in pluripotent hESCs versus ME cells. The similarity of SMAD2 binding may be due to the fact that SMAD2/3 co-binders, OCT4, NANOG and FOXH1, all have similar levels in these two cell types. This is different from the situation in definitive endoderm cells where OCT4 and NANOG are significantly downregulated, which likely allows SMAD2 to target other loci. Therefore, we propose that the induction and subsequent binding of EOMES provides transcriptional activation downstream of SMAD2/3 binding, potentially by EOMES-dependent recruitment of chromatin modifiers (Kartikasari et al., 2013). One potential caveat of our analysis is that we compared SMAD2-binding data from two studies/labs (ESC data: Kim et al., 2011; mesendoderm data: Tsankov et al., 2015). However, we believe this comparison is justified by the fact that the agreement between these datasets is better than the agreement among the three replicates of the ME condition generated by a single study and/or lab (Tsankov et al., 2015).

In summary, our findings reveal that one essential function of WNT signaling is the induction of EOMES, providing a simple mechanism of signaling integration that drives cellular differentiation. How induction of EOMES may impart a dose-dependent response to Activin/SMAD2 signaling that is responsible for further patterning the mesendoderm remains an open question.

Experiments were performed with the RUES2 hESC line (XX female; US National Institutes of Health, human ESC registry number 0013), which was authenticated by STR profiling. Cells were routinely tested for mycoplasma contamination. For maintenance, hESCs were grown in HUESM medium that was conditioned by mouse embryonic fibroblasts and supplemented with 20 ng/ml bFGF (MEF-CM). Cells were grown on tissue culture dishes coated with Geltrex (Thermo Fisher Scientific) and passaged as aggregates using Gentle Cell Dissociation Reagent (STEMCELL Technologies). For experiments, single cells were obtained by dissociation with Accutase (STEMCELL Technologies) and seeded at low density in optical-quality plastic tissue culture dishes (50 K per dish) or 24-well plates (5 K per well) (Ibidi) in TeSR-E7 medium (STEMCELL Technologies) supplemented with 10 µM Rock inhibitor (Y-27632, Abcam). Prior to seeding, dishes or plates were coated with 10 µg/ml laminin-521 (BioLamina) in PBS +Ca/+Mg for 2 h at 37°C or overnight at 4°C. For the 2-day protocol in which cells were switched from Wnt3a to activin A (R&D Systems), the samples were washed with PBS +/+ before adding fresh medium. 10 µM Rock inhibitor was maintained throughout the duration of the experiment. Other small molecules were used at the following concentrations and were replaced every 24 h: 2.5 µM CHIR99021 (EMD Millipore), 1 µM IWP-2 (Stemgent), 1 µM endo-IWR-1 (Tocris) and 10 µM SB431542 (Stemgent).

Cells were rinsed once with PBS, fixed with 4% paraformaldehyde (Alfa Aesar, Thermo Fisher Scientific) for 20 min at room temperature, and then rinsed twice and stored in PBS. Cells were blocked and permeabilized with blocking buffer (2% bovine serum albumin and 0.1% Triton X-100 in PBS) for 30 min at room temperature. Cells were incubated with primary antibodies in blocking buffer overnight at 4°C and then washed three times with 0.1% Tween-20 in PBS (PBST). The following primary antibodies and dilutions were used: SOX2 (rabbit monoclonal, Cell Signaling Technology, 3579, 1:200), brachyury (goat polyclonal, R&D Systems, AF2085, 1:150), brachyury (rabbit monoclonal, R&D Systems, MAB20851, 1:400), EOMES (mouse monoclonal, R&D Systems, MAB6166, 1:400), GSC (goat polyclonal, R&D Systems, AF4086, 1:200). Cells were incubated with secondary antibodies (diluted 1:1000): Alexa Fluor 488, 555 or 647-conjugated (Invitrogen Molecular Probes, Thermo Fisher Scientific) and DAPI nuclear stain in blocking buffer for 30 min at room temperature, and then washed twice with PBST and once with PBS.

Wide-field images were acquired on an Olympus IX-70 inverted microscope with a 10×/0.4 numerical aperture objective lens. Tiled image acquisition was used to acquire images of large areas in four channels corresponding to DAPI and Alexa Fluor 488, 555 and 647. Image analysis was performed using custom software in MATLAB as described previously (Yoney et al., 2018). Nuclei segmentation and signal quantification were performed on background-corrected images as follows. The DAPI image was thresholded to generate a binary image separating the foreground (nuclei) from the background. The DAPI image was then filtered with a median and sphere filter with parameters matching the expected size of individual nuclei. Local maxima corresponding to individual nuclei were detected using the MATLAB extended-maxima transform function. Maxima falling within the foreground were used as seeds for watershed segmentation, which was also restricted to the foreground and was used to obtain a labeled object corresponding to each nucleus within the image. The results of the segmentation were used as a mask to obtain the median per cell nuclear intensity in each channel for 5-10 K cells per condition from which histograms were generated.

Chromatin accessibility profiling was carried out using the Omni-ATAC-sequencing protocol as an attempt to reduce mitochondrial DNA in our samples (Corces et al., 2017). We obtained single cells using Gentle Cell Dissociation Reagent (STEMCELL Technologies), rather than Accutase or other enzyme-based reagents, which we found disrupted downstream processing steps. Cells (100,000 per condition) were resuspended in 50 µl lysis buffer (0.1% NP-40, 0.1% Tween-20 and 0.01% digitonin) prepared in resuspension buffer [10 mM Tris-HCl (pH 7.4), 10 mM NaCl and 3 mM MgCl₂ in water] for 3 min on ice. We optimized the lysis time so that the outer membrane was disrupted, as indicated by Trypan Blue staining of the nuclei, but that the nuclei remained intact. The lysis buffer was removed by washing with 0.1% Tween-20 in resuspension buffer, and the transposition was subsequently carried out in 50 µl reaction volume containing 25 µl 2× tagmentation DNA (TD) buffer and 2.5 µl Tn5 transposase (available as individual products upon request from Illuminia). The transposition reaction was incubated at 37°C for 30 min in a thermomixer set to 1000 r.p.m. The transposed DNA was isolated using the Qiagen MinElute kit and was eluted in 30 µl of water in the final step. Library preparation was carried out using 10 µl of transposed DNA per sample as described previously with 10 or 11 cycles of amplification (Buenrostro et al., 2015). Our DNA tended to be under digested, so we performed doubled-sided bead purification using AMPure XP beads (Beckman). Libraries were sequenced as paired-end 75 bp reads, multiplexing all samples per experiment (7) on one lane of the Illumina High-Seq 500 platform at The Rockefeller University Genomics Resource Center. Two biological replicates for each condition were collected, processed and sequenced from independent experiments.

For ChIP experiment in Fig. 2E, we used a SMAD2/3 antibody (R&D systems, AF3797) in 2×10⁶ hESCs that were treated with activin for 2 h or 24 h with or without Wnt priming. We used the Pierce magnetic ChIP kit and followed the procedure provided by the kit. Known SMAD2 binding sites near GSC, NODAL and CER1 genes were selected for subsequent qPCR analysis based on previously published SMAD2/3 ChIP-seq data (Kim et al., 2011; Tsankov et al., 2015). A region near CER1 gene that has low SMAD2/3 ChIP-seq signal was chosen as the negative control. The primer sequences used for ChIP are included in Table S4.

RNAseq data were aligned by Rsubread, and the read count associated with genes were analyzed by featureCounts. The differential fold changes and effective P-values were obtained using DESeq2 (Love et al., 2014). The data in Fig. 1C include those genes that show a significant change (adjusted P-value <0.01) in at least one of the four conditions relative to E7: (1) –/activin, (2) 24 h WNT, (3) WNT/– or (4) WNT/activin. The k-mean clustering was carried out using MATLAB kmeans function. We used a number of clusters (eight in total) such that the final clusters do not change based on the randomized starting conditions for each run of the algorithm. For easier description, we combined three of them into cluster 3 (genes more activated by WNT/activin than other three conditions) and another three into cluster 4 (repressed genes). We provided a list of genes belonging to each cluster (or sub-cluster) in Table S1.

The ChIP-seq and ATACseq data were aligned with bowtie2 for the hg19 reference genome and visualized using deepTools bamCoverage function. The peaks were identified using MACS2 narrow peak calling. The areas underneath the peak (as in Figs 2C,D and 3C,E) were calculated using MultiCovBed function without further normalization. The heatmaps in Fig. 3A,C were generated using deepTools computeMatrix and plotHeatmap. The differential peaks in Fig. 6B were identified using DESeq2 based on the areas underneath each peak (P<0.01). In Fig. 6C, we collected the transcription start site (TSS) and transcription end site (TES) of genes in different categories and examined whether there are any WNT/ACT enhanced ATAC-seq peaks that fell between 20 kb upstream of the TSS and 20 kb downstream the TES. The enriched motifs in the WNT/activin enhanced ATAC-seq peaks were identified by AME (McLeay and Bailey, 2010) with the HOCOMOCOv10 database (Kulakovskiy et al., 2016) against other ATAC-seq peaks that are not sensitive to the WNT/activin condition.

Published datasets used in this study include β-catenin ChIP data (GSM1579346), SMAD2/3 ChIP data from pluripotent hESCs (GSM727557), SMAD2/3 ChIP data from human mesendoderm cells differentiated from hESCs (GSM1505750) and EOMES ChIP data (GSM1505630 and GSM1505631).

We thank current and past members of our group for feedback on this work. We are especially grateful to Dr Riccardo De Santis, Shu Li and Peter Ingrassia for technical support in completing the revisions.

The peer review history is available online at https://journals.biologists.com/dev/article-lookup/doi/10.1242/dev.200335