Quantitative analysis of transcriptome dynamics provides novel insights into developmental state transitions

Abstract

Background

During embryogenesis, the developmental potential of initially pluripotent cells becomes progressively restricted as they transit to lineage restricted states. The pluripotent cells of Xenopus blastula-stage embryos are an ideal system in which to study cell state transitions during developmental decision-making, as gene expression dynamics can be followed at high temporal resolution.

Results

Here we use transcriptomics to interrogate the process by which pluripotent cells transit to four different lineage-restricted states: neural progenitors, epidermis, endoderm and ventral mesoderm, providing quantitative insights into the dynamics of Waddington’s landscape. Our findings provide novel insights into why the neural progenitor state is the default lineage state for pluripotent cells and uncover novel components of lineage-specific gene regulation. These data reveal an unexpected overlap in the transcriptional responses to BMP4/7 and Activin signaling and provide mechanistic insight into how the timing of signaling inputs such as BMP are temporally controlled to ensure correct lineage decisions.

Conclusions

Together these analyses provide quantitative insights into the logic and dynamics of developmental decision making in early embryos. They also provide valuable lineage-specific time series data following the acquisition of specific lineage states during development.

Introduction

How a single cell ultimately gives rise to a patterned, complex organism is a fundamental question in biology. Embryonic development can be generalized as a process of progressive restriction of cellular potential. In vertebrates, the zygote is totipotent, but by blastula stages the three primary germ layers, ectoderm, mesoderm and endoderm, have been specified. The fates of cells within these germ layers then become progressively restricted to single differentiated cell types characteristic of that germ layer. Conrad Waddington famously depicted this process as a topological landscape [1]. In his model, a ball positioned at the top of the landscape represents a cell with all developmental pathways open to it. As the ball progresses down the landscape, the paths it takes dictate which lineage states will remain accessible. Waddington noted that the valleys or channels of the landscape arise from the interactions between genes and from their interactions with the cell’s environment.

At blastula stages vertebrate embryos possess a transient population of pluripotent cells which, like the fertilized egg, occupy a position atop Waddington’s landscape. These cells- inner cell mass cells in mammals and naïve animal pole cells in amphibians- can give rise to the derivative cell types of all three germ layers and as such can recapitulate the path to different lineage states including the relevant gene regulatory network (GRN) topology and dynamics. Studies using explants of pluripotent cells from Xenopus blastulae (so called “animal caps”) have been central to our current understanding of the signals and transcriptional responses that direct these stem cells toward specific lineage states [2,3,4,5,6]. Some of these signals emanate from the blastopore lip, or the Spemann-Mangold organizer, and help to direct formation and patterning of the primary germ layers [7,8,9]. Exit from the pluripotent state also coincides with the loss of expression of many maternally provided pluripotency transcripts, such as Pou5f3.3 and Foxi2 [10,11,12,13]. As cells exit pluripotency, their potential becomes progressively restricted until their fate becomes specified and then determined.

Animal pole cells in Xenopus are fated to give rise to ectodermal derivatives. In situ these cells give rise to both epidermal and neural progenitor cells, as well as neural crest and cranial placodes, under the direction of signals from the organizer. Absent these signals, animal pole cells are directed by endogenous BMP signaling to become epidermis [14]. BMP2, 4, and 7 have all been shown to be potent epidermal inducers [15, 16] and BMP4/7 heterodimers have been identified as the most physiologically relevant ligands in early embryos [17]. The binding of these ligands to type I and II BMP receptors results in phosphorylation of Smad1/5/8 and the translocation of these phospho-Smads to the nucleus together with Smad4, resulting in transcription of target genes including Epidermal Keratin (EpK) and Dlx3 [18,19,20,21,22,23].

Although isolated animal pole cells will transit to an epidermal state absent additional instructions, it has been proposed that the default state of these initially pluripotent cells is a neural progenitor state [24,25,26]. This model arose from the findings that neural “inducing” factors secreted by the organizer, including Noggin, Chordin, Follistatin and Cerberus, function as BMP antagonists [27,28,29,30,31,32], as well as from the observation that cells neuralize upon dissociation [33]. The “neural default” model has not been without controversy, however. In chick embryos the expression patterns of BMPs and their antagonists do not fully fit the neural default model, misexpression of BMP antagonists does not induce neural progenitors, and ectopic expression of BMP fails to inhibit neural plate formation [34,35,36,37]. There is also evidence in Xenopus that FGF signaling may be required for competence to respond to BMP antagonists [38,39,40], although the roles of FGF and BMP signaling in neural induction have shown to be separable [41]. By contrast, Xenopus animal pole explants exposed to BMP antagonists such as Noggin adopt a neural progenitor state and express neural genes including Sox2/3 and Otx1/2. These transcription factors, and their homologs, play an important role in specification of the central nervous system (CNS) in vertebrates [42,43,44,45,46,47,48].

Explants of Xenopus pluripotent blastula cells have also played a central role in determining the signals that control formation of the other embryonic germ layers, mesoderm and endoderm [49,50,51,52]. Nodal, Activin, and Vg-1, ligands of the other branch of the TGF-beta family, have been shown to act as morphogens, with low levels of signaling inducing mesoderm and high levels inducing endoderm [52,53,54,55,56,57,58,59,60,61]. Treatment of these cells with exogenous Activin mimics the activity of Nodal/Vg-1 [62] promoting phosphorylation of receptor Smads 2/3 [63, 64]. High levels of Activin/Nodal signaling induce endoderm as evidenced by expression of key factors Sox17 and Endodermin [65, 66] whereas lower doses induce both dorsal and ventral mesoderm [67, 68]. BMP4/7 heterodimers also possess mesoderm inducing activity, but it is limited to ventral, not dorsal mesoderm [69,70,71]. Much still remains to be learned about the signaling dynamics and outputs of BMP4/7 versus Activin-mediated regulation of cell state transitions during germ layer formation.

Xenopus animal pole explants are an ideal system to probe the dynamics of developmental decision-making during lineage restriction. These cells can be cultured in a simple salt solution and will undergo lineage restriction on the same time scale as if they had remained in vivo. Given appropriate signals, animal pole cells can be directed to any lineage progenitor state within a time frame of approximately seven hours. Analogous experiments in cultured mouse or human embryonic stem cells (ESCs) take more than a week in culture, and unlike ESCs Xenopus cells do not need to be artificially retained in a pluripotent state that may be distinct from the transient pluripotency that exists in vivo. The unique features of the Xenopus system allow initially pluripotent cells to be followed at high temporal resolution as they progress to lineage restriction, providing insights into the dynamics of developmental decision making in early vertebrate embryos. Here we develop an experimental platform and quantitative framework in which pluripotent cells explanted from blastula stage Xenopus embryos can be used to study the transit of these cells to four different lineage states—epidermis, neural progenitor, endoderm and ventral mesoderm—by following changes in the transcriptome at six time points during this seven-hour process. These data provide quantitative insights into the dynamics of Waddington’s landscape and complement previous data sets looking at transcriptome dynamics in Xenopus [72, 73] by providing lineage specific dynamics, as well as transcriptome dynamics in response to BMP/Activin signaling. Our findings shed light on why a neural progenitor state is the default lineage state for pluripotent cells, uncover novel components of lineage specific GRNs, and provide insights into essential control of the timing of signaling inputs such as BMP for proper lineage decisions. These time-resolved data sets will serve as an important resource for future studies of developmental decision making in early vertebrate embryos.

Results

Naïve animal pole cells from Xenopus blastula can be programmed to any lineage state

To allow interrogation of transcriptome changes at high temporal resolution as pluripotent cells become lineage restricted, we established a highly regimented protocol for collecting synchronous populations of blastula explants across six points on the path towards lineage restriction. Late blastula stage explants (Nieuwkoop and Faber stage 9) were designated time zero and represent the pluripotent state atop Waddington’s landscape (Fig. 1A). In addition, explants were collected at 75, 150, 225, 315 and 435 min after stage 9 (Nieuwkoop and Faber stages 10, 10.5, 11, 12, 13), confirmed by stages of sibling embryos [74]. Stage 10, the onset of gastrulation, is marked by presence of the dorsal blastopore lip (Fig. 1B); stages 11 and 12 are mid and late gastrulae respectively. At stage 13, the neural plate stage, the developmental potential of embryonic cells has been restricted to a single lineage state, with the notable exception of neural crest cells [75]. Explants cultured without additional instructive cues transit to an epidermal state (Fig. 1C). To follow cells as they transited to an endodermal state, explants were treated with Activin [57]. Titration experiments determined that treatment with 160 ng/uL of Activin at stage 9 was the minimum concentration of Activin able to robustly induce endoderm without expression mesodermal markers, as assayed by qPCR. To direct cells to a neural progenitor state explants were exposed to Noggin [27] 100 ng/uL of Noggin was used for these experiments as it was the lowest concentration that induced Sox3 at stage 13 and effectively blocked BMP signaling as determined by qPCR and Western blot. Treatment with BMP 4/7 heterodimers was used to induce a ventral mesoderm state [69, 76]. 20 ng/uL of BMP 4/7 was found to induce pSmad 1/5/8 at near endogenous levels, as determined by western blot, shifting the timing but not amplitude of BMP signaling. Use of BMP4/7 heterodimers to induce ventral mesoderm allows the transcriptional responses to the two different branches of TGF-beta signaling to be compared. RNA was isolated from explants at each time interval and used to generate illumina libraries for transcriptome analysis.

Fig. 1

Xenopus blastula explants can be reprogrammed to adopt any lineage. A Schematic of Waddington’s Landscape portraying lineage specification process. B Embryos stage 9,10,10.5,11,12, and 13 used to confirm developmental stages of blastula explants. C Schematic showing signaling molecules used to generate all 4 lineages. D-I RNA Seq TPM expression over time of (D) maternally provided pluripotency marker Pou5f3.3, (E) maternally provided Foxi2, (F) epidermal marker EpK, (G) neural marker Otx1, (H) mesoderm marker Brachyury(T) (I) endoderm marker Sox17. H Graphs are sums of S + L allele. Width of lines represents SEM of three biological replicates

The transcript dynamics of maternally provided pluripotency factors Pou5f3.3 and Foxi2 across three biological replicates for all four lineage transitions (twelve independent experiments) demonstrates that this pipeline generated highly quantitative and highly reproducible data with minimal technical error (Fig. 1D, E). Pou5f3.3 and Foxi2 are representative of a class of 119 genes whose expression decreases monotonically by at least 25-fold between stage 9 and stage 13 (Fig. S1A, Supplementary Table S1) and includes maternal transcripts that characterize the pluripotent state [77, 78]. Normalized RNA read count data, quantified as transcripts per million (TPM), confirmed generation of each of the expected lineage states, with EpK and Dlx3 validating establishment of the epidermal state (Figs. 1F and S1B), Otx1 and Otx2 validating transit to a neural state (Figs. 1G and S1C). Brachyury(T) and Evx1 (Figs. 1H and S1D), and Sox17 and Endodermin (Figs. 1I, S1E) validated the mesoderm and endoderm data sets respectively.

Global analysis of transcription factors expressed at each stage revealed that the majority are expressed in all four lineages, although their expression levels or timing may vary between them. For example, Otx1 is unique to the endodermal lineage at stage 10 but to the neural lineage at stage 13 (Fig. S2A-E, Supplemental Table 2). This analysis also identified the transcription factors that at any given stage are expressed in only a single lineage and revealed that such genes are overrepresented in the endoderm lineage. We performed weighted gene correlation network analysis (WGCNA) on transcription factors for each lineage independently and with time as a continuous variable in order to identify modules of genes with similar dynamics over time in each lineage. This analysis showed that many known Activin induced transcription factors cluster together in the same module not only in the endoderm lineage, but also in the mesoderm lineage, with about half of these also clustering together in the neural lineage. Similarly, many transcription factors known to be induced by BMP cluster together in both mesoderm and epidermis, and surprisingly a small subset of these also cluster in the neural lineage, suggesting that there are subsets of transcription factors that will cluster together (i.e., exhibit similar dynamics to one-another) even in response to distinct inductive signals (Supplemental Table 3). Further study of these putative modular gene regulatory networks may provide novel insights into potential mechanisms of co-regulation of correlated factors.

PCA and time series analysis reveal novel lineage-specific dynamics

Together the transcriptomes of the four state transitions each across six time points yield 72 observations in a 45,661-dimensional gene expression space. Global insights into such high dimensional data require methods for dimensionality reduction. Principal Component Analysis (PCA) can provide key insights into the genes contributing most significantly to the variance between lineage states and developmental stages. We first used PCA to analyze each lineage individually, plotting the first two principal components against developmental time (Fig. 2A). For all four lineages the primary principal component (PC1) was found to be largely monotonic over time, suggesting that the majority of the gene expression variance is contributed by genes changing unidirectionally, such as pluripotency genes being turned off or lineage-specific genes being activated. Interestingly, when PC2 was plotted for each of the lineages it was found to exhibit temporal dynamics in the epidermal, ventral mesodermal and endodermal lineages that suggested genes exhibiting expression peaks at intermediate time points make a significant contribution to the variance, and potentially to these state transitions. By contrast, PC2 of the neural lineage shows no association with time, suggesting that the temporal dynamics in the neural lineage are primarily linear in time. This raises the possibility that the transition to a neural progenitor state follows a simpler trajectory than that of the other three lineages (Fig. 2A).

Fig. 2

PCA and Time Series Analysis Reveal Novel Lineage-Specific Dynamics. A PCA for each individual lineage with the coordinates of PC1 and PC2 for each lineage plotted against developmental time. B PCA performed on all four lineages simultaneously, with plot showing clustering of all lineages for PC1 vs PC2. C Scree plot of the variance explained by the top10 principal components for PCA done on all lineages. D Distances from stage 9 to 13 for all 74 PCs for each lineage, error bars are SEM of all 9 stage 9–13 distances for each lineage, (***P < 0.005, **P < .01). E Number of genes differentially expressed between successive stages in each lineage, p_adj ≤ 0.05

To gain further insights into these state transitions, PCA was carried out on all four lineages in concert. In this analysis, PC1 and PC2 together explain 75% of the variance across these data. (Fig. 2B, C). The distribution of these data across the PC1 axis corelates with developmental time, with the earliest (St.9) samples clustering at the top of the plot and the later samples progressively further down the y-axis (Fig. 2B). When examined in this context, the neural trajectory is striking as it extends a shorter distance along this axis, with the neural replicates for stages 11, 12 and 13 clustering closer to the stage 10.5 replicates for the other lineages. To quantify this, the distance between stages 9 and 13 in the full gene expression space was calculated and here too the neural lineage was found to move the shortest distance (Fig. 2D). Whereas PC1 correlates with developmental time, PC2 appears to distinguish the different lineage states, which after stage 10.5 show very distinct trajectories. Interestingly, endoderm and epidermis lie furthest from each other along PC2 with mesoderm lying between those states. Thus, PC2 captures the intermediate nature of the ventral mesoderm state which shares GRN features with endoderm but, like epidermis, is BMP-driven (Fig. 2B).

As some maternally provided transcripts persist through stage 10.5 and could bias transcriptome changes, we also performed PCA on zygotic transcripts only, removing maternally provided genes [79] with no zygotic transcription through state 13 from the analysis (Fig. S2F). The PCA was largely unchanged using this gene set, indicating that the observed dynamics are driven primarily by zygotic transcripts.

We next used differential expression analysis (DESeq2) to gain insights into the dynamics of gene expression changes across the four state transitions. Plotting the number of genes differentially expressed between successive developmental stages revealed that the number of genes whose expression changes significantly during these state transitions is relatively modest (Fig. 2E, Supplemental Table 4). For example, between two and three thousand genes are differentially expressed between stages 9 and 10 in each lineage, which represents four to six percent of the transcriptome. Between stages 10 and 10.5 the gene expression changes in the endodermal lineage are strikingly different from those of the other lineages. The number of differentially expressed genes increases almost 60% between these stages in the endoderm, whereas there is a significant decrease in differentially expressed genes in the epidermal and neural lineages and little change in the mesoderm. While there is a gradual decrease in differentially expressed genes in the endodermal lineage between stages 10.5 and 13, this state transition continues to exhibit the most dynamic changes in gene expression compared to the other lineages. After reaching a minimum between stages 10.5 and 11, both the epidermal and ventral mesodermal lineages exhibit increasing numbers of differentially expressed genes. By stages 12–13 the endoderm and ventral mesoderm exhibit comparable gene expression dynamics. Interestingly, almost a quarter of the genes changing in the endodermal and ventral mesodermal trajectories between stages 12 and 13 are shared between these lineages, likely reflecting the overlapping landscape of the combined mesendoderm GRN (Fig S3I-J). In contrast to the other three lineages, the neural trajectory exhibits very few differentially expressed genes after stage 11, providing additional evidence that this lineage reaches an early equilibrium (Figs. 2E, S3G-J, Supplemental Table 5).

Gene expression dynamics provide novel insights into the neural default state

The above analyses suggested that the neural progenitor state follows a simpler trajectory than that of the other three lineages. Comparison of gene expression dynamics during transit to an epidermal versus a neural state reveals that through stage 10.5 these two lineages share a remarkably similar trajectory (Fig. 3A). For each, the number of genes differentially expressed between successive stages decreases, the number of genes changing is highly similar, and there is significant overlap in the genes exhibiting differential expression. Between stages 9 and 10, for example, 60% of the genes differentially expressed in the neural linage are also differentially expressed in the epidermal lineage and that is true of 53% of genes differentially expressed between stages 10 and 10.5 (Fig. 3B, Supplemental Table 6). The majority of shared genes exhibit decreasing expression during these stages, in part reflecting the downregulation of pluripotency genes. Nevertheless, of the 1161 genes whose expression increases in the neural trajectory between stages 9 and 10, more than half also show increased expression in the epidermal trajectory.

Fig. 3

Gene Expression Dynamics Provide Novel Insights into the Neural Default State. A Number of differentially expressed genes between successive developmental stages of the epidermal and neural lineage, p_adj ≤ 0.05 B Venn Diagrams for the total number of DE genes between stages for the epidermal and neural lineages, as well as the number of genes increasing and decreasing over time. C TPM of Foxi1, Foxi2 and Grhl1, revealing epidermal onset dynamics. Graph shows sum of S + L allele. Width of the line represents SEM of three biological replicates. D Number of differentially expressed genes between lineages at each developmental stage (p_adj ≤ 0.05). E Graphs of the TPM of three pluripotency markers maintained in the neural lineage (E) Sox3, Sox11, and Zic1. Graphs are sums of S + L allele. Width of the line represents SEM of three biological replicates

Between stages 10.5 and 11, which corresponds to early gastrulation, there is a striking divergence of the epidermal and neural lineages. During these stages, gene expression dynamics largely cease in the neural lineage; fewer than 100 genes are differentially expressed between stages 10.5 and 11, and only 52 genes are differentially expressed between stages 12 and 13. By contrast, in the epidermal lineage temporal changes in gene expression begins to sharply increase and the overlap of these genes with those changing in the neural lineage is minimal (Fig. 3A, B). Of all genes differentially expressed between stages 12 and 13 in these two stage transitions, less than 1% are differentially expressed in both the neural and epidermal lineages, mainly because the neural lineage has ceased to change. Interestingly, the increasing gene expression dynamics that characterizes the epidermal lineage at stage 11 coincides with a loss of enrichment for neural GO terms (Fig. S4A, Supplemental Table 7). This suggests that pluripotent blastula cells possess neural-like features that begin to be lost around the onset of gastrulation as cells transit to an epidermal state but are retained and reinforced in neural progenitor cells.

The observed gene expression dynamics and GO term enrichment together point to the onset of gastrulation as a critical point on the landscape topology of early developmental decision making. To further explore this, we examined the expression of Foxi1 and Grhl1 which are key upstream components of the GRN mediating the formation of epidermis [80, 81]. Expression of Foxi1 has been shown to be activated by the pluripotency factor Foxi2 [82]. Examining the expression dynamics of these three genes across the epidermal trajectory reveals a sharp increase in the expression of Grhl1 that correlates with rapidly extinguishing expression of Foxi2 and the onset of gastrulation (stage 10) (Fig. 3C). As early gastrulation is also when gene expression dynamics virtually cease in the neural trajectory, this suggests systems dynamics that favor a neural progenitor state over that of other lineages and a network structure that requires cells to be actively propelled toward an alternative, in this case epidermal, state.

To further examine the genes that distinguish the neural and epidermal states we analyzed the genes differentially expressed between these states at each time point on their trajectories. At stage 10 the two lineages remain strikingly similar, with only 261 genes differentially expressed (Fig. 3D, Supplemental Table 8). The number of differentially expressed genes increases by more than 500% between stages 10.5 and 12, driven almost entirely by gene expression dynamics in the epidermal lineage. Interestingly, 13 of the top 20 most differentially expressed genes at stage 10.5 are known BMP responsive genes and all but three, Jun.L/S and Actc1.S, are more highly expressed in the epidermal lineage (Fig. S4B). Consistent with this, beginning at stage 10.5 genes associated with the TGF-beta pathway in the KEGG database show enrichment in the epidermal lineage (Fig. S4C). As stage 10.5 represents the time when the trajectories of the neural and epidermal lineages diverge after neural reaches early equilibrium, this enrichment is consistent with a model where BMP signals actively propel cells away from a neural well, the state favored by the systems dynamics absent perturbation of the landscape, and onto the path toward an epidermal state. Together, the early equilibrium reached by the neural lineage, combined with the neural features of the pluripotent state provide new context for the neural default hypothesis. Complementing prior experimental studies, our transcriptome data suggests that neural is the default state following exit from pluripotency because of the shorter and more linear path from the pluripotent state to the neural progenitor state. This is further evidenced by the expression dynamics of genes that play important roles in both pluripotent cells and neural progenitors such as Sox3, Sox11 and Zic1 [43, 44, 83, 84]. All three of these genes retain or increase their expression in the neural lineage but are rapidly down-regulated in the epidermal lineage after stage 10.5 (Fig. 3E).

Robust BMP signaling is initiated in explants around the onset of gastrulation

Consistent with a model where BMP signals actively propel cells away from the neural state, phosphorylation of BMP R-Smads is first detected in animal pole explants at stage 10.5 (Fig. 4A). Translocation of pSmad1/5/8 to the nucleus drives expression of BMP responsive genes, and its timing correlates with the divergence of the epidermal lineage from neural lineage. To gain further insights into the timing of BMP responsiveness we examined whether genes differentially expressed at successive developmental stages displayed over-representation of genes associated with BMP signaling [85] using the DESeq2 Wald test. We computed the significance at which these BMP associated genes comprised a larger fraction than would be expected by random chance via the hypergeometric p-value [86]. The greatest divergence in overrepresentation between the epidermal and neural lineages was seen between stages 10.5 and 11, which was also the maxima for overrepresentation in the epidermal data (Fig. 4B). Consistent with this finding, Ventx2.1 and Id3, which are both BMP targets genes, exhibit expression maxima in the epidermal lineage and minima in the neural lineages at these stages (Fig. 4C, D) [87, 88]. The expression of Id3 across these two state transitions is particularly noteworthy for its opposite intermediate non-monotonic dynamics at successive developmental stages despite comparable expression at the start and end of these lineage trajectories.

Fig. 4

Robust BMP Signaling is initiated in Explants Around the Onset of Gastrulation. A Western blot analysis of lysates of developing epidermal (WT) and neural (20uM K02288) explants for pSmad1/5/8 and Smad1 with Actin loading control. Gel was cropped to minimize empty space, uncropped gel is Supplemental Fig. 8A. B Significance of BMP overrepresentation (hypergeometric p-value) in temporally differentially expressed genes. C-D Graphs of BMP responsive genes in epidermal and neural lineages (C) Ventx2.1 (D) Id3. Graphs are sums of S + L allele. Width of the line represents SEM of three biological replicates

Early response to Activin and BMP4/7 reveals unexpected overlap

While BMP signaling plays a central role in instructing pluripotent cells to form epidermis, it is the other branch of the TGF-beta that directs mesendodermal fates. Members of the Activin/ Nodal/Vg1/GDF1/TGF-beta subfamily act via pSmad2/3 to activate mesodermal and endodermal target genes including Foxh1, Eomes, Mixer, Tcf3 (also known as E2a) and Tp53 [89,90,91,92,93,94]. While recent work has suggested that Vg1-Nodal heterodimers mediate this process endogenously [95], Activin has long been used for efficient mesendoderm induction in ex-vivo assays [52, 53, 68]. Like Activin, BMP4/7 heterodimers have been shown to be potent mesoderm inducers at physiological concentrations. However, BMP4/7 has been reported to induce only ventral mesoderm [69, 96], unlike Activin and Nodal which are able to induce both ventral and dorsal mesoderm at low concentrations, as well as endoderm at high concentrations. We focused on BMP4/7-mediated mesoderm induction for this analysis because the response to Activin/Nodal is a spectrum with no clear threshold cleanly distinguishing an endodermal versus mesodermal response. An additional advantage of focusing on BMP4/7-mediated mesoderm induction is that it provides an opportunity to directly compare the signaling dynamics and transcriptional responses to the two branches of TGF-beta signaling in the same quantitative framework.

Treatment with Activin at stage 9 leads to robust signaling at stage 10 as evidenced by Western detection of phosphorylated Smad2 (p-Smad2) (Fig. S5A) and robust induction of Sox17 beginning at stage 10 (Fig. 1I). Interestingly, transit to an endodermal state is distinguished from the other lineage transitions by its unique transcriptome dynamics. It is the only lineage in which there is a large increase in differentially expressed genes between stage 9 and stage 10.5, early gastrulation, after which the number of differentially expressed genes decreases (Figs. 2E, 5A). By contrast, treatment of stage 9 explants with BMP4/7 does not elicit an immediate increase in differentially expressed genes, distinguishing the responses to the two different arms of TGF-beta signaling between successive developmental stages (Fig. 5A). Instead, the number of genes that are differentially expressed between stages 9–10 and 10–10.5 following BMP4/7 treatment remains fairly constant, before decreasing between stages 10.5 and 11. Intriguingly, however, the genes that are differentially expressed between these early stages in response to Activin or BMP4/7 show significant overlap. For example, approximately 48% of genes differentially expressed between stages 9 and 10 in response to BMP4/7 are also differentially expressed between those stages in response to Activin, as are 71% of genes differentially expressed in response to BMP4/7 between stages 10 and 10.5 (Fig. 5B, Supplemental Table 9). These findings were unexpected as Smad1/5/8 and Smad2/3 generally regulate distinct target genes [97, 98].

Fig. 5

Early Response to Activin and BMP4/7 Reveals Unexpected Overlap. A Number of differentially expressed genes between successive developmental stages of the endoderm and ventral mesoderm lineage, p_adj ≤ 0.05. B Venn Diagrams for the total number of DE genes between stages for the epidermal and neural lineages, as well as the number of genes increasing and decreasing over time. C-F Heatmaps of the top 50 genes differentially increased in response to (C) BMP4/7 at stage 10 and/or 10.5 (D) BMP4/7 at stage 11 and/or 12. (E) Activin at stage 10 and/or 10.5 (F) Activin at stage 11 and/or 12. All genes in heatmaps are ranked by Log2FC of differential expression of either mesoderm (BMP4/7) or endoderm (Activin) compared to the wildtype epidermal lineage. Only genes increasing in response to BMP4/7 or Activin with a minimum expression of 10TPM at the relevant stages and higher expression at the relevant stage than at stage 9 were included in the heatmap. Colors represent z-scores of TPM

To further explore the unexpected overlap in transcriptional responses to the two different classes of TGF-beta ligands we used DESeq2 to examine the genes that exhibit the largest log2 fold change compared to untreated explants in response to BMP4/7 treatment. Figure 5C shows the top 50 genes exhibiting the largest expression increase at stages 10 or 10.5 in response to BMP4/7 treatment relative to untreated explants. This gene set, which captures the immediate response to this ligand, includes previously characterized BMP target genes such as Msx1 and Ventx1 [99, 100], as well as the pan-mesodermal gene Brachyury(T) [101]. Unexpectedly, it also includes a number of dorsal mesoderm/endoderm genes that have been characterized as targets of Activin/Nodal signaling including Bix1, Mix1, Mixer and Sox17 [65, 89, 102, 103]. Notably, activation of these genes occurs absent activation of pSmad2/3 (Fig. S5A). Using z-score scaling to visualize the expression dynamics of early responding genes across the time series revealed that genes generally associated with Activin/Nodal expression, including Bix1 and Mix1, displayed non-monotonic dynamics with expression peaks at intermediate stages, whereas the expression of ventral and pan-mesodermal genes increased monotonically (Fig. 5C). After reaching a minima between stages 10.5 and 11, the number of genes displaying dynamic expression changes in response to BMP4/7 greatly increased (Fig. 5A). The genes exhibiting the largest log2 fold change at stages 11 or 12 were therefore similarly examined. This gene set is more enriched for ventral mesoderm associated genes than the initially responding genes, suggesting that the Activin-like response to BMP4/7 is transitory and that ventral mesoderm character is stabilized secondarily (Fig. 5D). This is consistent with a role for BMP signaling in actively ventralizing mesoderm and other tissues [104]. We compared our endoderm and ventral mesoderm trajectories to a recently published Activin-induced mesoderm time series [6] and found that the Activin treated samples cluster more closely relative to BMP4/7 treated samples (Fig. S5B). To ask if this is due to a greater enrichment of dorsal genes, we used a previously curated set of dorsally and ventrally induced genes [105]. 31.0% of dorsally enriched genes were upregulated in the Activin-induced mesoderm samples as compared to 11.7% of ventrally enriched genes. By contrast, BMP4/7 treatment led to upregulation of 18.7% of dorsally enriched genes and 12.6% of ventrally enriched genes. Hierarchical clustering of stage 11 and 13 samples z-scored for dorsal and ventral genes identified by Ding et al. further highlights the more ventral character of BMP4/7-induced mesoderm (Fig. S5C).

While treatment with both Activin and BMP4/7 was initiated at stage 9, Activin-mediated phosphorylation of Smad2 was transient (Fig. S5A) whereas BMP4/7-mediated phosphorylation of Smad1/5/8 persisted through stage 13 (Fig. 7A), likely contributing to the distinct gene expression dynamics in these two lineages. For example, the genes exhibiting the largest log2 fold change compared to untreated explants at stages 10, 10.5 or 11, 12 in response to Activin treatment are enriched for those whose maximal expression occurs at intermediate stages of the lineage trajectory (Fig. 5E, F). Importantly, while the genes activated as an early response to Activin or BMP4/7 show significant overlap (Fig. 5B, Supplemental Table 9), the genes activated by these two classes of TGF-beta ligands nevertheless show significant differential expression with respect to one another (Fig. S6A, Supplemental Table 8). Over 1000 genes are differentially expressed between these trajectories as early as stage 10.5, and the number of genes differentially expressed between these lineages continues to increase over developmental time. Interestingly, analysis of KEGG pathway enrichment in these differentially expressed genes reveals that genes that are significantly higher in the endoderm lineage show enrichment for the TGF-beta pathway at stages 10 and 10.5, whereas genes that are significantly higher in the ventral mesoderm lineage are enriched for the TGF-beta pathway at stages 11 and 12 (Fig. S6B).

Time series data provides novel insights into mesendoderm GRN

The mesendoderm gene regulatory network (GRN) has been extensively studied and has yielded a significant “parts list” of genes that make significant contributions to the formation of these lineages [106, 107], however the ordering of the GRN components has lacked the temporal resolution that our time series data can provide. Accordingly, we examined the expression of forty-one validated mesendoderm GRN components across both the endoderm and ventral mesoderm lineage trajectories (Fig. 6A, B). Interestingly, many of these genes display non-monotonic expression, with their expression peaking at early time points before decreasing. These dynamics are particularly striking in the endoderm, and demonstrate that many of these GRN components respond to Activin rapidly and robustly, but transiently. Indeed, expression of twenty-one of these genes is induced in the endoderm by stage 10, which is 75 min after ligand exposure (Fig. 6A, B). Interestingly, sixteen of these genes are also activated by BMP4/7 by stage10, albeit less robustly, indicating that they are immediate responders to both classes of TGF-beta ligands. The Activin-induced expression of several GRN components, including Bix and Nodal family genes, Vegt, Eomes, Mix1 and Snai1, peaks by stage 10.5 and then declines. By contrast, a second group of GRN factors, including Sox17 and Gata2/6-related genes, Osr2 and Pitx2 exhibit sustained expression through stage 13. A third group of GRN components, including Hnf1b, Foxa2 and Gata5, are not robustly turned on until mid- to late-gastrula stages (Fig. 6A, B). The expression of Ventx genes, known for their strong ventralizing activity, in response to high levels of Activin signaling was unexpected, and correlates with the downregulation of Nodal and Bix family factors. Surprisingly, Ventx2.1 and Ventx2.2, which are among the first genes to respond robustly to BMP4/7, were induced as or more strongly by Activin, albeit with different temporal dynamics. Subsequently GRN factors expressed highly in mesoderm, including Bix1.3, Mix1 and Evx1, distinguish the BMP4/7 response from the Activin response. As with the Activin response, these BMP4/7 responding genes show distinct patterns of temporal dynamics (Fig. 6C).

Fig. 6

Time Series Data Provides Novel Insights into Mesendoderm GRN. A Heatmaps of genes in the published mesendoderm GRN across time in the endoderm and ventral mesoderm lineage with colors representing z-score of TPM expression in both lineages [106]. (B) Schematic of the timing of genes from published Mesendoderm GRN in both the endoderm and ventral mesoderm lineage, as defined by expression of at least 30 TPM in the L and S allele combined for the average of three biological replicates. C TPM of Bix1.3, Evx1, and Mix1 revealing mesoderm onset dynamics. D-G TPM of genes proposed as novel mesendoderm GRN members based on DESeq2 and Limma analysis (D) Tmcc1 (E) Nptx2 (F) Ca14 (G) Pygm. Graphs are sums of S + L allele. Width of the line represents SEM of three biological replicates

Given the distinct dynamics displayed by known mesendoderm GRN factors, we tested for differential linear and quadratic dynamics in the BMP4/7 and Activin responses using the R package limma [108]. This analysis identified genes displaying expression dynamics that position them as candidates for novel components the mesendoderm GRN. For example, one noted pattern was genes expressed rapidly and robustly in response to Activin but not BMP4/7, and downregulated after an early peak in expression. This pattern is exemplified by Tmcc1 in which is expressed non-monotonically in the endoderm trajectory (Fig. 6D). A second pattern described genes that respond to both Activin and BMP4/7 with transient non-monotonic expression, such as Nptx2 (Fig. 6E). A third pattern that emerged from this analysis was genes displaying a monotonic increase in expression only to Activin, such as Ca14, suggesting a role in the endoderm lineage specifically (Fig. 6F). Similarly, genes with monotonic and sustained expression in response to BMP4/7, such as Pygm may be novel mesoderm regulatory factors (Fig. 6G). While these genes are largely unstudied in mesendoderm formation, published transcriptome data sets provide further support for their involvement in mesendoderm formation. For example, Ca14, Pygm and Tmcc1 were among genes upregulated in stage 12 animal caps in response to Wnt and Nodal2 [109] and Ca14, Tmcc1, and Nptx2 were upregulated in stage 11 embryos in response to somatic cell nuclear transfer from an endoderm cell [110]. Similarly, Pygm has been identified as a target of Myod [111]. Thus, using limma analysis as an unbiased approach for detecting genes sets that share expression pattern dynamics allows identification of potential new members of developmental GRNs using our data sets.

To further examine the transition from a pluripotent state to endoderm versus mesoderm we performed WGCNA on these lineages using data from stages 10, 10.5 and 11 when these state transitions display very distinct dynamics. Two notable gene clusters emerged from WGCNA (blue and brown modules, Fig. S6C). These modules revealed correlations that increased over developmental time in either the endoderm (blue) or mesoderm (brown). Of the forty-one established mesendoderm GRN members, thirty were included in one of these two modules. Moreover, three of the novel GRN candidates identified by limma analysis were also found in one of these two clusters. The fourth of these was contained in a different module that also contained mesendoderm genes but had lower correlation values (turquoise). Identification by WGCNA analysis of a cohort of genes significantly correlated with known mesendoderm GRN members provides candidate genes that may also play important roles in these two lineages (Supplemental Table 10).

Early BMP signaling drives ventral mesoderm rather than early epidermal divergence

As discussed above, endogenous BMP4/7 signaling within animal pole cells will direct these cells to give rise to epidermis in the absence of BMP inhibitors, which in vivo are secreted by the organizer. BMP signaling is detectable in these cells by stage 10.5, as evidenced by detection of pSmad1/5/8 (Fig. 4A). When explants are treated with exogenous BMP4/7 at stage 9, pSmad1/5/8 is detected at stage 10 at levels comparable to those seen at stage 10.5 in untreated explants (Fig. 7A). This allows a quantitative comparison of the transcriptional responses to the same signal when presented with shifted developmental timing – a tilting of the landscape topology. One predicted outcome of such a heterochronic shift might have been an accelerated transit to the epidermal state rather than formation of ventral mesoderm. To examine how shifting the timing of BMP activity alters the transcriptional response we first compared the transcriptome dynamics. In this context it is interesting to note that premature BMP signaling actually dampened early gene expression changes – there is an ~ 32% reduction in the number of genes whose expression significantly changes between stages 9 and 10 (Fig. 7B). This is driven almost entirely by a reduction in the number of genes whose expression decreases during those initial stages (Fig. 7C, Supplemental Table 11). Between stages 10 and 11, the transcriptome dynamics are comparable in the two conditions (Fig. 7B, C). However, the number of genes displaying temporal differential expression that are shared by these two trajectories but not by the neural and endodermal state transitions (making them a general response) is remarkably low, ranging from 97 between stages 10 and 10.5, to 55 from stages 10.5 to 11 (Fig. S3C-F, Supplemental Table 5).

Fig. 7

Early BMP signaling drives ventral mesoderm rather than early epidermal divergence. A Western Blot Analysis of lysates for developing epidermal (WT) and ventral mesoderm (BMP4/7 20 ng/uL) explants for pSmad1/5/8 and Smad1 with Actin loading control. Gel was cropped to minimize empty space, uncropped gel is Supplemental Fig. 8B. B Number of differentially expressed genes between successive developmental stages of the epidermal and ventral mesoderm lineage, p_adj ≤ 0.05. C Venn Diagrams for the total number of DE genes between stages for the epidermal and ventral mesoderm lineages, as well as the number of genes increasing and decreasing over time. D-I TPM of genes representative of different expression dynamics in response to early BMP (D) Grhl1 (E) Ventx2.1 (F) Post (G) Vegt (H) Wnt8a (I) Mix1. Graphs are sums of S + L allele. Width of the line represents SEM of three biological replicates

We next used DESeq2 to compare the genes differentially expressed in response to early BMP signaling. At stage 10 there are only 332 genes differentially expressed between these two conditions. The number of differentially expressed genes remains relatively low until stage 11, when it begins to increase over time (Fig. S7A, Supplemental Table 8). However, beginning at stage 10 the genes whose expression is higher in response to early BMP are enriched for mesoderm GO terms (Fig. S7B, Supplemental Table 12). Because expression of Grhl1 is a key early driver of transit to the epidermal state (Tao et al. 2005), we examined its response to early BMP signaling. Strikingly, its expression fails to be robustly activated under this condition and the differential response is seen as early as stage 10.5 (Fig. 7D). Given this surprising finding we examined the expression of Ventx2.1, which is a direct target of BMP signaling [87]. Here there was a shift in transcription dynamics that reflected the earlier onset of BMP signaling; Ventx2.1 transcripts reach levels at stage 10 that would not be achieved until state 10.5 in response to endogenous BMP signals (Fig. 7E). Importantly, this demonstrates that BMP signaling is able to immediately elicit changes in gene expression when activated prematurely, and that this drives a heterochronic response in expression of some BMP target genes. Another BMP responsive gene, Post, which plays a role in conferring posterior/ventral attributes to both ectoderm and mesoderm [112], did not show a premature onset of expression, but instead displayed a significant increase in its amplitude of expression (Fig. 7F). Most striking, however, was the activation of genes categorized as Activin/Nodal target genes including Vegt, Wnt8a, and Mix1 (Fig. 7G-I) [103, 113, 114] Importantly, activation of these genes in response to early BMP is not due to the inappropriate activation of pSmad2/3 (Fig. S5A). We also confirmed that is the early timing and not the exogenous nature of BMP4/7 exposure that drives ventral mesoderm formation. While treatment of stage 10.5 explants with BMP4/7 leads to increased pSmad1/5/8 (Fig. S7C), it does not lead to expression of mesodermal genes (Fig. S7D).

BMP signaling is restrained until stage 10.5 by Dand5 activity

The striking finding that shifting the timing of BMP signaling leads to activation of Activin/Nodal-responsive genes, and the formation of mesoderm instead of epidermis, indicates that it is essential to tightly control when cells receive this signal endogenously. BMP2,4 and 7 ligands, as well as the receptor Smad1 are expressed at stage 9 and 10 (Fig. 8A) and cells are clearly competent to respond to early BMP signals as evidenced by the shifted activation of Ventx2.1 (Fig. 7E). Despite this, however, pSmad1/5/8 is not robustly detected in control explants until stage 10.5. We therefore investigated how BMP signaling is restrained in animal pole cells until the onset of gastrulation. We asked if there were BMP signaling antagonists expressed in blastula animal pole cells that were downregulated with dynamics consistent with the observed timing of pSmad1/5/8 accumulation. We identified two maternally provided BMP antagonists, Dand5 and Gtpbp2 [115,116,117,118,119], that are expressed at stage 9 but downregulated to significantly lower levels by stage 10.5, when BMP activity is observed (Fig. 8B). Dand5, in particular, is robustly expressed at stages 9 and 10 in animal pole cells. To determine if Dand5 plays a role in preventing premature BMP signaling we used a translation blocking morpholino to deplete it from early blastulae. We found that Dand5-depletion led to premature phosphorylation of Smad1/5/8, indicative of early BMP signaling (Fig. 8C), as well as expression of mesodermal marker Bra (T) at stage 13 (Fig. 8D). This suggests that a key role for Dand5 at these stages is preventing premature BMP signaling that would generate a mesodermal rather than ectodermal state transition. As Dand5 depletion does not increase pSmad1/5/8 or Brachyury levels to the same degree as BMP4/7 treatment, other BMP antagonists including Gtpbp2, likely cooperate in temporally constraining BMP signaling.

Fig. 8

BMP Signaling is Restrained Until stage 10.5 by Dand5 Activity. A-B TPM of genes involved in BMP signaling (A) BMP heterodimer ligands BMP7.2, BMP4, BMP2 and primary BMP target Smad, Smad1 (B) maternally provided BMP antagonists Dand5 and Gtpbp2. Graphs are sums of S + L allele. Width of the line represents SEM of three biological replicates. C Western Blot Analysis of lysates for stage 10 epidermal (WT), Dand5 MO injected (40 pmol/embryo), ventral mesoderm (BMP4/7 20 ng/uL), and endoderm (160 ng/uL Activin) explants for pSmad1/5/8 and pSmad2 with Actin loading control. Gels were cropped to minimize empty space, uncropped gel is Supplemental Fig. 8C. D qRT-PCR of animal pole explants examining the fold change from stage 9 to 13 of expression of mesodermal marker Bra(T) for epidermis(WT), ventral mesoderm(BMP4/7 20 ng/uL) and Dand5MO(300pmo/embryo) injected(***P < 0.005)

Discussion

How progenitor cells decide their fate is a fundamental question central to all of developmental biology. While in many cases the inductive cues that drive these decisions have been identified, less is understood about how the timing of these signals is controlled and the dynamics of the transcriptional circuitry they activate. A related and fascinating question is how a large and complex set of transcriptional responses is canalized into a discrete set of lineage trajectories.

Decades of research in Xenopus and other systems has shed important light not only on the signaling pathways driving lineage formation in early embryos, but also on many of the key transcriptional targets of these signals. Combined with gain and loss of function studies for individual factors, this work has allowed the construction of putative gene regulatory networks (GRNs) depicting how different lineage states are adopted. A powerful strength of the Xenopus system is the ability to easily isolate pluripotent cells from blastula embryos and culture them in simple saline with or without added inductive cues, and that these cells become lineage restricted over a period of approximately seven hours. This allows the transcriptional responses to inductive cues to be quantified with high temporal resolution, enabling the dynamics of individual genes as well as entire lineages to be followed. We used the pipeline we established here to study the transit of initially pluripotent cells to four distinct lineage states using bulk transcriptomics. However, this pipeline can be built upon to layer on additional analyses including ATAC to follow changes in chromatin accessibility and ChIP-Seq to follow changes in epigenetic marks dynamically during lineage restriction.

Neural as the default state

Our analysis of four different state transitions lends unexpected support for the neural default model and provides novel insights into why neural is the default state for pluripotent cells. Two types of analysis immediately distinguished the neural lineage from the other three state transitions. First, temporal DESeq reveals that this lineage reaches a steady state in gene expression dynamics by stage 10.5; after this time point the number of genes exhibiting differential expression is quite small (Fig. 2E). This is in marked contrast to the epidermal and mesodermal state transitions which become increasingly dynamic after stage 10.5. Similarly, principal component analysis reveals that the neural state lies closer to the pluripotent state than the other three lineages (Fig. 2B). This is true for both the PC1 and PC2 axes, which together explain 75% of the variation in gene expression. Since developmental times correlate along the PC1 manifold while state identities correlate along the PC2 manifold, this confirms that the time taken to transit from pluripotent to neural state is shorter than for the other lineage trajectories. Measuring the distance between stages 9 and 13 across all 74 principal components also indicates that the neural lineage is most positively correlated with pluripotency (Fig. 2D). Thus, neural progenitor cells occupy a position closest in state space to pluripotent cells relative to the other lineages.

The neural default model has not been without controversy. Studies in avian embryos have suggested that BMP inhibition may not be sufficient for transit to a neural progenitor state [38, 120]. While work in mESCs and hESCs provided additional support for the neural default model [121, 122], these cells do not necessarily recapitulate the in vivo state of inner cell mass cells. Our findings thus provide important validation of this model. Interestingly, when the neural lineage trajectory is compared to that of the epidermal trajectory they are highly correlated through stage 10.5, early gastrulation. This is true with respect to their dynamics, as evidenced by temporal DESeq (Fig. 3A), and also supported by the very small number of genes that exhibit differential expression between these states until stage 10.5 (Fig. 3D). The onset of gastrulation can thus be considered a point in time when a group of equipotent cells (neural/epidermal) diverge and either continue changing state to become epidermal or do not continue changing state and become neural. Indeed, individual cells at blastula stages and early gastrulation have been found to display multilineage gene expression [123].

Significantly, stage 10.5 is when robust BMP signaling is first detected in animal pole cells, as evidenced by pSmad1/5/8 detection (Fig. 4A), and this correlates with when the trajectory of the epidermal lineage diverges from that of neural (Fig. 3A). It is also the time when we observe a sharp increase in expression of Grhl1, a key upstream component of the epidermal GRN [81], significant down-regulation of the pluripotency factor Foxi2 [82](Fig. 3C) and a loss of enrichment for neural GO terms in the epidermal lineage (Fig. S4A). Neural features are retained and enhanced in the Noggin-treated explants, as exemplified by the expression dynamics of transcription factors Sox3, Sox11 and Zic1 (Fig. 3E). By contrast the BMP target genes Ventx2.1 and Id3 exhibit expression maxima in the epidermal lineage and minima in the neural lineages around stage 10.5 (Fig. 4 C, D). Given the distinct non-monotonic dynamics of Id3 in the epidermal and neural lineages it is tempting to speculate that this inhibitory bHLH factor may be playing a role in suppressing the function of neuralizing factors in the prospective epidermis thus helping to canalize this state transition.

Surprising overlap in transcription response to BMP4/7 and Activin

Our data sets allow direct comparison of the transcriptome changes driven by the two different branches of the TGF-beta superfamily, BMP and Activin. A striking feature of the Activin-driven endoderm trajectory is its distinct dynamics, characterized by a large increase in differentially expressed genes between stage 9 and stage 10.5, early gastrulation (Fig. 5A). This is not seen following treatment with BMP4/7. As expected, the characterized members of the mesendoderm GRN are induced in response to Activin and many of these respond rapidly and robustly, but transiently. Indeed, twenty-one of forty-one GRN factors examined are induced by stage 10, only 75 min after ligand exposure (Fig. 6A,B). PSmad2/3 is also robustly detected at this time, but is no longer detected by stage 13 indicating that the signaling response to Activin is transient (Fig. S5A). Distinct signaling dynamics for this branch of TGF-beta signaling may contribute to the markedly different transcriptome dynamics observed in the endodermal lineage as compared to the BMP-driven epidermal and mesodermal lineages (Fig. 2E). Among the genes that respond immediately to Activin are several members of the Bix/Mix family of transcription factors [124]. Somewhat unexpectedly, genes generally associated with ventral fates, including Ventx2.1, Ventx2.2 and Wnt8, were also induced by high levels of Activin. Induction of these factors occurred at later points in the trajectory, and their expression correlates with the downregulation of Bix/Mix family genes and other transiently responding factors, including Eomes, Nodal, and Vegt. Going forward it will be of interest to determine if these ventralizing factors play a direct role in downregulating the expression of the endodermal factors that are expressed only transiently. By contrast, Sox17 responds immediately to Activin and its expression increases linearly through stage 13. Endodermin also has a linear response to induction, although the increase in its expression does not commence until state 10.5, possibly reflecting a role for Sox17 in its activation. The distinct dynamics of early responding endoderm genes allowed the identification of putative new members of the endoderm GRN using limma analysis (Fig. 6D-G).

Among the most surprising findings emerging from these studies was the activation by BMP4/7 of genes that are generally characterized as Activin/Nodal targets. Indeed, among the earliest responses to BMP4/7 were Mix/Bix family genes, and similar to their response to Activin, their induction was transient (Fig. 5C). Other unexpected responding genes included Sox17, Eomes and Gsc. As Smad2/3 phosphorylation was not observed in response to BMP4/7, this suggests that the BMP R-Smads are capable of activating expression of these Activin/Nodal targets given a permissive cellular context. In this respect it is worth noting that Mix1.1 was previously identified in a screen for BMP4 responsive genes [125], supporting our current findings. It is intriguing that the genes exhibiting immediate responsiveness to BMP4/7 are dorsal mesendoderm factors, whereas pan and ventral mesoderm genes, including Bra (T), Wnt8, Post, Msx1 and Evx1, are turned on later in the trajectory. This suggests that the initial response to BMP signaling is “dorsal” as it is for Activin, and that more “ventral” attributes are a secondary response.

The timing of BMP signaling is critical for proper lineage segregation

The level of Bmp4/7 signaling utilized in these experiments was selected to match the level of pSmad1/5/8 levels present in untreated explants at stage 10.5 (Fig. 7A). This allows comparison of the response to the same signal and amplitude but with shifted developmental timing. Receiving the same level of BMP signaling at a slightly earlier time point might have been expected to accelerate transit to the epidermal state. Indeed, the shifted onset of Ventx2.1 expression is consistent with an accelerated response (Fig. 7E). However, explants also respond to BMP4/7 exposure at stage 9 by inducing expression of mesendodermal factors, including Vegt and Mix1(Fig. 7G, I), and by suppressing the endogenous BMP-mediated increase in expression of the epidermal regulatory factor Grhl1 (Fig. 7D). Thus, exposure to BMP4/7 at stage 9 elicits a fundamentally different transcriptional response than does exposure at stage 10.5. This was confirmed by treating explants with exogenous BMP4/7 at stage 10.5, which fails to elicit a mesendoderm response (Fig. S7D).

Together these findings indicate that it is critical to control the timing at which initially pluripotent cells are able to respond to endogenous BMP signaling. While we first detect pSmad1/5/8 at stage 10.5, it is likely that low levels of signaling are initiated by stage 10, as that is when increased expression of Ventx2.1 and Id3 is observed (Fig. 4C, D). Thus, cells undergo a fundamental change in competence between stages 9 and 10. Interestingly, expression of BMP inhibitors such as Noggin, Chordin, Follistatin, and Cerberus commences in the presumptive organizer region at late blastula stages [126], indicating that blocking BMP signaling in the marginal zone is also critical at these stages. This raised the question of how BMP signaling is restrained in blastula animal pole cells such that a mesendoderm response is prevented and ectodermal competence is established.

Using our data sets we identified two maternally provided BMP antagonists, Dand5 and Gtpbp2 that are expressed at stage 9 but are significantly downregulated by stage 10.5 (Fig. 8B). As Dand5 displayed significantly higher levels of expression, we examined the consequences of depleting it from initially pluripotent explants. Morpholino-mediated depletion of Dand5 depletion resulted in premature phosphorylation of Smad1/5/8 and expression of Bra (T) at stage 13 (Fig. 8C, D), indicating that a key role for Dand5 at these stages is preventing premature BMP signaling that would generate a mesodermal rather than ectodermal state transition. Interestingly, a role for Dand5 in controlling the spatial response to Activin/Nodal signaling had previously been suggested, restricting this response to the mesodermal mantle [115, 117, 118]. Our findings suggest a second centrally important role for this TGF-beta antagonist in the temporal control of BMP signaling as animal pole cells exit from pluripotency. Moreover, the data sets described here will facilitate future studies into the temporal control of transcriptional responses to inductive cues across multiple embryonic lineages.

Materials and methods

Embryological methods

Wild-type Xenopus laevis embryos were obtained using standard methods from a daily 2 pm fertilization from a single frog and placed into a 14C incubator at 2:45 pm until 8:30am. Ectodermal explants were manually dissected at early blastula (stage 8–9) from embryos cultured in 1X Marc’s Modified Ringer’s Solution (MMR) [0.1 M NaCL, 2 mM KCl, 1 mM MgSO₄, 2 mM CaCl₂, 5 mM HEPES (pH 7.8), 0.1 mM EDTA] from 8:30am-9:45am and then placed in a 20C incubator until 5 pm. Groups of 12–15 explants were collected at 9:45 am, 11am, 12:15 pm, 1:30 pm, 3 pm and 5 pm using sister embryos to confirm approximate stages of 9,10, 10.5, 11, 12, and 13 based on Nieuwkoop and Faber staging [74]. Explants for the neural progenitor lineage were generated using recombinant Noggin protein (R&D Systems) at a final concentration of 100 ng/mL in low calcium magnesium media supplemented with 0.1% bovine serum albumin (BSA) as a carrier for sequencing experiments, or using 20uM K02288 (Sigma) in 0.1X MMR for Westerns. The Noggin dose utilized was the minimum amount required to induce the neural progenitors Sox3 and Sox2 at stage 13, as determined by qPCR and in situ hybridization. The K02288 dose utilized was the minimum amount required to effectively block BMP signaling through stage 13, as determined by pSmad158 western blot. Endoderm lineage explants were generated using recombinant Activin protein (R&D Systems) at a final concentration of 160 ng/mL in 1XMMR supplemented with 0.1% BSA. This was the lowest dose that induced endoderm without mesoderm, as determined by qPCR with Bra(T) and Sox17 primers. Mesoderm lineage explants were generated using recombinant BMP4/7 heterodimer protein (R&D Systems) at a final concentration of 20 g/mL in 1XMMR supplemented with 0.1% BSA. This concentration of BMP4/7 was utilized as it induced levels of pSmad158 at stage 10 that were comparable to levels present at stage 10.5 in untreated caps. For morpholino experiments, a previously validated translation-blocking dand5 morpholino [116, 127]. (Gene Tools, Sequence: 'CTGGTGGCCTGGAACAACAGCATGT ') was injected in 4 cells at the eight-cell stage for a total of 40 pmol per embryo.

RNA isolation, cDNA synthesis, sequencing and qRT-PCR

RNA was isolated from blastula explants (12–15 explants) using Trizol (Life Technologies) followed by LiCl precipitation. 1 uG of purified RNA was used as a template for synthesizing cDNA using a High Capacity Reverse Transcription Kit (Life Technologies) Quantitative (q) RT-PCR was performed using SYBR Premix ExTaq 11 (Takara Bio) and detected using the Bio-Rad CFX96 Connect system. Brachyury primers used were Fwd: GAA GCG AAT GTT TCC AGT TC and Rev: ACA TAC TTC CAG CGG TGG TT. Expression was normalized to Ornitihine Decarboxylase (ODC) ODC primers used were Fwd: TGA AAA CAT GGG TGC CTA CA and Rev: TGC CAG TGT GGT CTT GAC AT. The fold change was calculated relative to stage nine samples from the same time course experiment. The results show the mean of three independent biological replicates, with error bars depicting the SEM. An unpaired, two-tailed t-test was utilized to determine significance. 500 ng of RNA was used for library prep with TruSeq mRNA library prep kit (Illumina) and sequenced using Next Seq 500 Sequencing (epi,neur, endo) or HiSeq4000(meso).

Western blot analysis

Blastula explants (20 explants/sample) were collected for specified stage / lineage and lysed in TNE phospho-lysis buffer [50 mM Tris–HCl (pH 7.4), 150 mM NaCL, 0.5 mM EDTA, and 0.5% Triton X-100, 2 mM Sodium Orthovanadate, 20 mM Sodium Fluoride, 10 mM B-Glycerophosphate, 1 MM Sodium Molybdate dihydrate] supplemented with protease inhibitors [Aprotinin, Leupeptin and phenylmethylsulfonyl fluoride (PMSF)] and a PhosStop phosphatase inhibitor and a complete Mini tablet (Roche). SDS-PAGE and western blot analyses were used to detect proteins and modifications using the following antibodies: anti-phospho SMAD2 (Ser465/467, Sigma, 1:500), Smad2 Polyclonal Antibody (Life Technologies 1:500) Anti-phospho Smad1/Smad5/Smad8 (Ser463/465, Sigma, 1:1000), Smad1 (D59D7, Cell Signal, 1:1000), Actin (A2066, Sigma-Aldrich, 1:5000). For chemiluminescence-based detection, horseradish peroxidase (HRP)-conjugated rabbit secondary antibodies were used (Vector Laboratories, 1:20,000). Blots were cut to the appropriate size based on the protein standard ladder before hybridization with antibodies in order to conserve antibody. Results shown are representative of at least three independent experiments. Raw images of blots in the main figures are included as Supplemental Fig. 8 and raw images of blots in supplemental figures are included as Supplemental Fig. 9.

RNA Seq processing and computational analysis

Read quality was evaluated using FastQC [128]. Mapping to X. laevis v9.2 genome downloaded from xenbase was performed using RSEM to get TPM values [129]. Alignment to X. laevis v9.2 genome was performed using STAR2.6.0 to get raw counts using standard parameters [130]. Computational Analysis of RNA Sequencing Data was performed using published R Packages. TPM data is an average of three biological replicates of combined data from S and L alleles for each gene, the width of each line represents SEM and graphs were plotted using ggplot2 [131]. A minimum raw read count of 15 was determined computationally using the voom function of the limma package. Based on the plot of the mean variance trend, filtering out genes with expression below 15 provided the best balance between filtering out lowly expressed genes without losing genes relevant to the transcriptome dynamics, thus all analyses are performed on genes with a minimum read count of 15 in any lineage at any stage [108].

Differential expression analysis was done between successive stages for each lineage and between pairs of lineages at corresponding states using DESeq2 with significance defined as p_adj ≤ 0.05 [132]. No fold change cut-off was used for examining temporal DESeq, in order to capture the greatest amount of dynamic change in the transcriptome over time. Overlapping DESeq genes were visualized using VennDiagram and UpsetR [133, 134]. GO and KEGG enrichment was calculated using the GOSeq R package [135]. Neural GO terms were identified based on their relevance to neural-related processes including, neuron, brain and nervous system development as well as axon formation and synapse firing. Mesoderm GO terms were identified based on their relevance to mesoderm-related processes including mesoderm formation, gastrulation and somitogenesis. DE genes for heatmaps were determined by ranking the genes DE between mesoderm and epidermis, and between endoderm and epidermis separately based on Log2 fold change and only genes with a minimum normalized expression of 10 TPM at the relevant stages were plotted using the pheatmap package [136], expression was depicted using z-scores. All lists of DE genes with padj and Log2FC values are provided as supplementary tables with transcription factors bolded, Gene lists from all Venn Diagrams and UpSet plots are also provided as supplementary tables. Hierarchical clustering was done using the pheatmap package, cluster rows = TRUE, cluster columns = TRUE.

Principal Component Analysis was done with the prcomp function in the default stats package, and visualized using ggfortify [137]. We assessed whether the observed patterns were statistically significant by testing whether the differences between samples of different types was significantly greater than the differences between replicates. We calculated pairwise distances between all pairs of samples in the PCA space and tested whether the distances for between- and within-type pairs differed using a nonparametric KS test. Overall, the between-type distances were far larger than the within-type distances (p = 6.27e-03). We also performed this same analysis within each pair of lineages (eg, testing whether the endo-epi distances were larger than the endo-endo and epi-epi distances) and found this pattern to be consistent though low numbers of sample pairs reduced the power of the KS test. The top 2 PCs were determined as the most significant based on the elbow of the PC plot and thus the only two that were plotted. Distances from stage 9 to 13 were calculated using the dist() function. Statistical significance of differences in distance were calculated using the Wilcoxon rank sum, wilcox.test() in R.

Pattern dynamics were determined for linear, quadratic and cubic patterns with limma analysis using the limma voom R package [108]. Potential mesendoderm GRN candidates were selected by examining the 10 genes with most similar differential quadratic and linear dynamics between epidermal and ventral mesoderm. Of these genes, those with minimum expression of 30 TPM, dynamics unique only to the endoderm and/or ventral mesoderm lineages and not already defined as mesendoderm GRN members were identified as possible novel GRN members. Genes graphed were those we were able to corroborate with other genomic studies.

WGCNA was performed using the WGCNA R package [138, 139]. Genes with a TPM > 15 in endoderm or mesoderm at stage 10, 10.5 or 11 were included in the analysis. Analysis was run on Log2 of the TPM. Power was determine using the sft function and TOMType = “signed”, minModuleSize = 30 and mergeCutHeight = 0.25.