Dissecting the dynamics of signaling events in the BMP, WNT, and NODAL cascade during self-organized fate patterning in human gastruloids

Abstract

During gastrulation, the pluripotent epiblast self-organizes into the 3 germ layers-endoderm, mesoderm and ectoderm, which eventually form the entire embryo. Decades of research in the mouse embryo have revealed that a signaling cascade involving the Bone Morphogenic Protein (BMP), WNT, and NODAL pathways is necessary for gastrulation. In vivo, WNT and NODAL ligands are expressed near the site of gastrulation in the posterior of the embryo, and knockout of these ligands leads to a failure to gastrulate. These data have led to the prevailing view that a signaling gradient in WNT and NODAL underlies patterning during gastrulation; however, the activities of these pathways in space and time have never been directly observed. In this study, we quantify BMP, WNT, and NODAL signaling dynamics in an in vitro model of human gastrulation. Our data suggest that BMP signaling initiates waves of WNT and NODAL signaling activity that move toward the colony center at a constant rate. Using a simple mathematical model, we show that this wave-like behavior is inconsistent with a reaction-diffusion-based Turing system, indicating that there is no stable signaling gradient of WNT/NODAL. Instead, the final signaling state is homogeneous, and spatial differences arise only from boundary effects. We further show that the durations of WNT and NODAL signaling control mesoderm differentiation, while the duration of BMP signaling controls differentiation of CDX2-positive extra-embryonic cells. The identity of these extra-embryonic cells has been controversial, and we use RNA sequencing (RNA-seq) to obtain their transcriptomes and show that they closely resemble human trophoblast cells in vivo. The domain of BMP signaling is identical to the domain of differentiation of these trophoblast-like cells; however, neither WNT nor NODAL forms a spatial pattern that maps directly to the mesodermal region, suggesting that mesoderm differentiation is controlled dynamically by the combinatorial effect of multiple signals. We synthesize our data into a mathematical model that accurately recapitulates signaling dynamics and predicts cell fate patterning upon chemical and physical perturbations. Taken together, our study shows that the dynamics of signaling events in the BMP, WNT, and NODAL cascade in the absence of a stable signaling gradient control fate patterning of human gastruloids.

Fig 1. WNT initiates and NODAL up-regulates mesodermal differentiation.

(A–C) Images of samples immunostained for the indicated markers at 44 h post BMP treatment in different conditions: control and NODAL −/− cells were treated with 50 ng/ml BMP4. BMP+IWP2 represents wild-type cells treated with 50 ng/ml BMP and 5 μM IWP2. Quantification represents intensity levels of indicated markers normalized to DAPI, averaged at different positions along the colony radii. Error bars represent standard error of the mean. N ≥ 10. Scale bars = 100 μm. Colony diameter = 700 μm. Underlying data can be found in S1 Data. Data structure explanation can be found in S1 Text. BMP, Bone Morphogenic Protein.

Fig 2. BMP controls the differentiation of CDX2+ edge cells, which transcriptionally resemble in vivo human trophoblast cells.

(A) Images of samples immunostained for the indicated markers at 44 h post BMP treatment in different conditions—BMP4 alone, BMP and IWP2, and BMP4 and SB. Quantification represents intensity levels of indicated markers normalized to DAPI, averaged at different positions along the colony radii. Error bars represent standard error of the mean. N ≥ 10. (B) Images of cells grown in standard conditions in sparse culture and immunostained for the indicated markers at 48 h post BMP4 treatment in different conditions. No BMP4 was added to the mTeSR sample. Quantification represents average mean intensity levels per cell of indicated markers normalized to DAPI. Cells N > 500. Error bars represent standard deviation across cells. (C) Number of differentially expressed genes (p ≤ 0.05, fold change ≥ 2, FPKM ≥ 1 in at least one sample of the 4 samples) between indicated samples. (D) (Left) Heatmap showing z-scores of 284 most differentially expressed genes as defined in the text. (Right) Venn diagram of differentially down-regulated/up-regulated genes between indicated samples as compared to mTeSR1. Red: BMP4, Green: BMP4 + IWP2, Blue: BMP4 + SB. (E) Pearson correlation coefficients for normalized read count values (FPKM) of 284 differentially expressed genes between indicated samples. (F) (Left) Pearson correlation coefficients between indicated samples for z-scores of 284 differentially expressed genes determined from hESC samples. (Right). Pearson correlation coefficients between indicated samples for z-scores of 174 lineage-specific genes determined from in vivo human embryo data. (G) Normalized read count values, expressed in FPKM, of common trophoblast markers in the indicated samples. (H) Images of samples immunostained for TFAP2A in micropatterns, and regular culture. Quantification as defined in (A). N ≥ 10 colonies for micropatterns, N > 500 cells for regular culture. Scale bar = 100 μm. Underlying data can be found in S2 Data. BMP, Bone Morphogenic Protein; E5/E6/E7, embryonic day 5/6/7 (respectively); EPI, epiblast; FPKM, fragments per kilobase of transcript per million mapped reads; hESC, human embryonic stem cell; PE, primitive endoderm; SB, SB431542; TE, trophectoderm.

Fig 3. WNT signaling dynamics lie outside the Turing instability regime.

(A) Schematic and equations of the model. (B) Evolution of activator levels from the initial state to the steady state. Inequalities define parameter regimes that display each of the 2 behaviors (S1 Model). In each case, the top row shows the model simulated on a circular colony, while the bottom shows it simulated in a domain with cells throughout. Simulation parameters: s_A = 0.01, k_I = 1, k_A = 0, kd_A = 0.001, D_I = 0.4, s_I = 0.01, kd_I = 0.008. (left) D_A = 0.014, (right) D_A = 0.0025 for outside and inside Turing instability regime, respectively. Simulation domain: 190 × 190 pixel square lattice with periodic boundary conditions and random distribution of activator/inhibitor as initial conditions. To simulate the model in circular colonies, a circle (radius: 25 pixels) is defined at the center of the lattice. Outside the circle, s_A = s_I = 0, and kd = 0.01. (C) Snapshots of GFP-β-catenin hESCs from time-lapse imaging at indicated time points post BMP treatment. Marked regions are magnified in images shown below. (D) Average nonmembrane β-catenin intensity levels as a function of distance from the colony edge. The legend indicates the time post BMP treatment represented by each curve. Green and blue dots represent the front and back edges of the active signaling domain, respectively. (E) The positions of the inner and outer edge of the signaling domain as a function of time post BMP addition. Magenta line shows linear fit with equation 6.02 μm/h × t + 13.57 μm, R² = 0.98. (F) Kymograph showing spatiotemporal evolution of nonmembrane β-catenin levels. At time points earlier than the first time plotted in D, signaling is below threshold signaling at all positions in the colony. Scale bar = 100 μm. Underlying data can be found in S3 Data. BMP; Bone Morphogenic Protein; GFP; Green Fluorescent Protein; hESC, human embryonic stem cell.

Fig 4. Inward movement of signaling activities is not due to cell movement.

(A) Snapshot from time-lapse imaging at 47 h post BMP treatment in the bright field channel. Green and red dots represent initial and final positions, respectively, of labeled cells that were correctly tracked throughout imaging. (B) Histogram of displacement of tracked cells. n = 84 cells. (C) Image of colony (shown in A) immunostained for CDX2, BRA post live cell imaging. Quantification represents average nuclear intensities of indicated markers normalized to DAPI as a function of radial position. Error bars represent standard error of the mean. N = 6. The lines indicate the inner and outer positions of a region of high BRA as a rough guide to the different fate territories (D) Radial displacement of cells as a function of their starting position. Radial displacement = distance of the cell from the center at end of imaging − distance of the cell from the center at the start of imaging. The vertical dotted lines are in the same positions as in (C). The horizontal dashed line separates the cells moving toward the edge from cells moving toward the center. Underlying data can be found in S4 Data. BMP, Bone Morphogenic Protein; BRA, BRACHYURY.

Fig 5. Regulation of WNT signaling activity by BMP, WNT, and NODAL.

(A, E) Snapshots of GFP-β-catenin hESCs from time-lapse imaging movies at 46 h post BMP treatment. Time between BMP4 and LDN or IWP2 addition is indicated above each image. Boxes show the region represented at higher magnification below each colony image. (B, F) Average nonmembrane β-catenin levels as a function of radial position. Control represents samples with no inhibitor addition. The timing of LDN or IWP2 addition after BMP4 treatment for each curve is shown in the legend, while the time being analyzed in shown above the plot. Error bars represent standard error. N = 3 colonies for the samples LDN at 0 h (LDN@0h), LDN at 11 h (LDN@11h); for all other samples N ≥ 5. (C, G) Temporal evolution of the position and intensity of peak signaling (defined by maximal nonmembrane β-catenin intensity). (D, I) Temporal evolution of the front of the domain of active signaling. In (C, D, G, H, and I) at time points earlier than the first one in each curve, signaling was below threshold signaling at all positions. (H) Temporal evolution of average nonmembrane β-catenin levels in a narrow ring inside the region of mesodermal differentiation (distance from edge: 134–150 μm) and at the center (distance from edge: 277–350 μm) of the colony. (J) Images of GFP-β-catenin hESCs at 46 h post treatment with either BMP4 (control) or BMP4 and SB. Average nonmembrane β-catenin levels as a function of radial position at 46 h. Error bars represent standard error. N ≥ 6. Underlying data can be found in S5 Data. BMP, Bone Morphogenic Protein; GFP, Green Fluorescent Protein; hESC, human embryonic stem cell; LDN, LDN193189.

Fig 6. NODAL signaling moves inward independently of BMP and WNT.

Images of samples immunostained for SMAD2/3 after 44 h of BMP treatment. The time between BMP4 and inhibitor (LDN or IWP2) addition is indicated above the image. No chemical inhibitor was added in the No-LDN sample. Quantification represents average nuclear intensities of indicated markers normalized to DAPI as a function of radial position. Error bars represent standard error. N ≥ 10. Scale bar = 100 μm. Underlying data can be found in S6 Data. BMP, Bone Morphogenic Protein; LDN, LDN193189.

Fig 7. WNT enhances and NODAL inhibits BMP signaling at colony edge.

Images of samples immunostained for the indicated markers after 44 h of treatment in different conditions: control (wild-type) hESCs and NODAL−/− cells were treated with 50 ng/ml BMP4. The sample labelled as IWP2 at 0 h represents wild-type cells treated with 50 ng/ml BMP and 5 μM IWP2. Quantification represents intensity levels of indicated markers normalized to DAPI, averaged at different positions along the colony radii. Error bars represent standard error of the mean. N ≥ 10. Scale bar = 100 μm. Underlying data can be found in S7 Data. BMP, Bone Morphogenic Protein; hESC, human embryonic stem cell.

Fig 8. Duration of BMP signaling controls extra-embryonic differentiation at colony edge.

Images of samples immunostained for fate markers after 44 h of BMP treatment. The time between BMP4 and LDN addition is indicated above the image. Quantification represents average nuclear intensities of indicated markers normalized to DAPI as a function of radial position. Error bars represent standard error. N ≥ 10. Scale bar = 100 μm. Underlying data can be found in S8 Data. BMP, Bone Morphogenic Protein.

Fig 9. Continuous WNT and NODAL signaling synergize to achieve maximal mesodermal differentiation.

(A, B, D) Images of samples immunostained for fate markers after 44 h of BMP treatment. The time between BMP4 and chemical inhibitor (LDN/IWP2/SB) addition is indicated above the image. Quantification represents average nuclear intensities of indicated markers normalized to DAPI as a function of radial position. Error bars represent standard error. N ≥ 10. (C) (Left image) Image of a colony immunostained for BRA and DAPI after 47 h time-lapse imaging of BMP-treated GFP-β-catenin hESCs. (Right image) Snapshot from time-lapse imaging for the same colony at 47 h. (Left plot) Quantification represents average intensity levels for BRA and nonmembrane β-catenin at 47 h, normalized by DAPI, plotted as a function of radial position. Error bars represent standard error of the mean. N = 9. (Right plot) Average BRA intensity levels as a function of β-catenin levels, color-coded by distance from the colony edge (μm). Scale bar: 100 μm. Underlying data can be found in S9 Data. BMP4, Bone Morphogenic Protein; BRA, BRACHYURY; GFP, Green Fluorescent Protein; hESC, human embryonic stem cell.

Fig 10. A simple mathematical model recapitulates signaling dynamics and predicts cell fate patterning.

(A) (Left) Model equations. u and v represent WNT and NODAL signaling, respectively (parameter descriptions and values in S2 Model). (Right) Time evolution of WNT and NODAL to steady state. (B) The position of inner and outer edge of simulated WNT signaling levels as a function of time. (C) WNT levels at steady state (46 h) as a function of edge distance for different simulations. LDN addition was simulated by setting WNT’s dependence on BMP to 0 at indicated time points. IWP2 addition was simulated by setting autoactivation of WNT to 0 at indicated time points. (D) NODAL levels at steady state (46 h) for different simulations. (E, F) Fates assigned by the model on colonies of different sizes (E), shapes (F). (G) hESCs immunostained for CDX2, BRA, SOX2 after 44 h of BMP treatment on triangular and pacman colonies. Immunostaining data from n = 18 colonies were used to calculate average fate territory maps shown adjacent to image. Scale bar = 100 μm. Underlying data can be found in S10 Data. BMP, Bone Morphogenic Protein; BRA, BRACHYURY; hESC, human embryonic stem cell; LDN, LDN193189.