Prediction and control of symmetry breaking in embryoid bodies by environment and signal integration

In the early mouse embryo, the blastocyst implants into the uterus wall on the fourth day of gestation. A series of events follows, including formation of an anteroposterior (AP) axis initiated by the localization of anterior visceral endoderm (AVE) cells at the prospective anterior side and emergence of the primitive streak at the posterior side (Arnold and Robertson, 2009). The primitive streak is formed through interaction between BMP, Nodal, FGF and Wnt pathways at the posterior side of the blastocyst and is limited to that side by inhibitory signals originating from the AVE at the anterior side (Arnold and Robertson, 2009; Robertson, 2005). Mechanical pressure, presumably stemming from the uterus wall, was suggested to be necessary for proper AP axis formation, including positioning of the AVE and elongation of the egg cylinder (Hiramatsu et al., 2013), although this dependence was later challenged (Bedzhov et al., 2015). Mechanical strain associated with gastrulation was also demonstrated to drive mesoderm specification in fish, fly and sea anemones through modification of β-catenin (Brunet et al., 2013; Pukhlyakova et al., 2018). Another factor shaping primitive streak and AP axis formation is the visceral endoderm layer, characterized at embryonic day (E)5.5 by the expression of Foxa2. As mechanical signals, localized biochemical signals and neighboring cell layer contexts all co-occur in a specific reproducible pattern during embryogenesis, it is very hard to dissect their respective roles and dependencies in specifying primitive streak location selection and axis formation.

In vitro multicellular models such as embryoid bodies (EBs), and gastruloids derived from EBs, have been used in recent years to dissect certain aspects of early embryonic development (Poh et al., 2014; Turner et al., 2014; van den Brink et al., 2014; Warmflash et al., 2014). These models are formed by 3D aggregation of pluripotent stem cells. Given proper conditions, these aggregates differentiate in a partially organized manner, break their radial symmetry and develop an AP axis (Beccari et al., 2018; ten Berge et al., 2008; van den Brink et al., 2014). Interestingly, when aggregated from a small number of cells and under specific signal conditions, EBs form ‘gastruloids’, which display not only initial symmetry breaking but multiple axes as well as elongation (Turner et al., 2017; van den Brink et al., 2014). Germ layer structure was also shown in confined 2D colonies (Etoc et al., 2016; Morgani et al., 2018; Warmflash et al., 2014). Furthermore, by mixing mouse embryonic stem cells (mESCs) with extra-embryonic cells it is possible to obtain in vitro structures that resemble early stage mouse embryos in terms of both gene expression profile and morphology, including dual axis formation (Harrison et al., 2017; Sozen et al., 2018). The success of these in vitro structures in capturing some of the spatiotemporal aspects of early embryogenesis make them valuable model systems for studying the basic principles of early development, in addition to strategies of in vitro differentiation and tissue engineering. One of the earliest fate decisions that can be observed in embryonic stem cells is their differentiation to mesendoderm progenitors, which are characteristic of the primitive streak and are marked by brachyury (Bra; also known as T) (Gadue et al., 2006; Kubo et al., 2004; Tada et al., 2005). In 2D patterned colonies, the biased onset pattern of Bra has been associated with colony geometry as well as cell density (Blin et al., 2018). In 3D models such as EBs and gastruloids, the polar onset pattern of Bra expression provides one of the first detectable events of symmetry breaking (Beccari et al., 2018; Blin et al., 2018; Boxman et al., 2016; ten Berge et al., 2008; van den Brink et al., 2014). The mechanisms leading to the spatial selection of this pole and axis formation in these 3D models are, however, not clear.

Here we show that the onset of Bra expression is influenced by contact in EBs. Although the entire EB outer shell has the potential to express Bra, expression of this regulator is initiated at the point of EB contact with its surrounding surface. By manipulating the EB environment, we were able to alter the onset of Bra expression to multiple loci, a wide contact surface or an isotropic pattern throughout the outer shell. By creating spatially separate sources of contact and biochemical signals, we could demonstrate how the tissue can integrate these two modalities to define a single Bra locus. Remarkably, Foxa2 appeared with an eccentric bias toward the future onset site of Bra, evident of an earlier symmetry-breaking event and revealing different dependencies on Wnt signaling.

Our findings provide insight into symmetry-breaking mechanisms in 3D tissue models and primitive streak spatial selection in utero, demonstrating how contact and biochemical signals could integrate to drive symmetry breaking and mesendoderm differentiation.

We have previously observed that in EBs differentiated within agarose microwells in the absence of strong Wnt activation, Bra expression predominantly starts from a single locus on the EB external shell and spreads in that layer toward the opposite pole, where a locus is defined as a cluster of Bra-positive cells (see Materials and Methods) (Boxman et al., 2016). However, it is unclear whether the spatial selection of the locus is an outcome of stochastic symmetry breaking within the EB and why Bra expresion starts solely on the outer shell and not in the volume. Mouse E14 Bra-GFP ESCs (Fehling et al., 2003) pre-grown in serum and leukemia inhibitory factor (serum-LIF) were suspended in differentiation medium, aggregated into EBs in hanging drops for 24 h and then seeded in agarose microwell arrays for live imaging up to the onset of Bra-GFP expression (∼72 h), or differentiated within the hanging drops for a similar period (Fig. 1A). In EBs differentiated in agarose microwells, Bra expression commenced from the contact point of EBs with the wall or floor of the microwell (Fig. 1B-F; Movies 1-4). Out of 47 successfully imaged EBs, 16 showed Bra onset at the bottom of the well, 12 showed side onset, 14 showed joint side and bottom onset (Fig. 1C) and five showed two disparate onset loci. We noted that side-onset or dual-onset EBs were generally larger than bottom onset EBs (mean EB radius 120 μm, σ=16 μm, versus mean radius 95 μm, σ=17 μm, P=0.001), raising the hypothesis that side onset is driven by contact with the side wall, which is more common in larger EBs. This was reinforced by examining EBs with two Bra loci, which touched the microwell boundaries at two points or a large surface area overlapping the onset loci (Fig. 1D). We therefore compared in a separate experiment the onset angles in bottom-only contact EBs (n=9) with EBs touching the side wall (n=10) and found a significant difference, indicating that onset location is affected by contact [P=3e–5, two-sample Kolmogorov–Smirnov (KS) test; Fig. 1F]. Interestingly, EBs differentiated in hanging drops also displayed a single locus Bra pattern (Fig. 2B). Because EBs grown in hanging drops lie at the drop's bottom boundary, we hypothesized that their onset is also biased toward this boundary. The alternative result would be random spatial onset, caused by either lack of an organizing event or random rotations of the EB in the drop. Analyzing the locus location in hanging drops (n=12, P=1e–5, KS test), we found that the locus of Bra onset was at the bottom, similar to the bottom-contact EBs grown in microwells (Fig. 1F). Taken together, the results indicate that both mechanical contact and confinement by the liquid-to-air interface can determine the onset location of Bra in EBs. Such contact dependence could be explained by signal consolidation around the contact interface and/or by mechanosensing pathways.

We next wanted to check whether the timing and location of Bra onset can be controlled by imposing contact at a specific locus at an early stage. To this end, we fabricated 200 μm wide rectangular ‘micropool’ structures to enforce early contact at two sides of the EB (Fig. 1G,H). Of the EBs grown in these micropools, 30/41 touched the micropool walls and 11/41 contacted only the well bottom. In the side-contacting group, 29/30 expressed Bra at their side-contact points, compared with 2/11 nontouching EBs with side onset (P<1e–5, Fisher's exact test; Movies 5 and 6). Out of the side-contact, side-onset EBs, 26/29 showed a single locus onset (Fig. 1G) and 3/29 showed simultaneous onset at two loci (Fig. 1H). The preference of one side might result from earlier Bra triggering on that side as a result of faster signal consolidation, allowing enough time for the signal to spread toward the other side and mask or prevent the formation of a second locus. EBs grown in wider micropools (500 μm width) and which did not touch the walls showed bottom Bra expression onset at the contact point with the plate floor (8/10), similar to the case in large microwells (n>50 in multiple experiments; Boxman et al., 2016). We then evaluated the effect of early contact on the timing of Bra onset. When grown in the narrow (200 μm) micropools, EBs created contact with the micropool sides earlier (50–60 h), but Bra onset did not occur earlier than for nontouching controls in the wide micropools (79 h, σ=2 h, n=11 versus 75 h, σ=5 h, n=5). This suggests that cells on the outer shell of the EB did not yet have the potential to express Bra when contact was initially applied to the EB, and more molecular changes had to take place before Bra onset could occur around the contact point.

The observation that EBs can form multiple contact loci or an entire cap contact locus when confined in micropools raised the question of whether the potential to express Bra is limited to specific areas on the outer layer of the EB, to the entire shell, or is uniform throughout the EB volume. Growing EBs in hanging drops, we first confirmed that a canonical Wnt activation reporter (SuTOP-CFP) (Serup et al., 2012) showed colocalization with Bra, in agreement with observations that Bra is activated via the canonical Wnt pathway (Fig. S1A) (Arnold et al., 2000; ten Berge et al., 2008). We then tested EB spatial limitations for Bra onset under external Wnt induction. To this end, we induced the canonical Wnt pathway by adding the small-molecule Wnt pathway agonist CHIR99021 (CHIR) bypassing the pathway receptor via Gsk3β inhibition, adding soluble Wnt3a ligand or embedding Wnt3a-coated beads within the entire EB volume during aggregation (Fig. S1B) (Habib et al., 2013). EBs under CHIR treatment showed isotropic onset with no distinct locus, suggesting that the entire outer shell was Bra competent (n=16/16; Fig. 2B, center). Activation by soluble Wnt3a resulted in a stronger, yet polarized, locus (n=10/10; Fig. 2B, right) (ten Berge et al., 2008). EBs embedded with Wnt3a beads from aggregation displayed isotropic or near isotropic Bra expression onset on the outer shell at 72 h, similar to CHIR-treated EBs, which can be attributed to multiple loci caused by different beads (Fig. S1C). As in our forced early contact results, Wnt3a beads did not lead to Bra expression until 72 h after aggregating the EBs, nor did they induce Bra under pluripotency conditions (Fig. S1B, left). These results confirm that, similar to their susceptibility to contact, outer-layer cells need to reach a certain maturation stage before becoming Wnt-susceptible to activate Bra.

As the entire outer shell is Bra competent and contact is associated with the Bra onset locus, we hypothesized that once given uniform contact conditions, outer layer cells will induce Bra in a spherically isotropic pattern. We therefore tested the Bra onset pattern when EBs were differentiated while embedded in rich ECM composition (Matrigel) or in a more biologically inert hydrogel (low gelling point agarose). When embedded in Matrigel, EBs in a range of sizes displayed isotropic onset patterns (13/15, mean radius at onset 95 μm, σ=50 μm), whereas two large EBs (radii of 120 and 180 μm) initiated Bra from a distinct locus (Fig. 2C; Movies 7 and 8). To test whether the effect of Matrigel was mainly a result of mechanical pressure or of biochemical interactions with its components, we tested the dynamics for EBs embedded in low melting point agarose of lower (0.5%), similar (1%) or higher (1.5%) stiffness than Matrigel. EBs embedded in agarose generated a single locus onset, despite having isotropic contact with the gel (0.5%, n=14/14; 1%, n=10/10; 1.5%, n=11/12; Fig. 2D). Small EBs differentiated in microwells (mean EB radius at onset 78 μm, σ =10 μm, n=9) showed uniform isotropic onset (7/8 EBs, average onset time 78 h, σ=9 h), as also observed in small EBs embedded in Matrigel (Fig. S1D). This suggests that an isotropic pattern can be the result of size, which may be explained by a minimal locus size that is larger than the small EB surface or a minimal signal gradient requirement for a polarized pattern. In summary, a polarized Bra pattern can be obtained at contact points given sufficient EB size, but also under uniform gel contact, suggesting that this is a robust EB behavior. However, as demonstrated by results from CHIR induction, Wnt3a bead embedding and growth in Matrigel, the entire outer shell is competent to activate Bra, given strong activation of the Wnt pathway.

It has been suggested that Oct4 (Pou5f1) and Sox2 expression levels in differentiating ESCs direct their differentiation toward mesendoderm or neural ectoderm, respectively (Thomson et al., 2011). To test whether the competence of the outer shell can be explained by the expression pattern of these early pluripotency factors, we immunostained for Oct4 and Sox2 at 17, 48, 76 and 111 h after the start of differentiation. We observed that expression of both genes receded from the center toward the outer shell with time; by 72 h both Sox2 and Oct4 were highly expressed only at the outer shell of the EB (Fig. 2E).

Which molecular components mediate the effect of contact on Bra onset in outer layer cells? Fibronectin has been shown to induce Bra expression, and a mechanism linking fibronectin to Wnt signaling was proposed (Cheng et al., 2013). We therefore hypothesized that such a link can explain both the relation between contact and Wnt signaling (leading to Bra onset) and the location choice for the Bra locus. We found that fibronectin was expressed isotropically on the outer shell of EBs, implying that its mere expression does not determine the location of Bra onset (Fig. 3A). We then inhibited fibronectin activity using an anti-Fn1 antibody, starting from 0, 24, 48 or 72 h after aggregation under differentiation conditions. Compared with the control group, which displayed normal Bra expression (16/16), fibronectin inhibition added at aggregation or at 24 h inhibited Bra expression (10/10 and 6/6 EBs, respectively, P<2e–5), had little to no effect when added at 72 h (4/4 EBs, P<3e–4) and a mixed effect when added at 48 h (8/8 EBs, P<2e–6) (Fig. 3B). Lack of Bra expression under fibronectin inhibition was also confirmed at the mRNA level (Fig. 3C). Our results suggest that the dependence of Bra-positive fate on fibronectin is limited to early stages, whereas the accumulation of Bra in cells and expansion of its expression to neighboring cells do not depend on fibronectin.

We found that both Wnt signaling and surface contact affect Bra expression onset, although our results suggest that contact bias is mediated through biochemical signaling. To further investigate how these two signal sources relate to each other in this context, we designed a system in which we can create a spatial separation between the surface contact point and a local source of Wnt signaling. To this end, we transfected HEK cells with a plasmid constitutively expressing either Wnt3a-P2A-H2B-CFP, Dkk1-P2A-H2B-Venus or, for the control group, H2B-CFP (Fig. 4A). A clump of HEK cells was then integrated into each EB at a random location shortly after aggregation. The EBs were then differentiated, as in the wild-type experiments, and the location of the Bra locus was analyzed with respect to the contact point and the locus of signaling cells using two-photon microscopy (Fig. 4B-D). In the EBs containing Dkk1 or Wnt3a cell clumps, the Bra locus overlapped with the plane defined by the EB centroid, HEKs and the contact point (n=16; average deviation from co-planarity 9.3°) (Fig. S2A). This allowed us to analyze the joint effects of contact and Wnt3a/Dkk1 signaling on Bra expression locus as angular deviations on a single plane (Fig. 4B). The amount of Wnt3a-secreting agents determined the effect on Bra locus bias from the contact. For small clumps of Wnt3a-expressing cells (3-5 cells), the source of Wnt3a led to a shift in the Bra locus, so that it was positioned between the contact point and the Wnt3a expressing locus (Fig. 4C,E-G). For larger clumps of Wnt3a-secreting cells (>10 cells), the Bra locus was pulled closer to the HEK cluster, maintaining a single locus (P<0.05 for all group pairs, two-sample KS-test; Fig. 4E,F). In one example, where the EB touched the microwell wall at the side opposite to the HEK clump, two Bra loci emerged, one next to the HEKs and one at the side contact point (Fig. S2B). In contrast, a source of Dkk1 shifted the Bra locus away from the contact point (Fig. 4D,E,G; P<0.05 for all group pairs, KS-test). This indicates that Dkk1 can exert its limiting effect over at least 150 μm. Moreover, the locus shift from the contact point, toward or away from the signaling source shows how Wnt3a and Dkk1 signal integration determines the locus location. In control EBs (harboring clumps of H2B-CFP cells) we saw no spatial effect on the locus of Bra onset, which was mostly at the contact point (Fig. 4E,F). These results suggest that contact and Wnt effects on Bra onset are integrated rather than override each other. One possible explanation is that contact mediates Wnt signaling (which is then spatially summed with the other source of Wnt signal). It is also possible, however, that contact acts through a different signaling pathway, which is summed with the Wnt pathway in triggering Bra onset.

To explain the choice of Bra onset location in the absence of external contact surfaces (i.e. in hanging drops), we examined additional genes that are predicted to be expressed at or near the primitive streak before gastrulation. We have previously suggested that once Bra expression begins, the EB progresses at a common developmental pace, delineating a differentiation plan associated with multiple early development genes (Boxman et al., 2016). We therefore hypothesized that the spatial expression pattern of key proteins in EBs before and during Bra onset shows similarity to their in utero pattern around the primitive streak. To test this, we immunostained EBs differentiated in hanging drops for a few early markers at different time points from aggregation, and checked how their spatial expression patterns related to that of Bra. We found that Eomes, a primitive streak marker regulated by Nodal and an EMT driver, was localized to the same pole as Bra although it expanded further on the outer shell (Fig. 5A). E-cadherin (Cdh1) was expressed throughout the EB prior to Bra expression onset, but then its expression declined where Bra was expressed, indicating that epithelial mesenchymal transition takes place in the EB (Fig. 5B) (ten Berge et al., 2008). Loss of E-cadherin enables cell movement, which is supportive of previous observations of gastrulation-like events in EBs and 2D models (Boxman et al., 2016; van den Brink et al., 2014; Warmflash et al., 2014) and is consistent with primitive streak behavior. Oct4, whose expression initially declines in the inner layers, gradually receded on the outer shell toward the opposite pole as Bra expression expanded, similar to the opposing gradients observed in late streak stage epiblasts (Figs 2E and 5C).

These examples demonstrate how the expression and spatial arrangement of different proteins in the EB during pre-streak and through streak phases resemble the equivalent in vivo stages, albeit with an altered topology. Because the primitive streak is initiated adjacent to the visceral endoderm layer in the embryo, we tested for spatiotemporal correlation of the visceral endoderm marker Foxa2 with Bra in the EB. Such a correlation could point to an additional driving factor in Bra locus determination. Foxa2 expression appeared up to 24 h prior to Bra at an internal area within the EBs, with Bra onset following on the shell juxtaposed to Foxa2-expressing cells (Fig. 6A), suggesting a common driver or interlayer dependency. This spatial relation between Foxa2 and Bra was accentuated in large EBs (∼700 μM in diameter at onset) (Fig. 6A, bottom; Fig. S3A), suggesting that these two cell populations have a regulatory relation. To better characterize Foxa2 cells around their onset, we examined the expression of Gata6 and Sox17 at 48 h. Sox17 was not expressed at this early time point. Similarly, no Gata6 expression was detected in the EBs (Fig. S4D), consistent with previous reports (Turner et al., 2017; van den Brink et al., 2014). This indicates that Foxa2-positive cells at this stage do not represent primitive endoderm or visceral endoderm, but probably posterior epiblast or early anterior streak identities (Arnold et al., 2000; Burtscher and Lickert, 2009; Cai et al., 2008; Choi et al., 2012).

To look for candidate genes that could convey the regulation of Bra next to Foxa2-positive cells, we analyzed a published single-cell expression dataset from E6.5 mouse embryos (Scialdone et al., 2016). We found that within epiblast cells, both Foxa2-expressing cells and Bra-expressing cells were significantly enriched for porcupine (Porcn) expression (P=2.6e–8 and P=8.83e–10, respectively). Porcn palmitoylation of Wnt ligands is required for their secretion, suggesting that Foxa2-positive cells in EBs secrete Wnt ligands to the adjacent perimeter, thus promoting the onset of Bra in an adjacent locus. We therefore immunostained for Porcn in differentiating EBs and found that Porcn expression was spatially limited to areas adjacent to and overlapping the Bra onset locus (Fig. 6B). Inhibiting Porcn activity with IWP-2 starting either at aggregation or 48 h after aggregation (when Foxa2 is already expressed in the volume) resulted in the abolishment of both Bra and Foxa2 (Fig. 6C; Fig. S4A). A similar result was obtained by selectively activating non-canonical Wnt with Wnt5a; however, blocking canonical Wnt using XAV-939 or Wnt secretion using IWP-2 resulted in no expression of Bra or Foxa2. Adding just Wnt5a had no significant effect on expression of Bra and Foxa2 compared with the control group (XAV-939+Wnt5a, 0/10, P=1.2e–4; IWP-2+Wnt5a, 0/10, P=1.2e–4; Wnt5a, 8/8; control, 6/6; Fig. S4C). This suggests that both genes require canonical Wnt signaling for their expression (Gadue et al., 2006; Sawada et al., 2005). We were able to rescue Bra by introducing IWP-2 with CHIR to the EBs at 48 h after aggregation; however, Foxa2 was fully downregulated. When Wnt3a was added to IWP-2, we obtained a partial rescue of both Bra and Foxa2. Interestingly, perturbing the EBs with just CHIR or Wnt3a gave a similar result; Bra was over-expressed and Foxa2 was either fully (in the case of CHIR) or partially downregulated and spatially decoupled from Bra (in the case of Wnt3a) (Fig. 6C,D; Fig. S3B). These results suggest that an intricate, probably dynamic, signaling balance involving canonical Wnt determines and maintains the fate of these cells. In small EBs (<65 μm radius) Foxa2 had small to negligible expression whereas Bra encompassed the entire EB (5/5; Fig. S4B). This suggests that the EB requires a certain volume (or length scale) to allow the development of two distinct cell populations, consistent with results in 2D colonies (Morgani et al., 2018; Warmflash et al., 2014). The threshold size might be required to establish sufficient differences, or gradients, in signaling activity. Furthermore, although Bra expression was abolished under fibronectin inhibition starting at aggregation, Foxa2 showed wild-type expression, pointing to a fate decision made at least 24 h before Bra onset and leading to Foxa2-positive cells (Fig. 6C, bottom).

Together, these results show how Bra and Foxa2 are spatially related under canonical Wnt signaling. Bra is monotonically upregulated by canonical Wnt whereas Foxa2 shows a more complex regulation nature, being upregulated by an intermediate level of canonical Wnt activation but downregulated by a strong level (Fig. 6E).

Here we dissect the factors and mechanisms affecting the polarized expression pattern of Bra expression in EBs. These patterns can be formed in EBs and gastruloids with no extra-embryonic tissue and no external Nodal or BMP sources (Boxman et al., 2016; Turner et al., 2017; van den Brink et al., 2014). We have shown that physical contact with surfaces, a localized source of Wnt signaling and proximity to Foxa2-expressing cells can all induce or predict the onset of Bra. Interestingly, in the embryo, these three conditions co-occur near the future location of the primitive streak, making it hard to dissect their respective roles. Bra expression begins during gastrulation near a source of Wnt signaling (arising from a positive feedback loop between Bmp4 sourced from the ExE and Nodal and Wnt at the posterior epiblast), adjacent to Foxa2-positive cells, in a region that might have experienced high mechanical pressure as a result of implantation. We analyzed some of the hierarchy between these three primitive streak predictors (Fig. 6E).

Bra onset polarization was robust to different external conditions and EB parameters, suggesting that it is the ‘default’ behavior in this system. A similar result was demonstrated in mouse embryos grown in vitro, where mechanical contact might direct AP axis formation and positioning of the AVE, but they can also robustly form in hanging drops (Bedzhov et al., 2015; Hiramatsu et al., 2013). The polar pattern in EBs can be obtained even in the absence of any chemical or physical cue directing its orientation (as seen in the agarose embedding experiments). Overriding this polarization process pattern in EBs required specific intervention, such as continuous canonical Wnt activation by CHIR or application of inducing sources (contact points or localized Wnt sources) at two distant locations on the EB surface (Fig. 1G; Fig. S2B). Contact with external surfaces, or with the hanging drop boundary surface, was sufficient to bias the location of Bra onset and orient the EB primary axis. We have shown how we can generate two Bra loci through controlled contact, depending on the geometry of the surrounding surfaces. The fact that the onset of Bra expression in hanging drops or large wells was similar to that in narrow micropools suggests that it is not the effect of mechanical pressure, but of the interface with contact surface or of geometrical confinement. These could act through limiting the diffusion of secreted signal molecules, such as Wnt ligands, resulting in locally higher signal concentration (Siggia and Warmflash, 2018). Another option is that the ligand activation mechanism is bypassed through interfacing with outside surfaces by integrin-mediated signaling (Cheng et al., 2013) or by mechanotransduced phosphorylation of β-catenin (Brunet et al., 2013), although the localization of the Bra locus in hanging drops to the bottom air–liquid interface suggests that diffusion confinement is sufficient. Interestingly, the factors that we found to affect the Bra locus in 3D were very different from those in the case of 2D patterned colonies, where a recent study related the percentage of Bra-positive cells to cell density, and their localization to high curvature regions at the tip of elongated colonies (Blin et al., 2018). In 2D colonies, there is no equivalent for limited ligand diffusion in 3D near surfaces, and alternative factors might bias the response potential in different locations within the colony (Etoc et al., 2016).

We found that the potential to initiate Bra expression was limited to cells at the outer layers of the EB, at least in large-enough EBs, and its expression did not expand much inwards. Whole-volume Bra expression was observed mostly in small EBs. The patterns of Sox2 and Oct4 expression suggest there could be multiple expression and epigenetic differences between different EB layers prior to Bra onset, which could explain the biased expansion pattern, despite positive feedback between Bra and Wnt signaling (Evans et al., 2012; Yamaguchi et al., 1999). The limitation of Foxa2 initiation to inner EB cells provides further support for the differential potentials in different EB layers.

We found that we can ‘pull’ or ‘push’ the Bra locus away from the contact point by using an external source of Wnt3a- or Dkk1-secreting cells, respectively. This demonstrates that contact and biochemical signal modalities can be summed by the EB, resulting in a single locus, as opposed to the two loci or one locus overriding the other that we observed in the forced contact experiments. This can be attributed to the gradient strength or to the surface distance between the loci, which might be too large for signal integration to have an effect. One possible mechanism for this integration is activation of Wnt signaling through mechanical contact. The dependence of Bra activation on fibronectin points to such a mechanism (Cheng et al., 2013). Interestingly, the mechanism proposed by Cheng et al. (2013) bypasses the Wnt pathway receptors and therefore cannot explain why Dkk1 diverts Bra onset away from the contact point, suggesting other mechanisms are at play. Another option is that Wnt ligands consolidate at the contact surfaces, or have limited diffusion range near them, leading to locally higher concentration (Siggia and Warmflash, 2018). This option is particularly reinforced by localization of the Bra locus to the air–water interface in hanging drops, where we do not expect to find integrin interactions. The isotropic peripheral onset we observed in small EBs could be attributed to a lack of Wnt gradient formation.

The relation of Bra onset to Foxa2-positive cells provides some interesting insights. First, from the differential response to fibronectin inhibition, we conclude that the decision junction of differentiating ESCs within the EB toward Bra and Foxa2 tissue layers takes place at 24-48 h prior to Bra onset. The lack of overlap with visceral or definitive endoderm markers suggests that Foxa2-positive cells initially represent a non-endodermal population, such as anterior primitive streak (Burtscher and Lickert, 2009; Gadue et al., 2006). Second, the observation that small EBs do not express Foxa2 suggests that a signal gradient controls segregation between these two cell populations, which potentially represent posterior and anterior primitive streak. This is reinforced by the dependency of Foxa2 expression on intermediate canonical Wnt signaling, as well as by the high spatial dependency seen between these two cell populations in large EBs (Fig. 6; Fig. S3). The facts that different perturbations can lead to Bra or Foxa2 expression while abolishing the other fate in the EB and that the spatial adjacency of these two cell populations can be decoupled using external Wnt signaling suggest that the observed spatial relation is a result of the dependence of both proteins on different levels of canonical Wnt signaling, which can be spatially coupled through a morphogen gradient. The earlier onset of Foxa2 relative to Bra might be because Foxa2-competent cells, which form within the inner EB layers, require lower levels of canonical Wnt (which can be reached earlier), whereas Bra requires higher levels of canonical Wnt. Another option is that the Foxa2 locus (or its relevant enhancer) becomes accessible earlier than that of Bra. Further molecular work is needed to decipher the Foxa2 onset mechanisms.

In contrast to its location, the timing of Bra expression onset in EBs was not affected by external or localized Wnt3a signaling, or by inducing early contact in narrow microwells. These results point to the maturation needed in outer shell cells before they become Wnt- and/or contact-responsive to induce Bra. Further characterization of the chromatin or expression changes defining that ‘mature’ state are needed.

Through different perturbations on a 3D in vitro system we were able to delineate the basic dependencies of an early fate decision on contact and biochemical factors that co-occur in the embryo. The insights we gained here improve our understanding on the forces shaping differentiation within in vitro systems such as EBs and gastruloids and shed more light on how they are integrated in vivo.

E14 Bra-GFP mESCs (kindly provided by Dr Gordon Keller, McEwen Stem Cell Institute, UHN, Toronto, Canada), R1 mESCs and SuTOP-CFP,AR8-RFP mESCs (kindly provided by Dr Palle Serup, University of Copenhagen, Denmark) were cultured using standard conditions on irradiated primary mouse embryonic fibroblasts and knockout DMEM containing 15% fetal bovine serum, 50 μg/ml penicillin/streptomycin, 2 mM L-glutamine, 100 μM nonessential amino acids, 80 μM β-mercaptoethanol and 10³ U/ml LIF.

Embryoid bodies were formed in hanging drops under differentiation conditions. Each drop contained 25 μl differentiation medium (IMEM containing 20% FBS, 50 μg/ml penicillin/streptomycin, 2 mM L-glutamine, 100 μM nonessential amino acids, 100 μM β-mercaptoethanol) with approximately 300 mESCs; 100 cells for small EB experiments, 2000 cells for large EB experiments. The number of cells in hanging drops was controlled using a cell counting slide and dilution of the cell-containing medium to obtain the desired starting cell number per single EB within each 25 μl hanging drop. Cells were aggregated for 24 or 48 h before being transferred to micropools or microwells for imaging in the same medium. Alternatively, the EBs were kept in hanging drops for the entire duration of the experiment before being imaged or fixed and immunostained. In some of the staining experiments (Figs 3 and 6A; Fig. S3), R1-derived EBs were used. In these EBs, Bra onset typically occurs about 12 h later than in E14-derived EBs. When adding inhibitors or activators to hanging drops, we aggregated the EBs in 20 μl hanging drops and then added at the appropriate time 5 μl of medium with the required substance. For signal perturbation experiments, differentiation medium was supplemented with CHIR99021 (5.35 μM) or recombinant Wnt3a (500 ng/ml). For fibronectin inhibition, the medium was supplemented with anti-Fn1 antibody (Abcam, ab6328). To incorporate Wnt3a-bound beads within the EB, prepared beads (kindly provided by Shukry J. Habib, King's College London, UK) were premixed with mESCs before distributing the suspension to 25 μl hanging drops for aggregation. We added 12,000 beads to 10 ml of medium with the mESCs, averaging 30 beads per EB. Bead integration into the EBs was confirmed by two-photon microscopy (Fig. S1B).

Low melting point agarose (Sigma-Aldrich, A9414) was mixed with preheated double-distilled water and vortexed to avoid clumps. It was then deposited on the glass bottom of a 24-well plate to allow two-photon imaging. EBs were injected into the gel when it had cooled down sufficiently, about 24 h after aggregation in hanging drops. Matrigel (Gibco Geltrex, A1413202) was thawed on ice in 4°C overnight. Matrigel drops (20 μl each) were deposited on the glass bottom of a 24-well plate. EBs were transferred into the drops at 48 h after aggregation; after gel solidification, the wells were filled with differentiation medium.

HEK 293 cells were split to three wells (12-well plate, 50% confluence). After 12 h, the cells were transfected with constitutive pCMV-H2B-Cerulean, pCMV-Wnt3A-P2A-H2B-Cerulean or pCMV-DKK1-P2A-H2BVenus plasmids, using the Xfect transfection reagent (631318). At 24 h post transfection, the transfected HEKs were trypsinized, suspended as single cells and injected into hanging drops containing 12-h-old E14 Bra-GFP EBs at different concentrations (10, 20 or 50 HEKs per EB). After 24 h, the injected EBs were transferred into microwells and imaged using a two-photon microscope.

To create the silicon template for the shaped micropools and microwells, we implemented a high aspect ratio etching system usually used for microelectronics fabrication. Using deep reactive-ion etching (DRIE), we were able to determine the micropattern's exact depth by selecting a wafer with a specific handle/device thickness. This suggests an advantage over the common photoresist techniques, which require fine-tuning of the process to obtain an exact pattern depth. We fabricated 200-500 μm diameter microwells and 100-500 μm wide micropools, with a depth of approximately 400 μm. The patterns were then coated with PTFE to minimize friction with the template when pulling the elastomer out. To create the microwell/pool negative, PDMS elastomer (Sylgard 184) mixed with a curing agent (1:10) was poured into the silicon wafer mold and incubated overnight at 50°C. Individual stamps were then cut out using a scalpel. To generate the positive agar microwells/pools we placed the stamps face down inside glass-bottom plates coated with poly-D-lysine and injected agarose under the stamps. The plates with the injected stamps were then placed in a vacuum chamber for 90 min, followed by another round of agar injection and an additional 10 min in the vacuum chamber. The stamps were then removed, leaving the micropatterns imprinted in the agar.

EBs were aggregated in hanging drops on plastics slides. The slides were mounted in the two-photon microscope and imaged with EC Plan-Neofluar 20×/0.50 dry objective, providing a working distance of 2 mm. The data were then segmented in Imaris and analyzed in MATLAB.

We sourced data from both 2D and 3D live imaging experiments. 2D epifluorescence experiments offer a short acquisition cycle for a large cohort of EBs, thus yielding accurate temporal tracking and statistics based on a large number of samples. On the other hand, they do not possess elevation distinction and the superior resolution and cell separation that more advanced microscopy techniques offer. For 3D imaging, we used a two-photon microscope with which we were able to have a full view of signal dynamics and track individual cells. However, this technique requires large data analysis, has a longer acquisition cycle and is less compatible with multiple EB imaging. Two-photon imaging was carried out using a Zeiss LSM7 inverted two-photon microscope with a Plan-Apochromat 20×/0.8 dry objective. Each EB was scanned at 3 μm intervals along the z-direction. Horizontal resolution was set to 512×512 pixels at approximately 0.6 μm per pixel, 8 bit depth. GFP was excited at 930 nm. RFP, Alexa Fluor 488, 405 and 594 were excited at 800 nm. CFP was excited at 860 nm. Epifluorescence imaging was carried out using a Nikon TiE epifluorescence microscope equipped with a motorized XY stage (Prior) at 10× magnification using NIS Elements software. Acquisitions were taken every 25 or 30 min for 14–72 h. For live imaging, in both systems we used Okolab incubation cages, maintaining 5% CO₂ and 37°C.

Three different primer pairs were designed for different Bra exons, and two for each of the housekeeping genes Ywhaz and Rock2. EBs were gathered from different pools, distinct from each other by experiment group. mRNA was extracted by disassembling by pipetting and incubating the pooled EBs for 5 min in 0.5 ml Trizol (Sigma-Aldrich T9424-25 ml), followed by adding 100 µl chloroform and another 5 min incubation. After centrifuging the samples for 10 min at 12,000 g at 2-8°C, the translucent phase was incubated in IPA 100% for 10 min at room temperature. Centrifuging the samples again in the same manner gave a pellet containing the mRNA, to which we added 100 μl 70% ethanol, vortexed the samples and then centrifuged the samples for 5 min at 10,000 g at 2-8°C. The ethanol was then removed and the rest vaporized within 15 min, followed by addition of double-distilled water. We then produced cDNA (Invitrogen, 18080-051) and performed RT-PCR (Thermo Fisher StepOnePlus Real-Time PCR system) using Fast SYBR Green (Thermo Fisher, 4385612). Each Bra primer pair was run in duplicate. Results were normalized by subtracting the average result for housekeeping genes (ΔCt) and then subtracting the control group (−ΔΔCt).

We fixed the embryoid bodies with 4% paraformaldehyde overnight at 4°C, followed by a rinse cycle and an overnight wash in PBS. We added primary antibodies against Eomes (Abcam, ab23345), Bra (R&D Systems, AF2085), Foxa2 (Abcam, ab108422-100), E-cadherin (Cell Signaling, 24E10), Oct4 (Santa Cruz Biotechnology, sc-8628) and Sox2 (Santa Cruz Biotechnology, sc-365964) with blocking buffer (PBS, 0.1% Triton X-100, 5% FBS) and incubated overnight at 4°C, followed by four wash cycles in blocking buffer for at least 30 min per cycle. We then added Alexa Fluor-conjugated secondary antibody in blocking buffer and incubated it overnight, followed by four wash cycles as described for the primary antibody. Antibody concentrations used were the same as recommended by the manufacturer for 2D staining.

The two-photon microscopy data were contrast enhanced in Fiji and then spot segmented in Imaris, where the integration radius was set to 7 μm and the filter was set to sum intensity with a threshold of 110/255. The background was subtracted prior to segmentation. The segmented data were read and analyzed in MATLAB. For Bra locus estimation in hanging drops or microwell experiments, locus elevation and azimuth angles were calculated from the segmented two-photon data by the sum vector of all Bra-positive cells. Bra onset was defined as ‘polarized’ when Bra-GFP occupied less than 25% of the shell, in one locus. The EB center and radius were estimated by Bra-positive cell curvature. The angle between the EB bottom point and Bra locus was then computed. For HEK incorporation experiments, we developed a MATLAB script that calculates geometrical properties between the HEK cluster, Bra locus and contact point. The script detects HEK cells at each time point by channel and intensity. We calculated the geometrical mean of the HEKs and set it as the HEK location point. Next we looked for the Bra locus, where a locus is defined as a minimum of five cells above the intensity threshold of 100 within a distance of 50 μm from each other, all in the green channel. Next, we defined the contact point of the EB with the well bottom as the bottom point of the sphere and calculated the geometric properties between all pairwise combinations of HEK cluster, Bra locus and contact point. The geometrical properties calculated were the spatial angles, Euclidean distances and geodesic distances. Side onset in microwells was manually detected by superimposing the GFP signal on brightfield images in which both the EB and the microwell wall are clearly visible, as shown in Fig. 1. Angle quantification was automatically performed using a MATLAB script on 3D segmented data.

We thank Shukri Habib for providing Wnt3a-coated beads, Palle Serup for the SuTOP-CFP mESC line, Amitai Menuchin for help with RT-PCR, and David Sprinzak, Omri Wurzel and the three anonymous reviewers for providing helpful feedback on the manuscript.