- Open Access
Quantitative analysis of transcriptome dynamics provides novel insights into developmental state transitions
BMC Genomics volume 23, Article number: 723 (2022)
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.
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.
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.
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 . 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 . 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 . 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 . 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 . 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  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.
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 . 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 . 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 . 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  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.
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).
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  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.
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 . 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  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 . 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.
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 , 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].
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) . 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 . We compared our endoderm and ventral mesoderm trajectories to a recently published Activin-induced mesoderm time series  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 . 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).
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 . 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  and Ca14, Tmcc1, and Nptx2 were upregulated in stage 11 embryos in response to somatic cell nuclear transfer from an endoderm cell . Similarly, Pygm has been identified as a target of Myod . 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).
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 . 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 , 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.
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 .
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 , significant down-regulation of the pluripotency factor Foxi2 (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 . 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 , 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 , 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
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 MgSO4, 2 mM CaCl2, 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 . 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 . Mapping to X. laevis v9.2 genome downloaded from xenbase was performed using RSEM to get TPM values . Alignment to X. laevis v9.2 genome was performed using STAR2.6.0 to get raw counts using standard parameters . 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 . 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 .
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 padj ≤ 0.05 . 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 . 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 , 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 . 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 . 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.
Availability of data and materials
The datasets generated and analyzed during the current study are available in the NCBI repository and are accessible through GEO Series accession number GSE198598.
Bone morphogenetic protein
Bovine serum albumin
Central nervous system
Differential expression analysis for sequencing count data
Embryonic stem cells
Human embryonic stem cells
Fibroblast growth factor
Gene regulatory network
Horse radish peroxidase
Kyoto encyclopedia of genes and genomes
Mouse embryonic stem cells
Marc’s modified ringer’s solution
Principal component analysis
Real time quantitative reverse transcription polymerase chain reaction
Transforming growth factor beta
Transcripts per million
Waddington CH. The strategy of the genes: a discussion of some aspects of theoretical biology. London: Gorge Allen and Unwin; 1957.
Snape A, Wylie CC, Smith JC, Heasman J. Changes in states of commitment of single animal pole blastomeres of Xenopus laevis. Dev Biol. 1987;119(2):503–10.
Ariizumi T, Asashima M. In vitro induction systems for analyses of amphibian organogenesis and body patterning. Int J Dev Biol. 2003;45(1):273–9.
Ariizumi T, et al. Isolation and differentiation of Xenopus animal cap cells. Curr Protoc Stem Cell Biol. 2009;9(1):1D – 5.
Ariizumi T, Michiue T, Asashima M. In vitro induction of Xenopus embryonic organs using animal cap cells. Cold Spring Harbor Protocols. 2017;2017(12):pdb-prot097410.
Satou-Kobayashi Y, Kim JD, Fukamizu A, Asashima M. Temporal transcriptomic profiling reveals dynamic changes in gene expression of Xenopus animal cap upon Activin treatment. Sci Rep. 2021;11(1):1–12.
Spemann H, Mangold H. über Induktion von Embryonalanlagen durch Implantation artfremder Organisatoren. Archiv für mikroskopische Anatomie und Entwicklungsmechanik. 1924;100(3):599–638.
Harland R, Gerhart J. Formation and function of Spemann’s organizer. Annu Rev Cell Dev Biol. 1997;13(1):611–67.
Niehrs C. Regionally specific induction by the Spemann-Mangold organizer. Nat Rev Genet. 2004;5(6):425–34.
Whitfield TT, Heasman J, Wylie CC. Early embryonic expression of XLPOU-60, a Xenopus POU-domain protein. Dev Biol. 1995;169(2):759–69.
Henig C, Elias S, Frank D. A POU protein regulates mesodermal competence to FGF in Xenopus. Mech Dev. 1998;71(1–2):131–42.
Lef J, Clement JH, Oschwald R, Köster M, Knöchel W. Spatial and temporal transcription patterns of the forkhead related XFD-2/XFD-2′ genes in Xenopus laevis embryos. Mech Dev. 1994;45(2):117–26.
Paraiso KD, Cho JS, Yong J, Cho KWY. Early Xenopus gene regulatory programs, chromatin states, and the role of maternal transcription factors. Curr Top Dev Biol. 2020;139:35–60.
Wilson PA, Lagna G, Suzuki A, Hemmati-Brivanlou A. Concentration-dependent patterning of the Xenopus ectoderm by BMP4 and its signal transducer Smad1. Development. 1997;124(16):3177–84.
Wilson PA, Hemmati-Brivanlou A. Induction of epidermis and inhibition of neural fate by Bmp-4. Nature. 1995;376(6538):331–3.
Suzuki A, Kaneko E, Ueno N, Hemmati-Brivanlou A. Regulation of epidermal induction by BMP2 and BMP7 signaling. Dev Biol. 1997;189(1):112–22.
Little SC, Mullins MC. Bone morphogenetic protein heterodimers assemble heteromeric type I receptor complexes to pattern the dorsoventral axis. Nat Cell Biol. 2009;11(5):637–43.
Kretzschmar M, Liu F, Hata A, Doody J, Massagué J. The TGF-beta family mediator Smad1 is phosphorylated directly and activated functionally by the BMP receptor kinase. Genes Dev. 1997;11(8):984–95.
Macı́as-Silva, M., Hoodless, P. A., Tang, S. J., Buchwald, M., & Wrana, J. L. Specific activation of Smad1 signaling pathways by the BMP7 type I receptor, ALK2. J Biol Chem. 1998;273(40):25628–36.
Shi Y, Massagué J. Mechanisms of TGF-β signaling from cell membrane to the nucleus. cell. 2003;113(6):685–700.
Jonas E, Sargent TD, Dawid IB. Epidermal keratin gene expressed in embryos of Xenopus laevis. Proc Natl Acad Sci. 1985;82(16):5413–7.
Dirksen ML, Morasso MI, Sargent TD, Jamrich M. Differential expression of a Distal-less homeobox gene Xdll-2 in ectodermal cell lineages. Mech Dev. 1994;46(1):63–70.
Luo T, Matsuo-Takasaki M, Sargent TD. Distinct roles for Distal-less genes Dlx3 and Dlx5 in regulating ectodermal development in Xenopus. Molecular Reproduction and Development: Incorporating Gamete Research. 2001;60(3):331–7.
Hemmati-Brivanlou A, Melton DA. Inhibition of Activin receptor signaling promotes neuralization in Xenopus. Cell. 1994;77(2):273–81.
Hemmati-Brivanlou A, Melton D. Vertebrate embryonic cells will become nerve cells unless told otherwise. Cell. 1997;88(1):13–7.
Muñoz-Sanjuán I, Brivanlou AH. Neural induction, the default model and embryonic stem cells. Nat Rev Neurosci. 2002;3(4):271–80.
Lamb TM, et al. Neural induction by the secreted polypeptide Noggin. Science. 1993;262(5134):713–8.
Re’em-Kalma Y, Lamb T, Frank D. Competition between Noggin and bone morphogenetic protein 4 activities may regulate dorsalization during Xenopus development. Proc Natl Acad Sci. 1995;92(26):12141–5.
Sasai Y, Lu B, Steinbeisser H, De Robertis EM. Regulation of neural induction by the Chd and Bmp-4 antagonistic patterning signals in Xenopus. Nature. 1995;376(6538):333–6.
Piccolo S, Sasai Y, Lu B, De Robertis EM. Dorsoventral patterning in Xenopus: inhibition of ventral signals by direct binding of chordin to BMP-4. Cell. 1996;86(4):589–98.
Piccolo S, et al. The head inducer Cerberus is a multifunctional antagonist of Nodal. BMP and Wnt signals Nature. 1999;397(6721):707–10.
Iemura, et al. Direct binding of follistatin to a complex of bone-morphogenetic protein and its receptor inhibits ventral and epidermal cell fates in early Xenopus embryo. Proc Natl Acad Sci. 1998;95(16):9337–42.
Grunz H, Tacke L. Neural differentiation of Xenopus laevis ectoderm takes place after disaggregation and delayed reaggregation without inducer. Cell Differ Dev. 1989;28(3):211–7.
Streit A, Lee, et al. Chordin regulates primitive streak development and the stability of induced neural cells, but is not sufficient for neural induction in the chick embryo. Development. 1998;125(3):507–19.
Streit A, Stern CD. Mesoderm patterning and somite formation during node regression: differential effects of chordin and Noggin. Mech Dev. 1999;85(1–2):85–96.
Streit A, Stern CD. Establishment and maintenance of the border of the neural plate in the chick: involvement of FGF and BMP activity. Mech Dev. 1999;82(1–2):51–66.
Stern CD. Neural induction: old problem, new findings, yet more questions. Development. 2005;132(9):2007–21.
Linker C, Stern CD. Neural induction requires BMP inhibition only as a late step, and involves signals other than FGF and Wnt antagonists. Development. 2004;131(22):5671–81.
Delaune E, Lemaire P, Kodjabachian L. Neural induction in Xenopus requires early FGF signalling in addition to BMP inhibition. Development (Cambridge, England). 2005;132(2):299–310.
Wawersik S, Evola C, Whitman M. Conditional BMP inhibition in Xenopus reveals stage-specific roles for BMPs in neural and neural crest induction. Dev Biol. 2005;277(2):425–42.
Wills AE, Choi VM, Bennett MJ, Khokha MK, Harland RM. BMP antagonists and FGF signaling contribute to different domains of the neural plate in Xenopus. Dev Biol. 2010;337(2):335–50.
Collignon J, et al. A comparison of the properties of Sox-3 with Sry and two related genes, Sox-1 and Sox-2. Development. 1996;122(2):509–20.
Mizuseki K, Kishi M, Matsui M, Nakanishi S, Sasai Y. Xenopus Zic-related-1 and Sox-2, two factors induced by chordin, have distinct activities in the initiation of neural induction. Development. 1998;125(4):579–87.
Penzel R, Oschwald R, Chen Y, Tacke L, Grunz H. Characterization and early embryonic expression of a neural specific transcription factor xSOX3 in Xenopus laevis. Int J Dev Biol. 2003;41(5):667–77.
Pannese M, et al. The Xenopus homologue of Otx2 is a maternal homeobox gene that demarcates and specifies anterior body regions. Development. 1995;121(3):707–20.
Kablar B, et al. Xotx genes in the developing brain of Xenopus laevis. Mech Dev. 1996;55(2):145–58.
Andreazzoli M, Pannese M, Boncinelli E. Activating and repressing signals in head development: the role of Xotx1 and Xotx2. Development. 1997;124(9):1733–43.
Plouhinec JL, et al. A molecular atlas of the developing ectoderm defines neural, neural crest, placode, and nonneural progenitor identity in vertebrates. PLoS Biol. 2017;15(10):e2004045.
Asashima M, et al. Mesodermal induction in early amphibian embryos by Activin A (erythroid differentiation factor). Rouxs Arch Dev Biol. 1990;198(6):330–5.
Asashima M, et al. The vegetalizing factor belongs to a family of mesoderm-inducing proteins related to erythroid differentiation factor. Naturwissenschaften. 1990;77(8):389–91.
Henry GL, Brivanlou IH, Kessler DS, Hemmati-Brivanlou A, Melton DA. TGF-beta signals and a pattern in Xenopus laevis endodermal development. Development. 1996;122(3):1007–15.
Smith JC, Price BMJ, Nimmen KV, Huylebroeck D. Identification of a potent Xenopus mesoderm-inducing factor as a homologue of Activin A. Nature. 1990;345(6277):729–31.
Hemmati-Brivanlou A, Melton DA. A truncated Activin receptor inhibits mesoderm induction and formation of axial structures in Xenopus embryos. Nature. 1992;359(6396):609–14.
Gurdon JB, Harger P, Mitchell A, Lemaire P. Activin signalling and response to a morphogen gradient. Nature. 1994;371(6497):487–92.
McDowell N, Zorn AM, Crease DJ, Gurdon JB. Activin has direct long-range signalling activity and can form a concentration gradient by diffusion. Curr Biol. 1997;7(9):671–81.
McDowell N, Gurdon JB. Activin as a morphogen in Xenopus mesoderm induction. Semin Cell Dev Biol. 1999;10(3):311–7.
Gamer LW, Wright CV. Autonomous endodermal determination in Xenopus: regulation of expression of the pancreatic gene XlHbox 8. Dev Biol. 1995;171(1):240–51.
Agius E, Oelgeschlager M, Wessely O, Kemp C, De Robertis EM. Endodermal Nodal-related signals and mesoderm induction in Xenopus. Development. 2000;127(6):1173–83.
Dale L, Matthews G, Colman A. Secretion and mesoderm-inducing activity of the TGF-beta-related domain of Xenopus Vg1. EMBO J. 1993;12(12):4471–80.
Thomsen GH, Melton DA. Processed Vg1 protein is an axial mesoderm inducer in Xenopus. Cell. 1993;74(3):433–41.
Kessler DS, Melton DA. Induction of dorsal mesoderm by soluble, mature Vg1 protein. Development. 1995;121(7):2155–64.
Jones CM, Kuehn MR, Hogan BL, Smith JC, Wright CV. Nodal-related signals induce axial mesoderm and dorsalize mesoderm during gastrulation. Development. 1995;121(11):3651–62.
Reissmann, et al. The orphan receptor ALK7 and the Activin receptor ALK4 mediate signaling by Nodal proteins during vertebrate development. Genes Dev. 2001;15(15):2010–22.
Kumar, et al. Nodal signaling uses Activin and transforming growth factor-β receptor-regulated Smads. J Biol Chem. 2001;276(1):656–61.
Hudson C, Clements D, Friday RV, Stott D, Woodland HR. Xsox17α and-β mediate endoderm formation in Xenopus. Cell. 1997;91(3):397–405.
Sasai Y, Lu B, Piccolo S, De Robertis EM. Endoderm induction by the organizer-secreted factors chordin and Noggin in Xenopus animal caps. EMBO J. 1996;15(17):4547–55.
Sokol SY, Melton DA. Interaction of Wnt and Activin in dorsal mesoderm induction in Xenopus. Dev Biol. 1992;154(2):348–55.
Green JB, New HV, Smith JC. Responses of embryonic Xenopus cells to Activin and FGF are separated by multiple dose thresholds and correspond to distinct axes of the mesoderm. Cell. 1992;71(5):731–9.
Nishimatsu SI, Thomsen GH. Ventral mesoderm induction and patterning by bone morphogenetic protein heterodimers in Xenopus embryos. Mech Dev. 1998;74(1–2):75–88.
Hemmati-Brivanlou A, Thomsen GH. Ventral mesodermal patterning in Xenopus embryos: expression patterns and activities of BMP-2 and BMP-4. Dev Genet. 1995;17(1):78–89.
Graff JM, Thies RS, Song JJ, Celeste AJ, Melton DA. Studies with a Xenopus BMP receptor suggest that ventral mesoderm-inducing signals override dorsal signals in vivo. Cell. 1994;79(1):169–79.
Owens ND, et al. Measuring absolute RNA copy numbers at high temporal resolution reveals transcriptome kinetics in development. Cell Rep. 2016;14(3):632–47.
Session AM, et al. Genome evolution in the allotetraploid frog Xenopus laevis. Nature. 2016;538(7625):336–43.
Nieuwkoop PD, Faber J, editors. Normal Table of Xenopus laevis (Daudin). New York: Garland Publishing; 1994.
Prasad MS, Sauka-Spengler T, LaBonne C. Induction of the neural crest state: control of stem cell attributes by gene regulatory, post-transcriptional and epigenetic interactions. Dev Biol. 2012;366(1):10–21.
Köster M, et al. Bone morphogenetic protein 4 (BMP-4), a member of the TGF-β family, in early embryos of Xenopus laevis: analysis of mesoderm inducing activity. Mech Dev. 1991;33(3):191–9.
Collart C, et al. High-resolution analysis of gene activity during the Xenopus mid-blastula transition. Development. 2014;141(9):1927–39.
Gentsch GE, Spruce T, Owens ND, Smith JC. Maternal pluripotency factors initiate extensive chromatin remodelling to predefine first response to inductive signals. Nat Commun. 2019;10(1):1–22.
Paranjpe SS, Jacobi UG, van Heeringen SJ, Veenstra GJC. A genome-wide survey of maternal and embryonic transcripts during Xenopus tropicalis development. BMC Genomics. 2013;14(1):1–17.
Mir A, et al. FoxI1e activates ectoderm formation and controls cell position in the Xenopus blastula. Development. 2007;134(4):779–88.
Tao J, et al. BMP4-dependent expression of Xenopus Grainyhead-like 1 is essential for epidermal differentiation. Development. 2005;132(5):1021–34.
Cha SW, McAdams M, Kormish J, Wylie C, Kofron M. Foxi2 is an animally localized maternal mRNA in Xenopus, and an activator of the zygotic ectoderm activator Foxi1e. PLoS ONE. 2012;7(7):e41782.
Hyodo-Miura J, et al. Involvement of NLK and Sox11 in neural induction in Xenopus development. Genes Cells. 2002;7(5):487–96.
Nakata K, Nagai T, Aruga J, Mikoshiba K. Xenopus Zic family and its role in neural and neural crest development. Mech Dev. 1998;75(1–2):43–51.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.
Virtanen P, et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature methods. 2020;17(3):261–72.
Onichtchouk D, et al. The Xvent-2 homeobox gene is part of the BMP-4 signalling pathway controlling [correction of controling] dorsoventral patterning of Xenopus mesoderm. Development. 1996;122(10):3045–53.
Hollnagel A, Oehlmann V, Heymer J, Rüther U, Nordheim A. Id genes are direct targets of bone morphogenetic protein induction in embryonic stem cells. J Biol Chem. 1999;274(28):19838–45.
Chen X, Rubock MJ, Whitman M. A transcriptional partner for MAD proteins in TGF-β signalling. Nature. 1996;383(6602):691–6.
Ryan K, Garrett N, Mitchell A, Gurdon JB. Eomesodermin, a key early gene in Xenopus mesoderm differentiation. Cell. 1996;87(6):989–1000.
Henry GL, Melton DA. Mixer, a homeobox gene required for endoderm development. Science. 1998;281(5373):91–6.
Rashbass J, Taylor MV, Gurdon JB. The DNA-binding protein E12 co-operates with XMyoD in the activation of muscle-specific gene expression in Xenopus embryos. EMBO J. 1992;11(8):2981–90.
Wills AE, Baker JC. E2a is necessary for Smad2/3-dependent transcription and the direct repression of lefty during gastrulation. Dev Cell. 2015;32(3):345–57.
Cordenonsi M, et al. Links between tumor suppressors: p53 is required for TGF-β gene responses by cooperating with Smads. Cell. 2003;113(3):301–14.
Montague TG, Schier AF. Vg1-Nodal heterodimers are the endogenous inducers of mesendoderm. Elife. 2017;6:e28183.
Suzuki A, Kaneko E, Maeda J, Ueno N. Mesoderm induction by BMP-4 and-7 heterodimers. Biochem Biophys Res Commun. 1997;232(1):153–6.
Massagué J, Wotton D. Transcriptional control by the TGF-β/Smad signaling system. EMBO J. 2000;19(8):1745–54.
Wardle FC, Smith JC. Transcriptional regulation of mesendoderm formation in Xenopus. Semin Cell Dev Biol. 2006;17(1):99–109.
Suzuki A, Ueno N, Hemmati-Brivanlou A. Xenopus msx1 mediates epidermal induction and neural inhibition by BMP4. Development. 1997;124(16):3037–44.
Rastegar S, Friedle H, Frommer G, Knöchel W. Transcriptional regulation of Xvent homeobox genes. Mech Dev. 1999;81(1–2):139–49.
Smith JC, Price BMJ, Green JBA, Weigel D, Herrmann BG. Expression of a Xenopus homolog of Brachyury (T) is an immediate-early response to mesoderm induction. Cell. 1991;67(1):79–87.
Tada M, Casey ES, Fairclough L, Smith JC. Bix1, a direct target of Xenopus T-box genes, causes formation of ventral mesoderm and endoderm. Development. 1998;125(20):3997–4006.
Rosa FM. Mix.1, a homeobox mRNA inducible by mesoderm inducers, is expressed mostly in the presumptive endodermal cells of Xenopus embryos. Cell. 1989;57(6):965–74.
Schmidt JE, Suzuki A, Ueno N, Kimelman D. Localized BMP-4 mediates dorsal/ventral patterning in the early Xenopus embryo. Dev Biol. 1995;169(1):37–50.
Ding Y, et al. Genome-wide analysis of dorsal and ventral transcriptomes of the Xenopus laevis gastrula. Dev Biol. 2017;426(2):176–87.
Charney RM, Paraiso KD, Blitz IL, Cho KWY. A gene regulatory program controlling early Xenopus mesendoderm formation: network conservation and motifs. Semin Cell Dev Biol. 2017;66:12–24.
Jansen C, et al. Uncovering the mesendoderm gene regulatory network through multi-omic data integration. Cell Rep. 2022;38(7):110364.
Law CW, Chen Y, Shi W, Smyth GK. voom: Precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biol. 2014;15(2):1–17.
Ding Y, et al. Bighead is a Wnt antagonist secreted by the Xenopus Spemann organizer that promotes Lrp6 endocytosis. Proc Natl Acad Sci. 2018;115(39):E9135–44.
Hörmanseder E, et al. H3K4 methylation-dependent memory of somatic cell identity inhibits reprogramming and development of nuclear transfer embryos. Cell Stem Cell. 2017;21(1):135–43.
McQueen C, Pownall ME. An analysis of MyoD-dependent transcription using CRISPR/Cas9 gene targeting in Xenopus tropicalis embryos. Mech Dev. 2017;146:1–9.
Sato SM, Sargent TD. Localized and inducible expression of Xenopus-posterior (Xpo), a novel gene active in early frog embryos, encoding a protein with a ‘CCHC’finger domain. Development. 1991;112(3):747–53.
Watanabe M, Whitman M. FAST-1 is a key maternal effector of mesoderm inducers in the early Xenopus embryo. Development. 1999;126(24):5621–34.
Christian JL, McMAHON JA, McMAHON AP, Moon RT. Xwnt-8, a Xenopus Wnt-1/int-1-related gene responsive to mesoderm-inducing growth factors, may play a role in ventral mesodermal patterning during embryogenesis. Development. 1991;111(4):1045–55.
Bell E, Muñoz-Sanjuán I, Altmann CR, Vonica A, Brivanlou AH. Cell fate specification and competence by Coco, a maternal BMP. TGFβand Wnt inhibitor Development. 2003;130(7):1381–9.
Vonica A, Brivanlou AH. The left–right axis is regulated by the interplay of Coco, Xnr1 and derrière in Xenopus embryos. Dev Biol. 2007;303(1):281–94.
Bates TJ, Vonica A, Heasman J, Brivanlou AH, Bell E. Coco regulates dorsoventral specification of germ layers via inhibition of TGFβ signalling. Development. 2013;140(20):4177–81.
Reich S, Weinstein DC. Repression of inappropriate gene expression in the vertebrate embryonic ectoderm. Genes. 2019;10(11):895.
Kirmizitas A, Gillis WQ, Zhu H, Thomsen GH. Gtpbp2 is required for BMP signaling and mesoderm patterning in Xenopus embryos. Dev Biol. 2014;392(2):358–67.
Streit A, Berliner AJ, Papanayotou C, Sirulnik A, Stern CD. Initiation of neural induction by FGF signalling before gastrulation. Nature. 2000;406(6791):74–8.
Smukler SR, Runciman SB, Xu S, van der Kooy D. Embryonic stem cells assume a primitive neural stem cell fate in the absence of extrinsic influences. J Cell Biol. 2006;172(1):79–90.
Vallier L, Reynolds D, Pedersen RA. Nodal inhibits differentiation of human embryonic stem cells along the neuroectodermal default pathway. Dev Biol. 2004;275(2):403–21.
Briggs JA, et al. The dynamics of gene expression in vertebrate embryogenesis at single-cell resolution. Science. 2018;360(6392):eaar5780.
Pereira LA, Wong MS, Mei Lim S, Stanley EG, Elefanty AG. The Mix family of homeobox genes--key regulators of mesendoderm formation during vertebrate development. Dev Biol. 2012;367(2):163–77.
Mead PE, Brivanlou IH, Kelley CM, Zon LI. BMP-4-responsive regulation of dorsal-ventral patterning by the homeobox protein Mix.1. Nature. 1996;382(6589):357–60.
Wessely O, Agius E, Oelgeschläger M, Pera EM, De Robertis EM. Neural induction in the absence of mesoderm: β-catenin-dependent expression of secreted BMP antagonists at the blastula stage in Xenopus. Dev Biol. 2001;234(1):161–73.
Maerker M, et al. Bicc1 and Dicer regulate left-right patterning through post-transcriptional control of the Nodal inhibitor Dand5. Nat Commun. 2021;12(1):1–15.
Andrews S. FastQC: a Quality Control Tool for High Throughput Sequence Data. 2010. Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12(1):1–16.
Dobin A, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21.
Wickham H. ggplot2: Elegant Graphics for Data Analysis. New York: Springer-Verlag; 2016. https://ggplot2.tidyverse.org. ISBN 978-3-319-24277-4.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):1–21.
Chen H, Boutros PC. VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinformatics. 2011;12(1):1–7.
Conway JR, Lex A, Gehlenborg N. UpSetR: an R package for the visualization of intersecting sets and their properties. Bioinformatics. 2017;33(18):2938–40.
Young MD, Wakefield MJ, Smyth GK, Oshlack A. goseq: Gene Ontology testing for RNA-seq datasets. R Bioconductor. 2012;8:1–25.
Kolde R, Kolde MR. Package ‘pheatmap.’ R package. 2015;1(7):790.
Tang Y, Horikoshi M, Li W. ggfortify: unified interface to visualize statistical results of popular R packages. R J. 2016;8(2):474.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9(1):1–13.
Langfelder P, Horvath S. Fast R Functions for Robust Correlations and Hierarchical Clustering. J Stat Softw. 2012;46(11):i11.
We thank Joe Nguyen and Rachel Schermerhorn for invaluable technical assistance and the members of the LaBonne laboratory, Madhav Mani, and Eliza Duvall for helpful discussions.
This work was supported by grants from the Simons Foundation/SFARI (597491-RWC) the National Science Foundation (1764421) and the NIH (R01GM116538).
Ethics approval and consent to participate
All animal procedures were approved by Northwestern University’s Institutional Animal Care and Use Committee, are in accordance with the National Institutes of Health’s Guide for the care and use of Laboratory Animals and are in accordance with ARRIVE guidelines. Wild-type Xenopus laevis embryos were utilized at embryonic stages that are not considered vertebrate animals for regulatory purposes. At these stages it is not possible to determine the biological sex of the embryos. Adult Xenopus laevis were used to obtain embryos only and were not subject to experimental procedures.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplemental Figure 1. Xenopus blastula explants can be reprogrammed to adopt any lineage. (A-E) RNA Seq TPM expression over time of (A) average of all monotonically decreasing genes with a >25fold decrease between stage 9 and 13 and maximum TPM of 20 at stage 13, (B) epidermal marker Dlx3, (C) neural marker Otx2, (D) mesoderm marker Evx1 (E) endoderm marker Endodermin. Graphs are sums of S+L allele. Width of lines represents SEM of three biological replicates.
Supplemental Figure 2. Transcription Factor UpSet Plots and Zygotic PCA. (A-E) UpSet plots of transcription factors expressed at a minimum of 10 TPM at (A) stage 10, (B) stage 10.5, (C) stage 11, (D) stage 12, (E) stage 13. X-axes show genes unique to each lineage and overlapping in all different combinations of lineages ordered from largest number of genes to smallest. (F) PCA on only zygotically expressed genes
Temporal DESeq UpSet Plots. (A-J) UpSet plots of Differentially Expressed Genes (A) increased between stages 9 and 10, (B) decreased between stages 9 and 10, (C) increased between stages 10 and 10.5, (D) decreased between stages 10 and 10.5, (E) increased between stages 10.5 and 11, (F) decreased between stages 10.5 and 11, (G) increased between stages 11 and 12, (H) decreased between stages 11 and 12, (I) increased between stages 12 and 13, (J) decreased between stages 12 and 13. X-axes show genes unique to each lineage and overlapping in all different combinations of lineages ordered from largest number of genes to smallest.
Epidermal and Neural Lineages Diverge at St 10.5 (A) Number of enriched Neural GO Terms in genes significantly higher in the neural lineage (blue) and epidermal lineage (red) at each developmental stage. (B) Heatmap of the top 20 DE genes by Log2FC with a minimum expression of 10TPM between the epidermal and neural lineages at stage 10.5 (C) KEGG enrichment analysis of genes differentially expressed between epidermal and neural lineage at each developmental stage. Genes significantly increased in the epidermal lineage are enriched for TGF-beta genes, as defined by KEGG database from stages 10.5-13.
BMP Ventralizes Mesoderm. (A)Western blot analysis of lysates of developing mesoderm (20ng/uL BMP4/7) and endoderm (160ng/uL Activin) explants for pSmad2 and Smad2 with Actin loading control. Gels were cropped to minimize empty space, uncropped gels are Supplemental Figure 9A. (B) PCA of published Activin induced mesoderm data (Satou-Kobayashi et al. 2021) with our ventral mesoderm and endoderm data shows that activin induced mesoderm clusters with our activin induced endoderm rather than BMP4/7 induced mesoderm. (C) Hierarchical clustering of Dorsal/Ventral genes (Ding et al. 2017) shows that endoderm and ventral mesoderm cluster on opposite sides of the heat map, suggesting BMP4/7 effectively ventralizes the mesoderm.
Mesendoderm Analysis. (A) Number of differentially expressed genes between the endoderm and ventral mesoderm lineages at each developmental stage (padj ≤ 0.05). (B) KEGG enrichment analysis of genes differentially expressed between ventral mesoderm and endoderm lineage at each developmental stage. Genes significantly increased in the endoderm lineage are enriched for TGF-beta genes, as defined by KEGG database from stages 10-10.5 and genes significantly higher in the ventral mesoderm lineage are enriched for TGFbeta genes for stages 11-12. (C) WGCNA on stages 10,10.5 and 11 in the mesoderm and endoderm lineages identifies 22 gene modules, the blue module demonstrating increasing correlation over time to the endoderm lineage and the brown module demonstrating increasing correlation over time to the mesoderm lineage.
Early not Exogenous BMP Causes Mesoderm Formation. (A) Number of differentially expressed genes between the epidermal and ventral mesoderm lineages at each developmental stage (padj ≤ 0.05). (B) Number of enriched Mesoderm GO Terms (solid line) and Epidermis GO Terms (dashed line) in genes significantly higher in the ventral mesoderm lineage (purple) and epidermal lineage (red) at each developmental stage. (C) Western Blot Analysis of lysates for epidermal (WT) and BMP4/7 treated at stage 9 (BMP4/7 20ng/uL) explants collected at stage 10 and epidermal (WT) and BMP4/7 treated at stage 10.5 (BMP4/7 20ng/uL) explants collected at stage 11 for pSmad1/5/8 with Actin loading control. Gel was cropped to show one replicate, uncropped gel is Supplemental Figure 9B. (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), treated with BMP4/7(20ng/uL) at stage 9 and treated with BMP4/7 (20ng/uL) at stage 10.5 (***P<0.005).
Raw Images of Main Figure Western Blots. (A) Figure 4A Western blot analysis of lysates of developing epidermal (WT) and neural (20uM K02288) explants for pSmad1/5/8 and Smad1 with Actin loading control. Blot was cut prior to antibody hybridization to conserve antibody, unedited scan showing blot edges provided. (B) Figure 7A Western Blot analysis of lysates for developing epidermal (WT) and ventral mesoderm (BMP4/7 20ng/uL) explants for pSmad1/5/8 and Smad1 with Actin loading control. Blot was cut prior to antibody hybridization to conserve antibody. Initial scan showing blot edges (top) and high resolution unedited scan of relevant size (bottom) both provided. (C) Figure 8C Western Blot of lysates for stage 10 epidermal (WT), Dand5 MO injected (40pmol/embryo), ventral mesoderm (BMP4/7 20ng/uL), and endoderm (160ng/uL Activin) explants for pSmad1/5/8 with Actin loading control. Blot was cut prior to antibody hybridization to conserve antibody. Initial scans showing blot edges and high resolution unedited scans of relevant size for all three replicates are provided. Figure 8C is replicate 1.
Raw Images of Supplemental Western Blots. (A) Supplemental Figure 5A Western blot analysis of lysates of developing mesoderm (20ng/uL BMP4/7) and endoderm (160ng/uL Activin) explants for pSmad2 and Smad2 with Actin loading control. Blot was cut prior to antibody hybridization to conserve antibody. Initial scan showing blot edges (top) and unedited high resolution scan of relevant size bands (bottom) both provided. (B) Supplemental Figure 7C Western Blot Analysis of lysates for epidermal (WT) and BMP4/7 treated at stage 9 (BMP4/7 20ng/uL) explants collected at stage 10 and epidermal (WT) and BMP4/7 treated at stage 10.5 (BMP4/7 20ng/uL) explants collected at stage 11 for pSmad1/5/8 with Actin loading control. Blot was cut prior to hybridization with antibody to conserve antibody. Initial scan showing blot edges (top) and unedited high resolution scan of revelant size bands (bottom) both provided. All three replicates are shown, Supplemental Figure 6A is replicate 1.
List of all Monotonically Decreasing Genes with >25 Fold Decrease from Stage 9 to 13 and their average TPM across 3 biological replicates in all conditions.
Gene Lists from all Overlapping Categories from UpSet Plots of all Transcription Factors expressed >10 TPM corresponding with Supplemental Figure 2.
Transcription Factor WGCNA Epidermal Modules.
List of Genes Differentially Expressed in Temporal DESeq for each lineage with Log2 Fold Change and pvalues, corresponding with Figure 2E.
List of Genes for all Overlapping Temporal DESeq Groups, corresponding with supplemental figure 3.
Temporal DESeq Epi Neur Venn Diagram Gene Lists corresponding with Figure 3B.
List of all Neural GO Terms and their over represented p value from GO Analysis on genes differentially expressed between epidermal and neural lineages corresponding with supplemental figure 4A.
List of all genes differnetially expressed between epidermal and neural lineages and their Log2 Fold Change and pvalue corresponding with Fig. 3D, genes differentially expressed between endoderm and and mesoderm lineages and their Log2 Fold Change and pvalue corresponding with Fig. 6A, and genes differentially expressed between epidermis and and mesoderm lineages and their Log2 Fold Change and pvalue corresponding with Fig. 7A.
List of all genes from the venn diagrams of genes temporally differentially expressed in endoderm and mesoderm lineages corresponding with Fig. 5B.
Endo-Meso St10-11 WGCNA Module Gene Lists.
List of Genes from Venn Diagrams of Genes Temporally Differentially Expressed in mesoderm and epidermal lineages corresponding with Fig. 7C.
About this article
Cite this article
Johnson, K., Freedman, S., Braun, R. et al. Quantitative analysis of transcriptome dynamics provides novel insights into developmental state transitions. BMC Genomics 23, 723 (2022). https://doi.org/10.1186/s12864-022-08953-3
- Neural default model