A wave of WNT signaling balanced by secreted inhibitors controls primitive streak formation in micropattern colonies of human embryonic stem cells

The WNT, BMP and activin/NODAL signaling pathways play a dominant and largely conserved role in vertebrate development. Despite extensive knowledge of the intracellular aspects of these pathways, however, many quantitative aspects of how they create and mediate large-scale patterns remain unknown, especially in the context of human development. A prime example of this is the patterning of the human epiblast, where in a period of approximately two days, and on a millimeter length scale, these three pathways overlap and create the necessary instructions to form and pattern the primitive streak (PS), from which the embryonic germ layers and future anterior-posterior axis of the human embryo emerges. Although mouse studies have shed some light on the involvement of these pathways in the patterning processes, differences in architecture, timing, and juxtaposition of extra-embryonic and embryonic tissues makes the direct comparison difficult (Arnold and Robertson, 2009; Tam and Loebel, 2007). Because of these differences and because the interactions of the aforementioned pathways are complex, with multiple morphogens and secreted inhibitors from overlapping regions acting at the same time, there is a need for an assay that allows for their control in a model human epiblast that also permits a dynamic readout on the scale of single cells.

We have previously shown that human embryonic stem cells (hESCs) confined to micropatterned colonies of 1 mm diameter can be used as an in vitro assay to model the human epiblast (Warmflash et al., 2014; Deglincerti et al., 2016). These micropatterns self-organize in response to WNT, BMP and activin/NODAL signaling and recapitulate the patterning of germ layers observed during mammalian gastrulation. For example, stimulation with BMP4 for 48 h results in concentric rings corresponding to ectoderm, mesoderm, endoderm and extra-embryonic tissue arranged from the center to the edge. As current guidelines prohibit the studies of human embryos after 14 days (the 14-day rule), these models currently remain the best alternative to direct in vivo studies. More importantly, these models allow single cell quantification and control over geometry, density, signaling strength, and genetics.

In subsequent work, we exploited this assay to uncover how the BMP pathway contributes to this patterning (Etoc et al., 2016). Briefly, cells in the colonies are apically-basally polarized and BMP4 receptors are located on the basolateral sides of the cells, restricting access to the apically supplied medium except near the edges. This receptor-mediated prepattern, already present in the pluripotent state, is reinforced by the secreted BMP inhibitor noggin, which in humans is directly induced by BMP4. Owing to diffusion and boundary conditions, noggin is highest in the center and lowest at the colony edges, resulting in an effective gradient of BMP response across the colony.

We have also shown that BMP4 induces WNT ligand and that this WNT is necessary for the mesoderm and endoderm part of the pattern (Martyn et al., 2018). Additionally, we have shown that WNT signaling is sufficient to induce the self-organization of a human PS at the edge of colonies, and that co-presentation of WNT3A and activin leads to the induction of functional human organizer cells that can induce an ectopic secondary axis in chick embryos. This demonstration of another self-organized patterned response to a uniformly presented ligand shows that our system offers an ideal environment in which to explore how WNT signaling leads to spatial organization, and specifically how the human PS forms and is spatially confined.

Here, we address the molecular mechanisms underlying WNT-mediated self-organization of human PS. We show that two primary factors control patterning: E-cadherin (E-CAD; CDH1) and DKK1. First, E-CAD establishes a pre-pattern by limiting the initial WNT response to the boundary. Second, and in parallel to the noggin dynamics in the BMP case, the secreted inhibitor DKK1 is upregulated by a combination of WNT and NODAL signaling and is required to ultimately confine the PS to the colony boundary. Multiple single and double combinations of homozygous CRISPR/Cas9 knockouts of secreted inhibitors of the WNT and NODAL pathways confirmed that only DKK1 plays a major role in the spatial restriction of the PS. We found that cerberus 1 (CER1) is also highly upregulated by a combination of WNT and NODAL signaling, but that in our cells it functions as a NODAL inhibitor rather than dual WNT/NODAL inhibitor. CER1 thus does not influence the size of the PS, but instead serves to bias the mesoderm versus endoderm fate decision in this region. We found also that in DKK1^−/− cells E-CAD not only establishes a pre-pattern, but, via its mutual antagonism with WNT, generates a cooperative epithelial-to-mesenchymal transition (EMT) wave that travels from the micropattern periphery to the center.

We previously showed that uniform application of WNT3A ligand to hESC micropatterns is sufficient for self-organization of a PS-like structure, with mesoderm and endoderm emerging from an EMT on the colony periphery after 48 h and with activin/NODAL level biasing the choice of endoderm versus mesoderm (Martyn et al., 2018) (Fig. 1A). During this time, the transcription factor SOX17 demarcates the endoderm and the transcription factor brachyury (BRA; also known as TBXT) demarcates mesoderm. Changes in the EMT markers SNAIL (SNAI1), E-CAD and N-CAD (CDH2) can also be used to identify the PS, but as these markers are more transient and harder to measure than BRA or SOX17 (Martyn et al., 2018), we decided to use the union of SOX17 and BRA to define the spatial extent of the induced PS. We also showed previously that despite the uniform application of WNT, the interior of the colony remains pluripotent, expressing both NANOG and SOX2 (Martyn et al., 2018). This pattern, with pluripotent cells on the interior and mesoderm and endoderm cells on the periphery, represents the terminal spatial pattern we seek to understand.

In order to decipher the molecular mechanism underlying this spatial pattern, we first attempted to use the detection of nuclear β-catenin (β-CAT; CTNNB1) as an early readout for canonical WNT signaling (Riggleman et al., 1990; Valenta et al., 2012; Cadigan and Peifer, 2009). However, commercially available antibodies did not have adequate resolution on our dense epithelia (Fig. S1A), so instead we used LEF1, a co-factor of β-CAT that is a direct target of WNT signaling with the same response profile as AXIN2 (Fig. S1B) and is localized to the nucleus (Cadigan and Waterman, 2012; Estarás et al., 2014). We found that the LEF1 response profile depended on colony density, with nuclear expression throughout the colonies at low density and restriction to the periphery at high density (Fig. 1B). Co-presentation of the SMAD2 pathway inhibitor SB-431542, together with WNT3A did not change the outcome, demonstrating that the density dependence of the LEF1 pattern is specifically due to WNT, and not caused by secondary activin/NODAL signaling.

Similar to previous work studying the effect of colony density on BMP signaling (Etoc et al., 2016), low density represents a pre-epithelial state before tight junctions have completely established themselves, whereas high density represents an epithelial state with complete tight junctions (Fig. 1C) where the SMAD1 response due to BMP stimulation would also be edge restricted (Fig. S2C). Because density is an important variable, for all experiments in this study we consistently use two defined values of seeding density: ‘low’ or ‘high’ (Fig. S2B).

As one of the factors involved in BMP4-induced self-organization was polarized signal reception, we first examined the localization of the WNT receptors (MacDonald and He, 2013) FZD7 and LRP5/6 in our micropatterns. We found that although some WNT receptors were detected on the apical side, they were predominantly and homogeneously located basolaterally underneath the tight junctions (Fig. 1D). We also found little distinction between edge and center. To test functionally for signal reception, we cultured cells on Transwell-filter culture dishes, in which cells can be selectively stimulated from the apical or basal side. Cells were cultured at the same density as the high-density micropatterns and were stimulated with WNT3A for 12 h. A stronger response was detected with basal than with apical stimulation, but the mean basal response on filters fell below the edge response on colonies (Fig. 1E,F), suggesting that additional factors were involved in setting up the WNT response on the colony boundary.

Given the fact that E-CAD is classically considered to antagonize WNT signaling via its binding and sequestration of β-CAT (Nelson and Nusse, 2004; Heuberger and Birchmeier, 2010; Fagotto et al., 1996; Ozawa et al., 1989), and the fact that cells on the periphery of our micropatterns have fewer neighboring cells and so presumably fewer E-CAD junctions, we hypothesized that E-CAD may contribute to the early WNT pattern. To test this hypothesis, we made clonal CRISPR/Cas9 E-CAD knockout cell lines. We used a guide RNA targeting a region present in all isoforms of the protein and confirmed the result with western blot and immunofluorescence using two separate antibodies (Fig. S2F-H). We found that E-CAD^−/− hESCs could be passaged and seeded as per normal hESC culture, grew at the same rate as wild-type cells, maintained pluripotency markers as single cells, unpatterned colonies, and micropatterned colonies, and still apical-basally polarized at high density to form intact epithelia (Fig. S2A-D). E-CAD^−/− WNT receptor localization was indistinguishable from that of wild type (Fig. S2E), and the response to BMP4 continued to be observed at the edge only. Interestingly, N-cadherin (N-CAD) protein was upregulated in the knockout (Fig. S3), and in other contexts can substitute for the loss of E-CAD (Libusova et al., 2010).

The most striking phenotype of the E-CAD^−/− cells, however, was that the early WNT pattern was abolished in E-CAD^−/− micropatterns (Fig. 2C). Cells at the center of high-density colonies now showed nuclear localization of LEF1. Quantifying our results over multiple micropatterns, we observed that the WNT response is modestly biased to the edge, and generally comparable to the level in low-density colonies (Fig. 2D). Our Transwell filter assay still showed an apical-basal asymmetry, but now the basal response was elevated to match the edge response of the parental cells (Fig. 1E,F). These results demonstrate that the early WNT pattern is primarily due to E-CAD activity, with a minor influence exerted by WNT receptor accessibility.

If all cells express E-CAD why are there spatial differences in the WNT response? We hypothesized that spatial differences in E-CAD localization or in the state of E-CAD junctions and their binding partners could account for spatial WNT signaling differences. E-CAD and β-CAT stains support this hypothesis as they show that in high-density wild-type micropatterns E-CAD is reduced in cells on the periphery and there is observable cytoplasmic β-CAT here as well (Fig. 2A,B). Actin stress fibers have also been observed on hESC micropattern boundaries (Rosowski et al., 2015; Przybyla et al., 2016) and have been implicated in E-CAD dysregulation (Kanellos et al., 2015). Phalloidin staining in our wild-type micropatterns revealed that there were indeed actin stress fibers and they were restricted to the region that is WNT responsive. Furthermore, these stress fibers were absent from low-density or E-CAD^−/− micropatterns (Fig. 2A), which suggests a connection between mechanics and WNT signaling.

To test further whether the classical connection between E-CAD, β-CAT and the actin cytoskeleton was responsible for our early WNT pattern, we performed two additional experiments. In the first, we inserted into our E-CAD^−/− cells either constitutively expressed full-length E-CAD or constitutively expressed E-CAD that lacked the β-CAT-binding domain (Gottardi et al., 2001) into the AAVS1 locus using transcription activator-like effector nucleases (TALENS). Clonal lines were cultured in micropatterns and stimulated for 12 h to examine the WNT response. We found that the constitutively expressed full-length E-CAD rescued the edge-restricted phenotype but that the E-CAD without the β-CAT-binding domain did not (Fig. 3A,B). This shows that E-CAD binding to β-CAT is essential for the early WNT pattern.

In the second experiment, we disrupted the actin cytoskeleton across the entire micropattern with the small molecule inhibitors blebbistatin or cytochalasin B while stimulating with WNT3A. Blebbistatin has been shown to dislodge E-CAD from the membrane into the cytoplasm in hESCs in pluripotency (Li et al., 2010), and cytochalasin B acts even more directly by dissociating the actin cytoskeleton. We found that both reduced the edge restriction, with blebbistatin broadening the size of the LEF1 band and cytochalasin B permitting a WNT response even in the center of the colony (Fig. 3C). Taken together, our results demonstrate that colony geometry acts via the cytoskeleton and E-CAD to bias WNT signaling to the colony boundary.

Having understood the reasons for the edge asymmetry in the initial response to WNT, we were surprised to see that the location of the PS was virtually the same in wild-type and E-CAD^−/− colonies at 48 h (Fig. 4A). Examination of LEF1, SOX2 and BRA expression at intermediate times showed that a homogeneous early expression of these markers gradually becomes restricted to the edge (Fig. 4B). Because cells in these colonies continue to grow and divide throughout this time course, we checked whether increasing density could be responsible for this effect. However, E-CAD^−/− colonies at a higher starting density (matching that of the Fig. 4A colonies after 36 h stimulation) also showed a WNT response in the center (Fig. S3C). Thus, we can rule out cell proliferation and increasing density as a major contributing factor for the progressive edge restriction of the WNT response, and instead note that these dynamics are suggestive of a WNT-induced secreted inhibitor of WNT that is highest in the center and progressively restricts WNT activity to the boundary.

We had previously shown that BMP4 directly induced the expression of its own inhibitor, noggin, which in turn was necessary and sufficient to restrict BMP signaling to the colony edge after 48 h of stimulation (Etoc et al., 2016). To determine whether a similar mechanism of WNT3A inducing its own inhibitor was involved in this case as well, we activated the pathway with CHIR-99021, a small molecule compound that acts cell-autonomously and will not interact with receptors and secreted inhibitors. After 48 h stimulation, the compound edge restriction was abolished (Fig. 4C). This result strongly suggests the involvement of secreted inhibitors in WNT-mediated self-organization and PS formation.

To test this hypothesis further, and identify the relevant inhibitors, we focused on WNT inhibitors for which loss of function leads to early gastrulation defects phenotypes in the mouse. These include SFRP1, SFRP2, DKK1 and DKK3 (Cruciat and Niehrs, 2013). Additionally, as WNT activates NODAL (Martyn et al., 2018), and these two ligands have been shown to act synergistically to induce mesendodermal genes (Funa et al., 2015), we also included LEFTY1, LEFTY2 and CER1 on our list. qPCR was used to assess the induction of these inhibitors when cells were treated with WNT3A alone, or WNT3A+SB-431542 (SB) to distinguish direct versus indirect induction. After 12 and 24 h of stimulation, expression of SFRP1, SFRP2 and DKK3 remained relatively unchanged regardless of WNT3A or WNT3A+SB treatment (Fig. 4D). DKK1 expression, however, was highly upregulated in response to WNT3A. Similar to what has been previously reported (Funa et al., 2015), this appears to depend on synergy between the WNT and NODAL pathways, as DKK1 induction is lower in WNT3A+SB conditions. A stronger dependence on SMAD2 signaling was observed for the WNT3A induction of CER1 expression. Finally, the expression of LEFTY1 and LEFTY2 depend even more on NODAL signaling as they were also downregulated in response to WNT3A+SB and were not activated with WNT3A alone. Thus, DKK1 and CER1 emerged as the leading candidates involved in WNT self-organization.

To test whether DKK1 and CER1 are required for the late time, 48 h WNT pattern, we used CRISPR/Cas9 to generate DKK1 and CER1 knockouts. Each clonal line was stimulated with WNT3A and compared with wild type. After 12 h of stimulation, DKK1^−/− colonies appeared similar to the control (Fig. 5A). After 48 h of treatment, however, the size of the PS was dramatically increased compared with the wild-type and E-CAD^−/− lines, with only a small center of SOX2⁺ undifferentiated cells remaining (Fig. 5B,C). This result, which was confirmed in two additional clonal DKK1^−/− lines (Fig. S4A), demonstrates that DKK1 activity is required for WNT-mediated patterning, and is consistent with a reaction-diffusion model.

To test this interpretation further and see whether human DKK1 protein can protect cells from WNT ligand at a distance in a non-cell-autonomous manner, we created a clonal RUES2 cell line that expresses human DKK1 tagged with V5 epitope under the control of a doxycycline promoter (Fig. S4G). When these RUES2-DOX:DKK1-V5 cells were seeded sparsely into E-CAD^−/− micropatterns and made to express DKK1, we observed that they could block BRA expression in cells up to approximately five cell lengths away from them, thus demonstrating that human DKK1 can act as a long-range diffusible WNT inhibitor (Fig. S4H,I).

To test for epistasis between WNT inhibition at early times mediated by E-CAD, and at late times by DKK1, we generated a double E-CAD^−/−DKK1^−/− knockout line. In response to WNT3A stimulation for 48 h, all cells in the micropatterned colonies differentiated, with no SOX2⁺ cells left in the center (Fig. 5B,C). This suggests that DKK1 and E-CAD are the two major players among the collection of WNT inhibitors that block differentiation in our micropatterns.

Comparison of the expanded PS in the single and double knockouts with wild type and classification of cell types (see Materials and Methods and Fig. S5A,B) also established that both the total number and the ratio of BRA to SOX17 cells changed (Fig. 5D). Whereas the proportion of BRA⁺/SOX17⁻ cells in the RUES2 wild-type line was ∼10%, in DKK1^−/− cells it doubled to 20%, and in E-CAD^−/−DKK1^−/− it tripled to ∼30%. The fraction of cells expressing both BRA and SOX17 also greatly increased in the mutant lines. This suggests that, in addition to determining the size of the PS, DKK1 may also be involved in the segregation of mesodermal and endodermal fates.

When stimulated with WNT3A, CER1^−/− cells did not show any change in the size of the PS domain in comparison with wild type (Fig. 5E). There was, however, a significant shift in the proportion of mesodermal versus endodermal fates. Unlike the DKK1^−/− or E-CAD^−/−DKK1^−/− cells, this time the shift was towards greater endoderm, with almost all differentiated cells expressing SOX17 and none expressing BRA (Fig. 5D). This represents a similar phenotype to the WNT3A+activin treatment (Fig. 1A), and prompted us to investigate the status of NODAL/activin signaling. We found that SMAD2 signaling is significantly enhanced in the CER1^−/− knockout line (Fig. 5F,G), penetrating farther into the colony from the edge and with a higher nuclear-to-cytoplasmic ratio. It is known that in mouse CER1 inhibits BMP and NODAL signaling but not WNT signaling (Belo et al., 2000) (at odds with other vertebrate systems in which its tri-functionality motivates its name; Cruciat and Niehrs, 2013). In human, however, it is only known that CER1 inhibits NODAL and a subset of BMP ligands, with a verdict on WNT inhibition still awaited (Aykul et al., 2015; Aykul and Martinez-Hackert, 2016). As the size of the PS remains unchanged while SMAD2 signaling increases along with the proportion of endodermal cells in the CER1^−/− colonies, our results suggest that in hESCs CER1 acts primarily as a NODAL inhibitor rather than a WNT inhibitor.

As it was previously shown in the mouse that the most dramatic CER1 phenotype is observed when it is doubly knocked out with LEFTY1 (Perea-Gomez et al., 2002), we also generated CER1^−/−LEFTY1^−/− and LEFTY1^−/− clonal cell lines. We detected no difference in the WNT-response phenotype between wild-type and LEFTY1^−/− cells, and no difference between CER1^−/− and CER1^−/− LEFTY1^−/− cells (Fig. S4B,C). In order to check for all other players identified in our RNA-seq and qPCR results, we generated DKK3^−/− and SFRP1^−/−SFRP2^−/− clonal cell lines. None of these lines displayed any phenotypic difference compared with wild type (Fig. S4D-F). We conclude that DKK1 and CER1 are the major secreted inhibitors that control WNT patterning in our model system.

The size of the PS region in the DKK1^−/− cell line at 48 h is intermediate between the smaller wild-type PS region and the fully converted PS region of the double E-CAD^−/−DKK1^−/− cell line. Given that in both RUES2 and DKK1^−/− colonies WNT signaling begins at the edge (Fig. 5A), an important and relevant question is whether the 48 h result is at steady state, or if given more time the PS would eventually expand inward and consume the entire colony. To address this question, we fixed and stained wild-type and DKK1^−/− micropatterns at 12, 24, 48 and 72 h. We found that although differentiation starts similarly for both, the wild-type micropatterns seemed to reach a steady state by 24-48 h whereas differentiation and EMT in the DKK1^−/− colonies continued to proceed inwards, eventually almost consuming the entire micropattern by 72 h (Fig. 6A). This is consistent with a wave of WNT differentiation proceeding from the outer edge to the center.

To confirm the existence of this wave and study it further, we knocked out DKK1 in a previously described RUES2-GLR (Germ Layer Reporter) cell line in which ectoderm (SOX2), mesoderm (BRA) and endoderm (SOX17) germ layers are tagged with three separate fluorescent markers (Martyn et al., 2018). This enabled us to evaluate the change in differentiation and fate acquisition in the same micropattern across different times. Stimulation of the control RUES2-GLR with WNT3A led to a downregulation of SOX2 and an upregulation of SOX17 that began at the edge but stopped a few cell layers in, maintaining the PS at the periphery. RUES2-GLR colonies in which the DKK1^−/− mutation had been introduced, however, displayed a wave of progressive downregulation of SOX2 and upregulation of SOX17, which, as in the stained time-course, began at the outer edge and did not stop (Fig. 6B,C).

What is the mechanism for such a wave? Given that E-CAD expression recedes as the differentiation front advances (Fig. 6A), we posit that the wave results from a positive EMT feedback loop. WNT3A is known to induce EMT and to downregulate E-CAD through the transcription factor SNAIL (Batlle et al., 2000; Cano et al., 2000). As illustrated in Fig. 6D, WNT3A first downregulates E-CAD in WNT-susceptible edge cells and causes them to go through EMT. In so doing, these cells destabilize their E-CAD junctions with their neighbors. This leads to a domino-like propagation of EMT from cell to cell via shared cell contacts. If the differentiated cells were induced to secrete a diffusing and thus long-range WNT inhibitor, i.e. DKK1, this would accumulate in the center more than at the edges and the wave would be expected to halt (Etoc et al., 2016). This is how the wild-type cells achieve the observed steady state, and why the DKK1^−/− and GLR:DKK1^−/− cells fail to do so.

In addition to these dynamics, previous results have shown that WNT ligand also upregulates WNT production in hESCs (Martyn et al., 2018). To test whether this endogenous WNT signaling contributed to the observed dynamics, we compared wild-type and DKK1^−/− micropatterns with and without IWP2, a small molecule that blocks all WNT ligand secretion. We found no significant differences (Fig. S6A), most likely due to the fact that we were already stimulating our micropatterns with a high dose of WNT in the media and thus created a saturated regime in which endogenous WNT would not make any significant contribution to the dynamics. We also tested the involvement of activin/NODAL signaling (which has a baseline activity in our media) in this wave by comparing DKK1^−/− micropatterns stimulated with either WNT3A or WNT3A+SB. Consistent with other studies of EMT (Derynck et al., 2014), our wave stopped when we blocked activin/NODAL signaling with SB (Fig. S6B).

The spatial pattern of WNT signaling in our colonies is defined by inhibition by E-CAD at the earlier stages and DKK1 at later stages. These inhibitors operate very differently in space: E-CAD bridges adjacent cells, whereas DKK1 diffuses across the colony and also leaks out at the edges (Etoc et al., 2016). Downstream of WNT, NODAL and CER1 are produced and, together with WNT, generate mesoderm and endoderm fates. To unravel the complexity involved in this process fully, we formulated a quantitative dynamic model.

A good model uses a portion of the data to fit parameters and then make testable informative predictions about the remainder of the data, with as few variables as possible. With these criteria in mind, we define a 2D partial differential equation (PDE) model in which the intracellular WNT signal, W(r,t), is normalized to [0,1] and a simple Michaelis–Menten system of equations links it to DKK1 and E-CAD. One advantage of our formulation is that it is separable, so we can fit the DKK1-specific parameters to the E-CAD^−/− data and vice versa (Fig. 7A). Because we cannot measure WNT levels directly, we use the immunofluorescence data for LEF1 at 12 h and the percentage loss of SOX2 at 48 h as surrogates, after normalizing to [0,1]. For the full list of equations, descriptions of the variables, and initial and boundary conditions, the reader is directed to the supplementary Materials and Methods.

Fig. 7C shows the quality of the fits to the immunofluorescence data at 12 and 48 h in the two knockout lines, which we consider to be acceptable (Figs S9 and S11). Our model with no additional adjustments is then able to predict and reproduce the 72 h data in Fig. 7C as well as the EMT wave (in the DKK1^−/− background with endogenous E-CAD), which is most visible in the GLR:DKK1^−/− cell line (Fig. 7D). This live data is somewhat variable owing to density variation and phototoxicity (see supplementary Materials and Methods) but the prediction of the shape of the wave front and where it should be at a specific time closely follows the observations. The collapse of the wave after it reaches a radius of about 200 μm is also predicted and the model agrees with both the live data and the immunofluorescence data in Fig. 6A.

The model demonstrates in an explicit and mathematically precise manner that the reciprocal inhibition between WNT and E-CAD (Fig. 7A) can give rise to bistability, and how, in the absence of the long-range secreted inhibitor DKK1, the bistability resolves by an inward-propagating wave that eliminates the epithelial state in favor of the PS mesenchyme. The model situates the WNT system within a general class of problems in which waves result from bistability. The interest for development, elaborated in the Discussion, is that waves propagate information faster than diffusion, and the proposed mechanism is generic and largely parameter independent.

Having captured the WNT response dynamics for the different cell lines, we can further model whether the differentiated cells become endoderm or mesoderm. To do so, we assume that NODAL favors endoderm over mesoderm (Martyn et al., 2018), and we describe this as a branching probability from an intermediate mesendoderm state (Fig. 7B), as we lack more quantitative data about the genetic network for mesendoderm specification. As a surrogate for NODAL signaling, we used measured SMAD2 profiles at 24 h for each cell line (Fig. 5G and Fig. S12 in supplementary Materials and Methods). Fig. 7E shows the comparison of model and data. The conversion of mesoderm to endoderm as a result of upregulation of NODAL is clearly visible in the comparison of CER1^−/− with wild type or E-CAD^−/−.

Untangling the dynamics of early cell fate specification, and pattern formation during gastrulation remains a formidable challenge, especially during early human development. Using highly quantitative and standardized human gastruloids to model human gastrulation, we had previously shown that the cascade of signaling begins with BMP4, which is sufficient to induce gastrulation and pattern all germ layers within 48 h. We subsequently established that the mechanism underlying this BMP self-organization of gastrulation is based on two components: TGFβ receptors are only exposed to apical signals in a pre-pattern around the edge of the colony, and a Turing-like reaction-diffusion mechanism mediated by noggin, which specifically in humans is induced directly by BMP4. We subsequently showed that the WNT pathway, operating downstream of BMP4, is necessary and sufficient for the mesoderm and endoderm layers and can self-organize its own pattern. In this work, we have unraveled the molecular mechanism underlying this WNT-mediated self-organization.

This mechanism includes several elements. The first element is the ‘pre-pattern’ of WNT sensitivity. Unlike the situation with BMP4, this pre-pattern is not imposed by polarity of receptor localization, but rather by β-CAT mechanosensation via E-CAD, and the cytoskeleton, which is governed by tissue geometry. Given the importance of this pre-pattern for WNT patterning in a model human epiblast, how applicable might such pre-patterns be to gastrulation more generally? Gastrulation takes many forms across all animals, but a common theme is the correlation of WNT signaling with invagination (Arendt, 2004) (such as through the PS in amniotes or through a blastopore in lower orders). Invagination necessarily involves a breaking of geometric symmetry and a concentration of mechanical forces in an area of high curvature, and it is tempting to think that β-CAT mechanosensation is also involved here in either pre-patterning or reinforcing a WNT signal. Other notable examples of β-CAT mechanosensation as a crucial element for WNT patterning include Drosophila gastrulation (Farge, 2011), Xenopus embryonic explants (Farge, 2011) and avian follicles (Shyer et al., 2018).

The second element in the WNT patterning process is the Turing-like activator-inhibitor pair of WNT and DKK1. In a manner similar to BMP4, which in a species-specific manner induces noggin, the WNT ligand also directly induces the expression of its own inhibitor: DKK1. In the mouse, the anterior visceral endoderm (AVE) on the opposite side of the embryo from the PS is traditionally thought of as the major source of DKK1, with this inhibitor being produced in a WNT-independent manner (Arnold and Robertson, 2009; Kimura-Yoshida et al., 2005). The expression of DKK1 by the AVE makes disentangling any intrinsic contribution of the epiblast to the position and control of the PS in the embryo difficult. Given that our synthetic epiblast achieves a stable pattern via the induction of DKK1 by WNT without external sources, however, we might expect to see DKK1 produced by cells in the PS. This is exactly what recent single cell and FACS-sorted RNA-seq data show (Peng et al., 2016; Simunovic et al., 2018 preprint). Similarly, results from studies in rabbit gastrulation show that DKK1 is also expressed in cells in the PS, especially in epiblast cells undergoing EMT (Idkowiak et al., 2004). Although expression of activator and inhibitor on opposite sides of an embryo (such as in PS and AVE) seems logical and necessary, the Turing mechanism requires production of the inhibitor where the activator is highest, because the inhibitor must also spread rapidly to confine the activator.

We discovered that the same is true for CER1, the third major element in the WNT patterning process. Despite also being thought of as an AVE product, single cell RNA-seq studies have shown that CER1 is also expressed in the mouse and rabbit PS (Peng et al., 2016; Idkowiak et al., 2004). In our human model, we find that CER1 acts as a NODAL inhibitor and controls the balance between mesodermal and endodermal fates emerging from the back of the EMT wave. Thus E-CAD, DKK1 and CER1 are the primary determinants of WNT-induced PS patterning in artificial human epiblasts.

Our knockouts of CER1 and DKK1 in hESCs also reveal possible species-specific differences between human, mouse, chick and rabbit. For example, a difference between human and other studied species is that the knockout of CER1, or even the double knockout of CER1 and LEFTY1, does not lead to an expansion of the PS in our model system. This is surprising as in mouse it leads to multiple streaks (Perea-Gomez et al., 2002) and in the chick NODAL repression by cerberus is thought to be the dominant repressor of ectopic streaks (Bertocchini and Stern, 2002). Another difference is that the mouse Dkk1 knockout does not broaden the PS as defined by brachyury (Bra) expression, though it does broaden the expression of a synthetic Wnt reporter (Lewis et al., 2008). This could suggest that PS formation in human is dominated more by WNT signals whereas in mouse the PS and Bra require other signals (such as Nodal) as well. Indeed, recent experiments with mouse epiblast-like cells on micropatterns (Morgani et al., 2018) suggest that there is a stronger dependence on Nodal signaling to achieve a PS-like region and Bra-positive cells than in hESCs. We hope that future comparative experiments will shed light on these apparent species-specific differences, and will also test for possible confounding in vitro versus in vivo differences as well.

Our results for the DKK1 and CER1 knockouts have implications not just for the spatial control of the streak, but also for mesoderm versus endoderm cell fate decisions within the streak. We found a higher ratio of mesoderm to endoderm in DKK1^−/− cells, which is consistent with the reduction in the mouse of anterior endoderm in favor of mesoderm following the analogous knockout (Lewis et al., 2008). Conversely, removal of CER1 favors endoderm over mesoderm in our system, with upregulation of nuclear SMAD2 compared with wild type. Although less well understood than the appearance of multiple streaks, the elimination of Cer1 and Lefty in the mouse embryo also results in a marked conversion of Bra cells to more Nodal-regulated fates, such as anterior endoderm and axial mesoderm (Perea-Gomez et al., 2002). Based on these observations, as well as the result that WNT induces both DKK1 and NODAL in hESCs, we speculate that mammalian anterior endoderm is specified by a transient WNT3 pulse terminated by self-induced DKK1 and followed by upregulation of NODAL signaling. In fact, a transient application of WNT followed by its removal and the addition of activin is precisely the most efficient and commonly used endoderm differentiation protocol used for hESCs (Deglincerti et al., 2016; D'Amour et al., 2006; Teo et al., 2014; Yoney et al., 2018).

One of the more intriguing results to emerge from our investigation was the presence of a propagating EMT wave or front from the edge of our micropatterns towards the center. This wave is generated by downregulation of E-CAD by EMT and the negative regulation of WNT signaling by E-CAD, but also requires activin/NODAL activity as SB treatment blocks the wave. As DKK1 acts as a negative feedback to halt the wave, the wave is much longer and more apparent in the DKK1^−/− micropatterns than in wild-type micropatterns. This is not to say, however, that wild-type cells could not undergo such a wave in vivo. For example, removal of the DKK1-expressing cells from the wave front, as happens in migration of the mesenchymal cells out of the PS in vivo, would allow the wave to continue progressing. This cannot take place in our system because cells are confined to the surface of the micropattern. We note that although it is harder to cleanly isolate, in mouse and chick embryos there is evidence for cooperative EMT in maintaining and extending the PS (Voiculescu et al., 2014; Williams et al., 2012).

Traveling waves of activity have been seen in several other developmental contexts, but it is important to distinguish propagation of a front separating two states with phase waves, such as in somitogenesis, that spatially modulate the phase of an underlying autonomous cellular oscillator (Aulehla and Pourquié, 2010). Waves with a propagating front have been seen in several other developmental contexts, such as calcium waves following fertilization or in large embryos that presumably function to synchronize tissues (Jaffe, 2008; Vergassola et al., 2018). A wave of mitotic activity in frog extracts has also recently been observed and linked to bistability in the CDK system (Chang and Ferrell, 2013). Why might a wave be useful for patterning? A wave is rationalized as a way to spread information more rapidly than diffusion in a large system, and the hundred-micron-scale disk-shaped epiblast in rabbit and humans may require such a solution. Patterning via a wave may be a widespread mechanism for tissue patterning as it only requires a bistable system and some means for the favored state to spread between cells.

Our hESC micropatterned PS model provides a platform in which complex and developmentally relevant patterning dynamics can be followed, quantified and, ultimately, deconstructed. It thus provides a means of understanding how embryonic development can be both robust and flexible, and to ultimately apply this understanding to engineer systems to further regenerative medicine.

All experiments used either the RUES2 hESC cell line, CRISPR/Cas9-edited cell lines based on this line, or CRISPR/Cas9-edited versions of the previously described RUES-GLR cell line (Martyn et al., 2018). hESCs were maintained in HUESM medium conditioned by mouse embryonic fibroblasts (MEF-CM) with additional 20 ng/ml bFGF. Mycoplasma testing was carried out before beginning experiments and again at regular 2-month intervals. Cells were grown on GelTrex (Invitrogen; 1:40 dilution)-coated tissue culture dishes (BD Biosciences) for maintenance and expansion. Dishes were coated overnight at 4°C and then incubated at 37°C for at least 20 min before the cells were seeded onto the surface. Passaging was performed using Gentle Cell Dissociation Reagent (Stem Cell Technologies, 07174).

Micropatterns were created using glass coverslips from CYTOO. The coverslips were first coated with 10 μg/ml laminin 521 (Biolamina) then diluted in PBS with calcium and magnesium (PBS++) for 3 h at 37°C. hESCs were dissociated with StemPro Accutase (Life Technologies) for 7 min and then washed once with growth media, washed again with PBS, and then re-suspended in growth media with 10 μM ROCK inhibitor (Y-27632, Abcam) in 35 mm tissue-culture plastic dishes. For each coverslip, 1×10⁶ cells in 2 ml of media were used. After 1 h, the ROCK inhibitor was removed and was replaced with standard growth media supplemented with Pen-Strep (Life Technologies). Cells were stimulated with the following ligands or small molecules 12 h after seeding: 100 ng/ml WNT3A, 100 ng/ml activin A, 10 μM SB-431542 (Stemgent), 6 μM CHIR-99021 (EMD Millipore).

Cells were fixed in 4% paraformaldehyde for 20 min at room temperature, washed twice and stored in PBS. Permeabilization (30 min at room temperature), primary antibody incubation (overnight at 4°C), and secondary antibody incubation (30 min at room temperature) were carried out with blocking buffer (3% donkey serum and 0.1% Triton X-100 in PBS). After the primary incubation, cells were washed three times with PBST (PBS+0.1% Tween-20) for 30 min each. After secondary incubation, cells were washed twice more with PBS, and then mounted on glass slides for imaging. All of the primary antibodies used are listed in Table S1. Secondary antibodies were donkey anti-rabbit, -mouse or -goat antibodies conjugated with Alexa Fluor 488, 555 or 647 (Life Technologies).

All images used were acquired with a Zeiss Axio Observer using ZEN software and a 20×/0.8 numerical aperture (NA) lens. All image analysis was conducted using custom software written in Matlab. Images were first background subtracted and normalized and then stitched on a colony-by-colony basis. Images for background subtraction and normalization were acquired in the spaces between colonies where no cells were present. For segmentation of individual cells, we first used Ilastik classification to separate foreground from background. The classifier was trained for each experiment on the DAPI images of four randomly chosen stitched colonies from that experiment. Once foreground and background were obtained, the DAPI channel was then filtered with a median and h-max filter and subtracted against a gradient of the image in order to identify the nuclei centers. These centers were then used as seeds for a watershed, against which the background mask was applied to obtain the final segmentation. Using this segmentation mask, we then obtained average intensities for each cell of the nuclear markers in the other channels. For radial plots, the intensity of immunofluorescence signal for each marker was normalized to the DAPI intensity, and these corrected single cell expressions were then radially binned and averaged. The final radial profile represents the average of the indicated number of colonies. For clustering and classification, single cell intensities were log transformed and then clusters were fitted with a Gaussian mixture model (see Fig. S5 for example).

RNA was collected in Trizol at the indicated time points from either micropatterned colonies or from small unpatterned colonies. Total RNA was purified using the RNeasy mini kit (Qiagen). qPCR was performed as described previously (Etoc et al., 2016), and primer sequences are listed in Table S2. RNA-seq data were taken from a previous publication (Etoc et al., 2016), and all raw data are available from the GEO database, accession number GSE77057.

We used Costar Transwell 24-well plates with 0.4 μm pore-sized clear polycarbonate membrane inserts (Fisher Scientific, 07-200-147). Membranes were coated with 10 μg/ml laminin 521 (Biolamina) diluted in PBS++ for 3 h, followed by washing three times with PBS++. Single cells were collected and seeded as per micropattern protocol. To image the membrane, the Transwell was removed from the multi-well plate after fixing and staining and placed on top of a coverslip.

Full-length E-CAD and E-CAD lacking the β-CAT binding domain were amplified by PCR from hE-cadherin-pcDNA3 and hE-cadherin/Δβ-catenin-pcDNA3 plasmids (Addgene plasmids #45769 and #45772) and inserted downstream of a pCAG promoter and puromycin resistance cassette flanked by 1 kb homology arms for the AAVS1 safe harbor locus. This plasmid and TALENS targeting the AAVS1 site were then nucleofected into 1×10⁶ pluripotent E-CAD^−/− cells using the B-016 setting on an Amaxa Nucleofector II (Lonza). Nucleofected cells were then plated as per maintenance conditions, but supplemented with 10 μM ROCK inhibitor. Selection for puromycin commenced after 2 days, and ROCK inhibitor was maintained until colonies reached adequate size (typically 8-16 cells per colony). To derive pure clones, individual colonies were picked in an IVF hood with a 20 μl pipette tip and seeded into separate wells with growth media and ROCK inhibitor. Once successfully established, each clone was assayed functionally for brightness and homogeneity of the overexpressed E-CAD. The RUES2-DOX:DKK1-V5 cell line was obtained in an identical fashion, except for the use of a TRE promotor instead of a pCAG promoter and an rTA-T2A-puromycin element instead of only puromycin. The V5 tag was attached to the C-terminal end of DKK1, and DKK1 itself was obtained from a commercially available cDNA plasmid (Thermo Fisher Scientific, MHS6278-202801665).

The CRISPR/Cas9 system was used to generate all knockout lines. Three rules of reproducibility and quality control were applied: (1) at least two independent clones for each gene clonal lines were isolated and studied in parallel; (2) lack of off-target effects was assessed by qPCR; and (3) their ability to maintain the pluripotent state as assessed by expression of NANOG, OCT4 (POU5F1) and SOX2. For the E-CAD, DKK1 and CER1 knockouts, one clone of each was also assessed for chromosomal integrity by karyotyping. The sgRNA target for each gene is listed in Table S3. A list of all knockout lines is given in Table S4. The sgRNAs were cloned into a pX330 plasmid (Cong et al., 2013) that we modified to co-express a puromycin-2A-EGFP cassette, as this strategy gave us a higher percentage of successfully targeted clones. Transfection was carried out using the B-016 setting of a Nucleofector II instrument and using the Cell Line Nucelofector Kit L (Lonza). Transfected cells were immediately seeded in ROCK inhibitor on GelTrex-coated culture dishes, and puromycin was added after 24 h for 24 h. Cells were then passaged as single cells using Accutase (Stem Cell Technologies) and sparsely seeded to facilitate picking individual clones. Clones were handpicked with a 20 μl pipette tip and, once expanded, genomic DNA from each clone was extracted with the DNeasy Blood & Tissue kit (Qiagen). The locus for each targeted gene was then PCR-amplified using the primers listed in Table S4, and submitted for Sanger sequencing. The resulting chromatograms for each clone were decomposed using the TIDE webtool (Brinkman et al., 2014; http://tide.nki.nl). Only clones that showed a high probability for both alleles of the gene of interest having a missense mutation leading to a premature stop codon were kept. Examples of these mutations are shown in Fig. S7A. As a check on the integrity of our most critical knockout cell lines, the E-CAD^−/− clone 1, DKK1^−/− clone 1 and CER1^−/− clone 1 cell lines were additionally sent for karyotyping and were found to be karyotypically normal (Fig. S7B). In the case of E-CAD, we note that although there are four recorded isoforms (https://www.ncbi.nlm.nih.gov/gene/999) our sgRNA target is after the start codon of each and is present in all of these isoforms. As an additional check to ensure that our CRISPR/Cas9 mutation results in a complete knockdown instead of a hypomorph, we note that the E-CAD antibody we use to check for E-CAD presence using immunofluorescence (Fig. 2A) was produced by immunization with a synthetic peptide corresponding to the sequence surrounding Pro780 of human E-cadherin [see datasheet for Cell Signaling Technology E-Cadherin (24E10) rabbit mAb; https://media.cellsignal.com/pdf/3195.pdf]. This corresponds to exon 15, which is near the C terminus of the final protein and is shared between all the E-CAD isoforms present in our cells (Fig. S2F). Immunofluorescence tests with an antibody from a different supplier that targets a different conserved region (produced by immunization with synthetic peptide within human E-CAD), aa 600-700 [see datasheet for Abcam anti-E-Cadherin (EP700Y) rabbit mAb; https://www.abcam.com/e-cadherin-antibody-ep700y-intercellular-junction-marker-ab40772.html], also gave the same result (Fig. S2G).

See supplementary Materials and Methods and Figs S8-S14 for details of the mathematical modeling.

hE-cadherin-pcDNA3 and hE-cadherin/Δβ-catenin-pcDNA3 were gifts from Barry Gumbiner (Addgene plasmids #45769 and #45772). We thank members of our lab for constructive comments, especially Sam Khodursky and Fred Etoc, and we thank the staff of the Rockefeller Bio-Imaging Resource Center (BIRC) for their guidance in live imaging of the GLR reporter cell lines.