Introduction

During gastrulation, instructed by genetic and epigenetic signals, naïve cells from blastula progressively acquire distinct cell fates and are allocated to specific domains within the embryo^1,2. Nodal, a member of the transforming growth factor β (TGF-β) superfamily, acts as a morphogen^1,3–5 and plays highly conserved roles in mesendoderm induction and patterning^6–8. It is widely accepted that Nodal induces endodermal and mesodermal cell fates via its concentration gradient^8–10. Higher Nodal signaling activity bias cells toward endoderm, while lower levels are thought to promote mesoderm induction^11,12. Mesodermal and endodermal common progenitors are induced by Nodal and mixed together at the onset of gastrulation¹, while it remains elusive how the common progenitors of mesendoderm took on distinct lineages. Several recent studies reported that multiple important signaling pathways interacted with Nodal to drive mesendoderm separation^2,13–15. Notably, lateral endoderm specification may comply with a stochastic cell fate switching model in zebrafish¹⁶.

The specification of dorsal mesendoderm requires the highest level of Nodal signaling, as evidenced by that this cell population was first to be depleted when Nodal signaling was gradually inhibited¹¹. Both the PP and anterior Endo originate from anterior mesendoderm (derived from the dorsal mesendoderm), and the commitment of these cell fates is regulated by similar levels of Nodal signaling. However, the mechanisms governing the separation of cell fates within the anterior mesendoderm remain unclear. Previous studies show that Nodal signaling enhances the duration of cell-cell contact in PP cells, creating a positive feedback loop between cell-cell contact duration and cell fate specification¹⁷. However, what factors downstream of Nodal pathway play major roles in segregating these two cell fates, when and how PP and Endo cells can be distinguished transcriptionally are the questions remaining to be investigated.

Previous research has shown that Gsc, whose expression can be induced by Nodal, acts as a repressor to suppress anterior endoderm specification by binding to the promoter region of sox17¹⁸. However, investigations into loss-of-function of Gsc in zebrafish have not yielded conclusive evidence supporting its significant roles in mesoderm and/or endoderm specification^19,20. These results imply that there might be additional, unidentified transcriptional repressors that work in collaboration with Gsc and possess a redundant function in suppressing the specification of endoderm.

In this study, we constructed a single-cell transcriptional trajectory to delineate the cell state segregation between PP and anterior Endo in zebrafish embryos and Nodal-injected zebrafish explants^21,22. Single-cell transcriptional trajectory and live imaging analyses demonstrated that anterior Endo cells originated and were specified from PP progenitors. Interestingly, this lineage separation can be traced back to as early as the onset of gastrulation. To achieve a more comprehensive understanding, we integrated single-cell datasets from wild-type, ndr1 (Nodal ligand) knockdown and lft1 (Nodal inhibitor) knockout Nodal explants at shield stage to construct an extended single-cell trajectory, which fully captured the branching event between PP and anterior Endo fates and unveiled a potential involvement of epigenetic regulation in this process. This finding was further supported by the experimental evidence that perturbation of the chromatin remodeler srcap could affect the cell specification of PP and anterior Endo. Furthermore, combining the multi-omic analysis and gain-/loss-of-function experiments, we found that a transcriptional repressor, Ripply1, cooperated with Gsc to supresse endoderm specification by directly binding to the regulatory elements of sox17 and sox32.

Results

Anterior Endo originates from PP progenitors in zebrafish

Both anterior Endo and PP are derived from embryonic shield and require high Nodal signaling to be specified^11,18,23. To investigate how these two cell lineages are further segregated, we collected the single-cell RNA-seq datasets from the published studies^24,25, identified the subpopulations representing the anterior Endo and PP cells, and reconstructed their lineage branching tree (Figures 1A-1C and S1A). We found that these two lineages started to be separated at shield stage, and their common progenitors highly expressed PP marker genes, like gsc, frzb, etc. (Figures 1C and S1A). Then, we further explored the dynamic molecular signatures of these two lineages during their fate determination. Pseudotime analysis revealed that PP and Endo exhibited distinct differentiation states by the end of gastrulation (Figure S1B). Of note, compared to sox32, gsc displayed earlier transcriptional differences in these two lineages (Figure S1C). The Uniform Manifold Approximation and Projection (UMAP) analysis²⁶ and PCA analysis showed that the common progenitors of PP and Endo clustered closely to PP and away from Endo (Figures S1D and S1E), suggesting that anterior Endo lineage is branched out from PP lineage. The cell lineage prediction analysis²⁷ also revealed that the progenitor cells resembled more of the PP lineage (Figures 1D, 1E and S1F), and similar result was obtained using another single-cell RNA-seq dataset (Figures S1G-S1I). Therefore, both UMAP and the cell lineage prediction analysis suggest that anterior Endo and PP originate from the same progenitor with transcriptome profile more similar to PP lineage. Anterior Endo becomes transcriptionally distinguishable from PP as early as at the shield stage.

Figure 1.

Cells from hatching gland (derived from PP) and pharyngeal pouch (derived from anterior Endo) are distant on the transcriptome-based phylogenetic tree of zebrafish embryos²⁵. Thus, the resolution of this dataset may not be sufficient to identify the exact lineage separation point between PP and Endo. To investigate the key molecules that regulate the PP and Endo separation, we employed our previously published Nodal induced embryonic explant system (referred to Nodal explant hereafter), which encompasses the axial mesendoderm, along with the PP and anterior Endo lineages²². Notably, the single-cell trajectory tree of Nodal explants showed a closer transcriptomic relationship between the PP and anterior Endo (Figure 1F). Thus, we reconciled the single-cell transcriptome dynamics to look for the branching point of PP and anterior Endo in the Nodal explants (Figures 1G, 1H, S2A and S2B). Consistent with the observations in wild-type embryos, the common progenitors of these two lineages highly expressed the PP marker gene (gsc), and the segregation between PP and Endo had already occurred by 6 hpf (hours post fertilization) (Figures 1G, 1H and S2B). Cell lineage prediction analysis on the Nodal explant datasets further indicated that the common progenitor cells were transcriptionally similar to PP (Figures 1I and S2I), which aligned with our previous observations in wild-type embryos (Figures 1D, 1E, S1H and S1I).

To explore the molecular signatures during the specification of PP and anterior Endo, we performed differential expression analysis for these two lineages in both wild-type embryos and Nodal explants. 18 genes were highly enriched along the trajectory of Endo specification in both wild-type embryos and Nodal explants (Figure S2C). These genes enriched GO terms such as “gland development”, “negative regulation of transcription, DNA-templated” and “negative regulation of RNA biosynthetic process” (Figure S2D). While for PP, 49 highly expressed genes were identified (Figure S2C), which enriched GO terms including “proteolysis involved in cellular protein catabolic process”, “cellular protein catabolic process”, “protein catabolic process”, etc. (Figure S2D). Cell fate specification requires the activation of a series of genes, and these gene sets can be regarded as gene modules, which are probably responsible for specific lineage formation²⁸. We identified two gene modules associated with Endo and PP lineages respectively (Figure S2E). Gene Ontology enrichment analyses²⁹ showed that terms like “endoderm development” and “cell motility” were enriched for Endo (gene module 1), while “embryo development ending in birth or egg hatching” was enriched for PP (gene module 2) (Figures S2F-S2H). These enriched GO terms were highly correlated with the developmental processes and functions of these two cell types.

To further validate this observation and investigate the dynamics of lineage segregation between anterior Endo and PP during development, we performed live imaging analysis on embryos of Tg(gsc:EGFP;sox17:DsRed) (Figures 2A and 2B). Notably, we observed a gradual emergence of sox17-positive cells from a subset of the gsc-positive cells (Figures 2A-2D, S3A-S3D and Movies S2, S3). And this cell fate transition was also evident in Nodal explants constructed from Tg(gsc:EGFP;sox17:DsRed) embryos (Figures S3E-S3H and Movie S1). Furthermore, we generated sox17:Cre and gsc:loxP-STOP-loxP-mCherry constructs (Figure S4D) and co-injected them with transposase mRNA into one-cell stage embryos (Figure S4A). Interestingly, distinct mCherry-positive cells were observed in the injected embryos, providing clear evidence that gsc-positive cells can give rise to sox17-positive cells (Figures S4B and S4C).

Figure 2.

We co-stained gsc and sox17 in embryos and in Nodal explants respectively at different time points during gastrulation (Figures 2E and S3I). By assessing the expression levels of sox17 and gsc in sox17-positive anterior mesendoderm cells, we observed a positive correlation between sox17 and gsc expression at the onset of gastrulation. However, with advancing development, sox17 and gsc expression exhibited a negative correaltion in both embryos and Nodal explants (Figures 2F and S3J). This suggests that the increase of sox17 expression during anterior Endo differentiation is accompanied by the simultaneous reduction of gsc expression. These findings align with our live imaging experiments (Figures 2A, S3B and S3F). In summary, combing single-cell tanscriptomic trajectory analyses and live imaging experiments using embryos and explants, we demonstrated that anterior Endo cells were derived from PP progenitor cells in zebrafish.

Nodal-Lefty regulatory loop is needed for PP and anterior Endo fate specification

In zebrafish, the PP and anterior Endo lineages are physically in close proximity and direct contact¹⁷. During development, interactions through cell-cell contacts between neighboring cell types often play a crucial role in promoting cell fate segregation^17,30. We hypothesized that communications between PP and anterior Endo cells may influence each other’s cell fate specification. We employed LIANA³¹, a ligand-receptor analysis framework, to investigate the molecular signals involved in these interactions using single-cell RNA-seq datasets. Our analysis revealed strong interactions between PP and Endo cells in both wild-type embryos and Nodal explants (Figures S5C-S5F). Interestingly, key molecular factors in the Nodal-Lefty regulatory loops (Figure 3D), such as ndr1, lft1 and lft2, were significantly enriched when setting PP as source cells and Endo as target cells (Figures 3A, S5A, S5B and S5G-S5I), while Nodal-Lefty interactions were reduced when setting the cell types conversely (Figures 3B and S5G-S5I). In addition to Nodal-Lefty signaling, Fgf and Wnt signalings were also significantly enriched in cell-cell communications between these two cell types in both embryos and Nodal explants (Figures 3A, 3B and S5G-S5I). This observation is consistent with previous findings that Fgf signaling interacts with Nodal to regulate mesoderm and endoderm separation^13,32.

Figure 3

The results presented above suggest that differences in the Nodal-Lefty signaling may exist between PP and anterior Endo cells. It is widely accepted that Nodal-Lefty functions as reaction-diffusion system to facilitate Nodal gradient formation and cell fate determination^33–37 (Figure 3D). As zebrafish possess redundant Nodal (Ndr1, Ndr2 and Ndr3) and Lefty (Lefty1 and Lefty2) paralogs, we firstly investigated the perturbation effects of individual Nodal or Lefty paralog on Nodal activity and its gradient. We generated an ectopic Nodal gradient in zebrafish animal pole cells by injection of ndr2 mRNA at 128-cell stage in wild-type embryos, ndr1 morphants and lft1 mutants (Figures 3E, S6A-S6F and S6H-S6M). By immunostaining and quantification of the pSmad2 signaling, we observed a significant enhancement of Nodal activity in lft1 mutants, while its level decreased upon ndr1 knockdown (Figures 3E, 3F and S6A-S6M).

To further determine how Nodal-Lefty is involved in the regulation of PP and anterior Endo separation, we performed single-cell transcriptomic sequencing for Nodal explants constructed from ndr1 morphants and lft1 mutants at shield stage (Figure 3G). We integrated these two single-cell RNA-seq datasets with the dataset of wild-type Nodal explants²². After applying unsupervised clustering, we identified a total of 9 different cell clusters (Figures 3H, S7A and S7B). Among these clusters, two mesodermal clusters named mesoderm_PP and mesoderm_NT (Figure 3H) were obtained. Mesoderm_PP cells transcriptionally resembled the progenitors of PP evidenced by the high expression of gsc and chrd, while mesoderm_NT cells were the progenitors of notochord characterized by the high expression of noto and tbxta (Figures S7A and S7B). We observed that a small cluster of sox32 and sox17 positive anterior Endo cells was close to the PP cells on the UMAP (Figure 3H), further supporting the transcriptional similarity between anterior Endo and the PP lineage. We then calculated the proportions of PP and Endo cells in the three datasets, and found that the proportions of both cell types increased in lft1-mutant explants, while decreased in ndr1-morphant explants compared to wild-type explants (Figures 3I and S7C).

To better estimate Nodal activity during PP and Endo cell fate segregation, we defined a Nodal score by regressing out the expression level of all the Nodal downstream targets, which were differentially expressed upon Nodal treatment²² and also possessed pSmad2 binding peaks near their genetic loci^22,38 (see Methods). A total of 29 Nodal downstream targets, whose expression exhibited a high correlation with Nodal activity, were used to calculate Nodal score (Figure 3J). We found Nodal score was dramatically decreased in ndr1 knockdown samples and slightly increased in lft1 mutant samples compared to that in wild-type samples (Figures 3K and S7D), which is consistent with the results of pSmad2 immunostaining (Figures 3E and 3F). We then examined the levels of Nodal score across all cell types in Nodal explants at 6 hpf, and found that axial mesendodermal cells, including PP, notochord and anterior Endo, displayed the highest Nodal scores (Figure S7E), which was consistent with previous knowledge that Nodal activity was the highest in the dorsal margin of zebrafish embryos at the onset of gastrulation^1,39,40. A further comparsion of Nodel scores between PP and Endo revealed that PP had a sightly higher Nodal score than Endo (Figure S7E). These analyses underscored the importance of Nodal signaling levels in PP and anterior Endo lineage specification.

A meticulous single-cell trajectory tree reveals the role of chromatin remodeling in PP and anterior Endo segregation

To thoroughly investigate the dynamic transcriptome changes during PP and anterior Endo separation, we reconstructed a single-cell trajectory tree using all PP and anterior Endo cells from three Nodal explant datasets at 6 hpf which was a crucial time point previously identified for the separation of these two lineages (Figure 4A). This trajectory tree captured a continuous cell-state-transition process from PP progenitors to anterior Endo, and a key branching point in this process was identified (Figures 4A and S7F). To elucidate the molecular dynamics in this branching process, we performed differential expression analysis for these two lineages around the branching point (Figure 4B). Three clusters of genes were identified, which were highly expressed in pre-branch (progenitors), left-branch (Endo) and right-branch (PP) respectively. GO enrichment analysis using these three gene sets identified terms related to mesendodermal cell fate specification (Figure 4B). For instance, “regulation of endodermal cell fate specification” was observed in Endo lineage, while “negative regulation of transcription by RNA polymerase II” was found in PP lineage. Interestingly, “chromatin organization” was significantly enriched in the pre-branch cells, and maintained high enrichment during PP specification, but was down-regulated upon endodermal lineage formation (Figure 4B). Furthermore, the GO term “chromatin organization” was also enriched in the differentially expressed genes of the PP cells when comparing lft1 mutant or ndr1 knockdown samples to wild-type samples (Figures 4C-4G and S7G).

Figure 4.

These analyses unveiled a potential role of chromatin regulators in the regulation of PP and anterior Endo segregation, prompting us to identify the specific factors involved in this process. Among the candidates within GO term of “chromatin organization” (Figure 4B-4E), we identified a key member of SWI/SNF complex, srcap, which was differentially expressed between PP and Endo cells in both embryos and Nodal explants (Figures 4J, 4K and S7H). Previous studies have demonstrated that the SWI/SNF nucleosome remodeling complex is necessary for TGF-β-induced transcription of numerous target genes^41,42. To investigate the role of the SWI/SNF complex in the separation of PP and Endo cell fates, we treated embryos with AU15330 (Figure 4H), a known degrader of SWI/SNF ATPase components⁴³. Subsequently, we assessed the cell fate specification of PP and anterior Endo lineages by HCR co-staining for gsc, frzb and sox32. We found that the numbers of anterior Endo cells were significantly increased in the embryos treated with AU15330 compared to DMSO-treated embryos (Figures 4H and 4I). Similar effects were observed in srcap morphants, where the cell number of anterior Endo increased compared to wild-type embryos (Figures 4L, 4M, and S7I).

Collectively, these findings underscore the importance of SWI/SNF complexes in ensuring accurate cell fate segregation within the anterior mesendoderm.

Integrative analysis of single-nucleus RNA-Seq and ATAC-Seq reveals distinct chromatin states underlying differential expression profiling in PP and Endo cells

To delve deeper into the involvement of chromatin states in the separation of PP and Endo, we conducted single-cell assays for transposase-accessible chromatin (ATAC) and nuclear RNA sequencing in zebrafish embryos at 6 hpf (Figure 5A). Recognizing that accurate cell clustering cannot solely rely on chromatin accessibility dataset, we integrated this dataset with RNA expression data to achieve precise cell type identification. The resulting cells were clustered into 10 different groups on the UMAP (Figures 5B, S8A, and S8B). Consistent with the earlier single-cell RNA atlas of Nodal explants, Endo was found to locate closely to the axis mesoderm (PP and notochord progenitors) on the joint UMAP (Figure 5B).

Figure 5.

Using the marker gene gsc, we identified PP progenitor cells within the axis mesoderm, and then Endo and PP progenitor cells were selected for downstream analysis (Figure 5C). A small cluster of cells expressing sox32, sox17 and foxj1a, which were located distantly from PP and Endo, were classified as dorsal forerunner cells (DFCs) and excluded from further analysis (Figures 5C and S8C).

Our data unveiled differential chromatin accessibility between PP and Endo cells (Figures 5G, S8K, and S8L). Correlation analysis revealed a strong correspondence between chromatin accessibility and RNA expression levels for marker genes in these two cell types (Figure S8D). To identify master regulators of PP and Endo separation, we performed differential analysis for both chromatin accessibility and gene expression (Figures 5D and S8E). A total of 170 up-regulated genes in PP, including gsc and ripply1, and 261 up-regulated genes in Endo, such as sox32 and sox17, were identified. Most of these genes also exhibited differential chromatin opening states in PP and anterior Endo lineages (Figures 5E-5G and S8F-S8I).

Previous study reported that Nodal signaling can boost chromatin accessibility both in vivo and ex vivo^44,45. And we observed that PP displayed a higher Nodal levels compared to Endo cells (Figures S2J and S7E). In this ATAC and RNA-seq dataset, we also observed a higher Nodal score in PP cells (Figures 5H and 5I). Thus, it is possible that the differences in Nodal activity may lead to differential chromatin accessibility states and gene expression level in PP and Endo.

To further determine how Nodal siganling levels may differentially influence the gene expression dynamics, we first explored the expression profiles of key markers of anterior mesendoderm in the bulk RNA-seq data of Nodal explants injected with different dosages of Nodal mRNA²². We observed that the highest Nodal dosages did not result in further elevation of the expression of the endodermal marker sox17 (Figure 5J), while higher concentrations of Nodal promoted the expression of gsc and ripply1 (Figure 5J). Importantly, these patterns of differential marker expression in response to Nodal/Activin concentration were not only observed in our Nodal explant model but also validated in a human embryonic stem cell system (Figure 5K). These findings suggest that Nodal signaling levels play a conserved role in determining the cell fate specification, which influences the differentiation towards either Endo or PP lineages. To further validate this finding, we injected different dosages of Nodal mRNA into one zebrafish animal pole blastomere at the 128-cell stage, and performed HCR co-staining for sox32 and frzb at 6 hpf (Figure S9A). By quantifying the cell numbers of Endo and the area sizes of PP cluster, we found that moderate concentration of Nodal levels promoted Endo cell specification (Figures S9B-S9E).

Based on our aforementioned findings, it becomes evident that subtle variation in Nodal signaling levels may regulate PP and Endo gene expression. This variation is likely regulated by Nodal-Lefty regulatory loops, and contributes to the establishment of distinct chromatin states and transcriptional response of key transcriptional regulators within these two lineages.

Gsc and Ripply1 collaborate to suppress the specification of anterior Endo in zebrafish

It has been reported that endoderm differentiation could be suppressed by upregulation of gsc¹⁸, while the role of ripply1 in endoderm development remains unknown, though it has been characterized as a transcriptional repressor ^46,47. Moreover, both Gsc and Ripply1 are capable of inducing a secondary axis formation when overexpressed in the ventral side of zebrafish embryos^19,22. Combining these findings from literature with our observation of high expression of ripply1 in PP of Nodal explants and zebrafish embryos at 6 hpf (Figures 5D and S8I), we hypothesized that Ripply1 may also play a role in suppressing Endo specification. To test this hypothesis, we generated mutant lines for gsc and ripply1 in zebrafish (Figure 6B). In order to comprehensively analyze the regulatory role of Ripply1 in cooperation with Gsc in the specification of cell fates between PP and Endo cells, we conducted a self-cross of the Tg(gsc^+/-;ripply1^+/-) lines (Figure 6A). Subsequently, we collected the descendant embryos to perform HCR co-staining for frzb and sox32, along with genotyping (Figure 6A). From the descendant embryos, we obtained a total of 9 genotypes with different phenotypes displaying variable defects in PP and anterior Endo separation (Figures 6C and S10A). Interestingly, we observed an increase in the number of anterior Endo cells with the loss of wild-type alleles of gsc and ripply1 (Figures 6D, S10B, and S10C). To further investigate the mechanisms of Ripply1 suppressing anterior Endo specification in zebrafish, we injected ripply1 mRNA at 1-cell stage in zebrafish embryos and found that Endo cell number was dramatically decreased as shown by whole-mount in situ hybridization (WISH) of sox32 (Figure 6F). Additionally, we constructed a plasmid that contained ripply1 cDNA downstream of the sox17 promoter. Injection of this plasmid at 1-cell stage also resulted in a decrease in Endo specification shown by EGFP expression (Figure 6E). Taken together, these results demonstrated that Ripply1 acted as an Endo repressor in zebrafish.

Figure 6.

It is well established that Gsc functions as an Endo suppressor by directly binding to the promoter region of sox17¹⁸. However, the mechanism by which Ripply1 suppresses Endo specification remains unknown. To address this question, we injected HA-tagged ripply1 mRNA into zebrafish embryos for conducting CUT&Tag experiments (Figures 6G, S11A and S11B). The CUT&Tag sequencing dataset was analysed following the standard process (Figures S11C-S11H), and 72,683 peaks were enriched in HA-ripply1 injected samples (Figure 6H). Genes nearby the differentially enriched peaks were used to perform GO enrichment analysis (Figure 6I). Notably, we observed significant enrichment of terms related to “muscle structure development”, “skeletal system development” and “connective tissue development” (Figure 6I), which were consistent with the well-established functions of Ripply1 reported in the published studies^46–48. Additionally, we also observed significant enrichment of GO terms such as “regulation of cell differentiation”, “negative regulation of signaling” and “negative regulation of transcription by RNA polymerase II” (Figure 6I), further supporting the transcriptional repressive function of Ripply1. More interestingly, Ripply1 binding peaks were found to be enriched in the immediate upstream region of the gene body of sox17 and sox32 (Figures 6J and 6K), suggesting a direct function of Ripply1 in the transcriptional regulation of these two endodermal markers.

In summary, our data proposes a model within PP and anterior Endo lineage separation. Variations of Nodal activities promote the establishment of distinct chromatin states in PP and anterior Endo lineages, which is facilitated by SWI/SNF complexes. These differential chromatin accessibility patterns, in turn, drive the differential expression of key regulators such as gsc and ripply1 in PP and Endo lineages. Collectively, these regulators collaborate to finely regulate the segregation of mesendoderm cell fates (Figure 7).

Figure 7.

Discussion

In this study, we delved into the intricacies of specifying and segregating PP and anterior Endo lineages in zebrafish through single-cell multi-omics and live imaging analyses. Our findings revealed the origin of anterior Endo from PP progenitors and unveiled a correlation between subtle differences of Nodal signaling levels and lineage commitment. These signaling variations were identified as potential regulators influencing the expression of cell-fate-specific genes at both epigenetic and transcriptional levels, with the collaboration of SWI/SNF complexes. Through gain-/loss-of-function studies, we demonstrated the collaborative suppression of Gsc and Ripply1 in anterior Endo specification.

While the induction of mesendodermal cell fates has been extensively studied^9,12,49,50, the mechanisms underlying the separation of meso-and endo-lineages from common progenitors remain a subject of ongoing research. The widely accepted notion posits that high Nodal signaling levels promote endoderm specification, while lower levels induce mesoderm^11–13. These quantitative effects of Nodal signaling lead to differential activation of transcriptional targets which have differential sensitivity to Nodal activation^2,38. However, recent studies challenge this perspective, proposing that Nodal signal level may not strictly determine cell fate separation of mesendoderm. Instead, Nodal signaling provides competency for mesendodermal progenitors to stochastically switch and commit to lineage-specific fates in zebrafish¹⁶. Moreover, the duration of Nodal signaling has been implicated in specifying different cell fates of mesendoderm, with prolonged signaling promoting prechordal plate specification and suppressing endoderm induction¹⁸. Besides, other important signaling pathways in early embryonic development also interact with Nodal signaling to regulate cell fate separation^13,32. It has been well established that Fgf signaling plays a crucial role in regulating cell fate switching between endoderm and mesoderm in lateral margin of zebrafish embryo. In our current work, we preliminarily explored whether Fgf siganling plays roles in regulating cell fate segregation between PP and anterior Endo cells (Figure S12). We observed that several genes related to Fgf siganling were highly enriched in PP and Endo cells (Figure S12B), and inhibition of Fgf signaling could obviously increase the cell number of anterior Endo (Figures S12A and S12C). The findings above suggest that Fgf signaling inhibits the specification of anterior Endo from PP progenitors in zebrafish, which is a similar process to the suppression observed during lateral endoderm specification. However, further investigation is required to understand the specific interactions between Fgf signaling⁵¹ and factors such as Gsc, Ripply1, Sox32^52,53 and the SWI/SNF complexes, as well as their roles in regulating the fate diversification of PP and Endo lineages.

Our study contributed a new layer of understanding to the regulation of mesendodermal cell fate segregation by Nodal signaling. We proposed that Nodal signaling may influence chromatin remodeling factors, thereby affecting the chromatin accessibility states of key regulators such as gsc and ripply1. It is worth noting that gsc and ripply1 have been identified as direct targets of Nodal signaling²². This suggests that the transcription of these genes may directly and sensitively respond to variations in Nodal signaling levels. Based on this understanding, we hypothesize that Nodal signaling can achieve more precise and robust control over the fate segregation of closely related cell types by introducing an additional layer of epigenetic regulatory mechanisms. Additionally, our findings underscored the pivotal role of chromatin openness in morphogen interpretation during development, emphasizing the need to consider chromatin structure when studying how morphogen gradients regulate patterning formation across different species.

In our work, combining single-cell multi-omics analyses with perturbation experiments, we revealed that chromatin remodeling together with Nodal morphogen promoted the PP and Endo separation in zebrafish. However, which and how chromatin remodeling factors respond to Nodal siganling to facilitate fate diversification remains largely unknown. Our current study identified a key member of SWI/SNF complexes, srcap, which is differentially expressed in PP and Endo cells. What triggers the transcriptional difference of srcap in these two lineages is still a open question. One of the hypotheses is that the transcription of srcap is very sensitive to Nodal signaling (Nodal concentration or Nodal duration), and a little variance in Nodal signaling can drive differential activation of srcap. Another speculation is that the activation of srcap needs the assistance of pioneer factors, like forkhead family^41,44, whose binding motifs are particularly accessible in PP cells compared to Endo cells (Figures S8K and S8L). All these hypotheses need to be further determined in future works.

Despite these insights, several questions remain unanswered. Firstly, we discovered a potential role of Nodal activity levels in the separation of PP and Endo cells. However, further investigation is required to understand how Nodal activity levels and signaling duration are integrated to regulate this process. Secondly, although our study determined the cooperative role of Gsc and Ripply1 in suppressing endoderm specification, the potential involvement of other unidentified transcriptional repressors, as suggested in previous work¹⁸, warrants further investigation. Notably, even in double homozygous mutants of gsc and ripply1, the specification of PP clusters was still detected, and the enrichment of another transcriptional repressor, osr1, in PP cells (Figure S8J) suggested a complex network of cooperative endoderm repressors in preventing cells from differentiating into Endo. Thus, it would be worth studying how Gsc, Ripply1, Osr1 and other potential factors cooperatively regulate PP and Endo specification in the future. Lastly, the anterior endoderm and ventral-lateral endoderm were mixed together in our single-cell ATAC-seq dataset, which brought in interferences when we interpreted the mechanisms of the separation between PP and anterior Endo. More replicates of embryonic single-cell multi-omics or performing single-cell multi-omics on Nodal explants may help to address this issue in future studies.

Methods

Animal ethics

Zebrafish strains were conducted following standard procedures, and experimental procedures were approved by the Institutional Review Board of Zhejiang University. The protocol number is ZJU20220375.

Generation of Nodal explants

Nodal explant was generated as our previous work described⁵⁴. 10 pg ndr2 mRNA was injected into one blastomere of zebrafish embryonic animal pole at the 128-cell stage. About half of the animal pole region of the embryo was cut off from the blastula at the 1k-cell stage, and was cultured in a Petri dish coated with 1.5% agarose gel filled with Dulbecco’s modified Eagle’s medium/nutrient mixture F-12 (DMEM/F-12) with 10 mM HEPES, 1x minimum essential medium (MEM) containing non-essential amino acids and supplemented with 7 mM CaCl₂, 1 mM sodium pyruvate, 50 μg/ml gentamycin, 100 μM 2-mercaptoethanol, 1× antibiotic-antimycotic (15240062, Thermo Fisher) and 10% serum replacement. For generating ndr1-morphant Nodal explants, ndr1 morpholino (0.5 mM, 2 nL) was injected into yolk of wild-type zebrafish embryos at 1-cell stage; while for generating lft1-mutant Nodal explants, the embryos of lft1 mutants were used to construct Nodal explants. The sequence of ndr1 morpholino was 5’ – ATGTCAAATCAAGGTAATAATCCAC - 3’.

Live imaging

Tg(gsc:EGFP;sox17:dsRed) embryos and Nodal explants generated from those embryos were used to perform live imaging analysis. 0.5% low-melt agarose was used to mount the embryos or Nodal explants in glass-bottom dishes (Cellvis, D35-20-0-N). The embryos or explants were live-imaged by confocal laser scanning microscopy (OLYMPUS FV12-IXCOV) using a 20X objective lens or a 40X oil immersion lens.

In situ hybridization chain reaction (HCR)

HCR co-staining of frzb, sox32 and gsc was performed as previous studies described²¹. The probes of those genes were ordered and produced from Molecular Technologies (http://www.moleculartechnologies.org/). Embryos were fixed in DEPC-treated PBS with 4%(w/v) paraformaldehyde at 4°C overnight. The embryos were incubated in 500 μL 30% formamide probe hybridization with 2 pmol of each HCR probe set at 37°C overnight. Next day, 30% formamide probe wash buffer was used to stop hybridization by repeated washing of the embryos at 37°C. Fluorescent signals were generated and amplified by probes that bound to fluorescent HCR amplifiers in an amplification buffer overnight at room temperature. To halt this process, the samples were washed several times using 5× SSC with 0.001% Tween 20. The reaction buffers, fluorescent hairpins and probes were manufactured by Molecular Technologies. OLYMPUS FV12-IXCOV confocal microscope was used to photograph those HCR co-staining samples.

Double-color RNA-fluorescence in situ hybridization (FISH)

Double-color FISH of gsc and sox17 was conducted according to the protocol from manufactory (http://www.pinpoease.com/, GD Pinpoease Biotech Co., Ltd.). Wild-type embryos or Nodal explants (6 hpf, 7 hpf and 8 hpf) were fixed in DEPC-treated PBS with 4%(w/v) paraformaldehyde at 4°C overnight. Pre-A solution was used to inhibit the peroxidase activity. Target RNA molecules were exposed through protease treatment, followed by hybridization with probes at 4°C for 2 hours. The fluorescent signal was amplified through sequential reactions 1, 2, and 3. Lastly, a tyramide fluorescent substrate (OpalTM520, Akoya Biosciences) was used to incubate the embryos, enabling fluorescent labeling of the target RNA through the Tyramide Signal Amplification (TSA) assay.

Immunostaining and quantification of pSmad2

Samples were fixed at 4% paraformaldehyde overnight for pSmad2/3 immunostaining. The experiment was performed as previously described⁵⁴. Confocal laser scanning microscopy (OLYMPUS FV12-IXCOV) was used to scan the immunostained samples. Imaris (Version 9.7) software was used to quantify the signaling activity of pSmad2 (488 nm), RFP (561 nm) and DAPI (405 nm). Firstly, the raw image data was loaded in Imaris and visualized in a three-dimensional (3D) view. The Spot plugin from Imaris was used to detect spots with pSmad2 signal. The diameter of spot detection was set to 5-6 μm for GFP channel. The spots were automatically identified with default parameters. And then, as the identification of some spots was inaccurate, we manually added or removed some spots. Lastly, the coordinates and signal intensities of each spot were exported for next-step analyses. To make all images comparable, we used DAPI signal as an internal control to normalize the signal intensity in other channels. When plotting a pSmad2 signal intensity gradient along a distance, the cell with the highest pSmad2 signal was assigned as the start point, and then the distance of other cells was calculated by comparing them with the start point cell. The fitting curve was plotted by ggplot2⁵⁵ in R (https://www.r-project.org).

Whole-mount in situ hybridization (WISH)

The processes of WISH and the construction of sox32 probes were described in our previous work⁵⁴. Embryos were fixed in DEPC-treated PBS with 4%(w/v) paraformaldehyde at 4°C overnight. To enable long-term storage, the embryos were dehydrated using 100% methanol. For WISH experiments, probes were labeled with digoxigenin-11-UTP (Roche Diagnostics), and the substrate was NBT/BCIP.

Activation of Activin/Nodal signaling in human embryonic stem cell system

Human embryonic stem cell (hESC) lines were used between passage 45-50. Cells were seeded equally on matrigel-coated plate and cultured in mTeSR medium. Before the experiment, each group was pretreated in N2B27 medium for 6 hours. Then the control group was maintained in N2B27, while the experiment groups were cultured in N2B27 medium supplemented with different concentrations of Activn-A (0, 10, 40, 70 and 100 ng/mL). The cellular RNA of each group was extracted with Trizol after 12 hours of treatment.

Overexpression of gsc and ripply1 in zebrafish embryos

The mRNA of gsc and ripply1 (synthesized in vitro) and sox17 promotor-driven ripply1 plasmid were co-injected into zebrafish embryos at 1-cell stage. The ddH₂O was injected as control. The embryos were imaged or fixed at 8 hpf. The sox32 probe was used for WISH.

Generating heterozygous embryos of gsc and ripply1

The ripply1 mutants were obtained from Dr. Ming Shao’s lab at Shandong University. We knockout gsc in ripply1-mutant embryos by CRISPR-Cas9 system. F0 embryos were raised up, and then crossed with wild-type (AB). The descendant embryos (F1) were raised up and genotyped. The adult zebrafish of the target genotype were selected.

Quantifying the cell number of anterior Endo

Imaris (version 9.7) was used to identify the Endo cells based on the HCR images. Briefly, the raw image data was imported into Imaris, and the object was visualized in a 3D view. The Spot plugin from Imaris was used to detect each cell. The diameter of spot detection was set to 10-12 μm for GFP channel (staining for sox32 or sox17). The candidate spots were selected by setting an appropriate quality threshold. Finally, the coordinates and signal intensities of each spot were exported from Imaris and used for further analysis.

CUT&Tag sample preparation

HA-ripply1 sequence with homologous sequences was generated through two PCR steps using two separate sets of primers (Table S1, reverse primer was common, ripply1-F1 primer: tcaaggcctctcgagcctctagaATGTACCCATACGATGTTCCAGATTACGCT; ripply1-F2 primer: ACGATGTTCCAGATTACGCTGGCAGCATGAATTCTGTGTGCTTTGCCACT; ripply1-R primer: TACGACTCACTATAGTTCTAGAtcagttgaaagctgtgaagtga). The details of generating mRNA in vitro was described in our previous protocol²¹. Briefly, the HA-ripply1 PCR products were inserted into pCS²⁺ plasmid through homologous recombination using ClonExpress II One Step Cloning Kit. The plasmid containing HA-ripply1 was linearized using Not1 cleavage, and HA-ripply1 mRNA were generated through mRNA in vitro synthesis using mMESSAGE mMACHINETM SP6 transcription kit. Wild-type embryos were injected with HA-ripply1 mRNA or ddH₂O at one-cell stage and harvested at 6 hpf. Yolk was removed through pipetting, and the cells were collected. The cell viability was assessed by Trypan blue staining and the cell number was then counted.

CUT&Tag experiment

CUT&Tag was performed on samples of wildtype, HA-ripply1-injected and HA-ripply1-injected without primary antibody using Hyperactive Universal CUT&Tag Assay Kit for Illumina Pro (Vazyme, TD904), and each sample had two replicates. Briefly, around 100,000 cells were collected and immobilized on Concanavalin A-coat Magnetic Beads. Cells were then incubated with primary antibody (HA-Tag Rabbit mAb, Cell Signaling Technology, 3724S, 1:50 dilution) at 4°C for 12 h-16 h. Tethered cells were washed thoroughly to remove excessive primary antibody and incubated with secondary antibody (Goat Anti-Rabbit IgG H&L, Vazyme, Ab207-01, 1:100 dilution) for 60 minutes at room temperature. Then samples were washed thoroughly to remove unbound secondary antibody and incubated with pA/G-Tnp for 1 hour at room temperature. DNA was then fragmentated through Trueprep Tagment Buffer L (TTBL). Fragmented DNA was extracted using DNA extract beads and amplified through PCR. PCR products were purified using VAHTS DNA CLEAN Beads (Vazyme, N411) and then sent to Novogene Co., Ltd. (Beijing, China) for sequencing. Paired-end sequencing at 150 bp read length was performed on an NovaSeq 6000 System.

CUT&Tag analysis

Sequencing data of CUT&Tag samples were first treated by Cutadapt v4.5⁵⁶ to remove adapters and trim reads using the following parameters “-m 18 -q 30,30 --max-n=0.05 -e 0.2 -n 2”. The trimmed reads were then aligned to zebrafish genome (danRer11) through Bowtie2 v2.5.2⁵⁷ using the following parameters “--end-to-end--very-sensitive--no-mixed--no-discordant -- phred33 -I 10 -X 700”. Samtools v1.18⁵⁸ and BEDtools v2.31.0⁵⁹ were then used for post-alignment processing to covert file formats in preparation for the peak calling and visualization. BigWig tracks were generated for visualization of chromatin landscapes in interested regions using the Integrative Genomics Viewer (IGV)⁶⁰. To assess the reproducibility between replicates, the genome was spilt into 500 bp bins, and a Pearson correlation of the log₂-transformed read counts in each bin was calculated between datasets.

SEACR⁶¹ was used to call “non stringent” peaks for wild-type and HA-ripply1-injected samples using merged HA-ripply1-injected samples without primary antibody as control. Additionally, top 1% of enrich regions were selected as peaks by AUC without controls. Peaks were annotated by ChIPpeakAnno v3.34.1⁶² using UCSC annotation from TxDb.Drerio.UCSC.danRer11.refGene v3.4.6⁶³. Significantly enriched peaks in HA-ripply1-injected samples were identified using DESeq2 v1.40.2⁶⁴, and enriched GO terms were identified using clusterProfiler v4.9.0.2⁶⁵. Signal intensity heatmaps were generated using deeptools v3.5.4⁶⁶. The peaks of HA-ripply1-injected rep1 were used as the reference to generate computeMatrix using the reference-point model (--skipZeros -- afterRegionStartLength 3000--beforeRegionStartLength 3000--referencePoint center), and PlotHeatmap function was used to plot signal intensity.

Processing and analyzing single-cell RNA-seq data

Illumina sequencing reads were aligned to the reference genome (GRCz11) of zebrafish through 10× Genomics CellRanger pipeline (version 3.0.2) with default parameters. Nodal explants generated from lft1 mutants and ndr1 morphants yielded 5,419 and 5,195 cells respectively. The expression matrix of each sample was obtained via CellRanger pipeline analysis and used for further analysis. The single-cell RNA-seq data of wild-type Nodal explants was reused from our previous study⁵⁴.

Processing and analyzing single-cell multi-omics data

Illumina sequencing reads were aligned to the reference genome (GRCz11) of zebrafish using 10× Genomics CellRanger-arc pipeline (version 2.0.0) with default parameters. Expression matrix and fragments files were obtained after running the pipeline. 4,879 cells were obtained. Signac⁶⁷ was used for downstream analyses with default parameters. Low-quality cells were filtered out with nCount_ATAC < 100000 & nCount_RNA < 50000 & nCount_ATAC > 1000 & nCount_RNA > 1000 & nucleosome_signal < 2 & TSS.enrichment > 1 & percent.mt < 4. And then, gene expression data and DNA accessibility data were processed respectively by Signac with default parameters. The RNA expression level was utilized to identify cluster markers for cell annotation. Gene activity matrix was obtained by assessing the chromatin accessibility associated with each gene through GeneActivity function.

Differential expression analysis of PP cells between wild-type, lft1-mutant and ndr1-morphant Nodal explants

The three sing-cell RNA-seq datasets of Nodal explants constructed from wild-type embryos, lft1 mutants and ndr1 morphants were integrated by Seurat⁶⁸. To investigate the functional and regulatory effects of Nodal signals on cell differentiation, we compared the wild-type mesoderm PP cells to the mesoderm PP cells in lft1 mutants and ndr1 morphants respectively, and differentially-expressed (DE) genes were identified through Seurat v4.0.2 (https://github.com/satijalab/seurat) with a Bonferroni adjusted p-value < 0.05. As it is known that EVL is not affected by Nodal signaling^54,69, DE genes identified in EVL were regarded as background noises and were removed from the mesoderm PP DE genes. The mesoderm PP DE genes positively correlated with Nodal concentration were then selected, constructing two gene sets: upregulated DE genes in lft1 mutants and downregulated DE genes in ndr1 morphants. GO enrichment analysis was performed on these two gene sets respectively using clusterProfilter v4.2.2 (https://github.com/YuLab-SMU/clusterProfiler). The enriched GO terms in “biological process” subontology were identified using a Benjamini-Hochberg adjusted p-value < 0.05.

Cell-cell communication analysis by LIANA

R package, LIANA³¹ v0.1.5 (https://github.com/saezlab/liana/), was used for the investigation of ligand-receptor interactions among different cell types. The ligand-receptor analysis was performed against the ligand-receptor reference achieved from Omnipath resource using 5 different methods implemented in LIANA, including Natmi, Connectome, SingleCellSingalR, iTALK and CellPhoneDBv2. The consensus rank aggregated from the results of these methods was used for the identification of enriched ligand-receptor interactions.

Inferring cell-cell communication by CellChat

CellChat⁷⁰ was employed to systematically analyze cell-cell communication based on prior known zebrafish ligand-receptor interaction database CellChatDB. As some known interactions were missing in the database, we manually added some interactions of interest into the database. The expression data of scRNA-seq was preprocessed to identify over-expressed ligands and receptors for each group. CellChat inferred communication probability between two interacting cell groups based on the average gene expression of a ligand in one cell group and the average gene expression of a receptor in another cell group. The communication probabilities of signaling pathways were then calculated by summarizing the communication probabilities of all ligands-receptors interactions associated with each pathway.

Gene set enrichment analysis

The wild-type mesoderm PP gene expression was compared to the mesoderm PP gene expression of lft1 mutants and ndr1 morphants respectively using Seurat v4.0.2 (https://github.com/satijalab/seurat). The genes expressed in more than 10% cells in either of the two compared populations were retained. A ranked list was formed on the retained genes using sign(log₂FC) * (-log₁₀PValue) as the ranking statistic. The mapping between GO terms and zebrafish genes was achieved through org.Dr.eg.db v3.14.0 (https://bioconductor.org/packages/org.Dr.eg.db/). GSEA was performed on the selected GO terms by clusterProfiler v4.2.2 (https://github.com/YuLab-SMU/clusterProfiler), using fgsea (https://github.com/ctlab/fgsea) with 10⁵ interactions. R package, enrichplot v1.14.2 (https://github.com/YuLab-SMU/enrichplot), was then used to visualize the GSEA results.

Single-cell RNA-seq trajectory analysis

Three different methods were used to perform single-cell RNA-seq trajectory analysis. The developmental trajectory tree of zebrafish embryos and Nodal explants was constructed by URD⁷¹ with default parameters as described in previous studies⁷¹. The branchpoint preference plot of PP and Endo was constructed by URD with default parameters. Firstly, branchpointPreferenceLayout function was used to define the preference layout for the branchpoint. And then stage, pseudotime information and gene expression were plotted on the branchpoint using plotBranchpoint function. The expression matrix of PP and Endo cells was also used to construct trajectory tree by monocle2²⁸ and monocle3⁷² with default parameters. And differential expression analysis was performed at the branching point of PP and Endo separation on trajectory tree constructed by BEAM function in monocle2. Modules of co-regulated genes along PP or Endo differentiation lineage were identified by find_gene_modules function in monocle3 with default parameters. The cells contributed to the preference layout for the Endo and PP brachpoint from URD analysis were also selected as inputs for Palantir²⁷. We used Harmony⁷³ to calculate the augmented affinity matrix for data of all time points. Data was visualized using force directed layouts. Palantir was used for downstream analysis, and a nanog^high cell was used as the start cell. A frzb^high cell (PP) and a sox17^high cell (Endo) were selected as the terminal cells. Endo and PP cell differential trajectory was detected by palantir.core.run_palantir function from Palantir. And finally, gene expression trend along Palantir inferred pseudotime was calculated by palantir.presults.compute_gene_trends function.

Calculating Nodal score in single-cell RNA-seq data

Nodal functions through activating a set of key downstream genes in a concentration-dependent manner. Our previous study identified 105 Nodal immediate target genes by analyzing bulk RNA-seq datasets of Nodal mRNA injected explants²². We overlapped these 105 genes with 61 Nodal direct targets, which were previously identified through ChIP-seq analysis of pSmad2³⁸. This analysis allowed us to pinpoint 29 Nodal downstream genes whose expression was highly sensitive to Nodal activity. To quantify Nodal activity, we defined a Nodal score based on its downstream transcriptional response. This score was calculated by regressing out the expression levels of these 29 Nodal downstream genes. To accomplish this, the AddModuleScore function in Seurat was employed.

Data availability

The data generated in this study is deposited in the Gene Expression Omnibus (GEO). The single-cell RNA sequencing and single-cell multi-omics data is under accession number GSE223636, and the CUT&Tag data is under the accession number GSE249292.

Code availability

Any custom code and data are available from the authors upon request. All analyses are based on previously published code and software.

Acknowledgements

We would like to thank Dr. Alex Schier for his kind gift of the lft1 mutant. We thank Dr. Ming Shao from Shandong University for his kind gift of the ripply1 mutant. We thank Dr. Peng Xia from Zhejiang University for his kind gift of the Tg(gsc:EGFP) embryos. We also thank Dr. Nai-He Jing at University of Chinese Academy of Sciences, Dr. Jun Ma, Dr. Xiao-hang Yang, Dr. Min-Xin Guan, Dr. Feng He and the members of Laboratory of Development and Organogenesis (LDO) at Zhejiang University for helpful suggestions and discussions. We thank Mr. Guang-Xu Zhang from Biostar Technology for the technical help on constructing 10x single-cell multi-omics sequencing library. We thank Shuang-Shuang Liu from the Imaging Platform and Ying-Niang Li from the zebrafish core facility at Zhejiang University School of Medicine for their technical support. This work was supported by grants from the Chinese National Key Research and Development Project (2022YFA1103100, 2019YFA0802402) and the National Scientific Foundation of China (32050109, 32300688, 32300677).

Conflict of Interest

The authors declare no competing interests.

Author Contributions

PFX, TC and HQL conceived and designed the research; TC, XL, YD, YMT, YYX, CYC, CL, YFL, YH, DHZ, XX, XXF, ZXJ, JYL and JM performed research; TC and YD analyzed the single-cell RNA-seq and single-cell multi-omics data; YD performed the CUT&Tag experiment and analyzed the data; TC, YD, HQL and PFX wrote the manuscript; all authors reviewed and approved the manuscript. PFX and HQL supervised this study. Fundings for the study were provided by PFX, TC and YD.