Integrated multi-omic atlas reveals the hierarchy of spatiotemporal regulatory networks of mouse gastrulation

Domain: biology
Published: 2026-01-12
Editorial outcome: major_revision

Abstract
Spatiotemporal coordination of cellular and molecular events is crucial for cell fate commitment during mouse gastrulation. However, the high-precision mechanisms governing the timing and spatial dynamics remain poorly understood. Here, we present a time-series single-cell multi-omic dataset of the gastrulating mouse embryos and construct a hierarchical gene regulatory landscape. Integrating this with real three-dimensional transcriptomic coordinate, we created ST-MAGIC and ST-MAGIC(+) atlas, dissecting the spatiotemporal logics of regulatory networks and signaling responsiveness underpinning the lineage commitment at gastrulation. Specifically, we delineated the multi-omic basis for left-right symmetry breaking events in the gastrula and also revealed the spatiotemporal molecular relay for axial mesendoderm lineage, where early and intermediate transcription factors first open the chromatin regions and setup the responsiveness to signaling, followed by terminal factors to consolidate the transcriptomic architecture. In summary, our study presents a spatiotemporal regulatory logic framework of mouse gastrulation for advancing our understanding of mammalian embryogenesis.

Introduction
Gastrulation is a pivotal phase in embryonic development, at which the multipotent epiblast is allocated to diverse tissue lineages of the three primary germ layers: ectoderm, mesoderm, and endoderm1,2. This stage is also crucial for the spatial organization of the embryo in anterior-posterior, dorsal-ventral and left-right axes1. Gastrulation, therefore, lays the blueprint of embryogenesis, making it a focal point for understanding the molecular and cellular mechanisms that govern embryonic patterning.

The dynamic nature of gastrulation involves rapid changes in cell composition, differentiation and proliferation, all of which are intricately regulated by a combination of epigenetic modifications, transcriptional networks, and intercellular communication. Recent advancements using limited cell numbers or single-cell omics significantly enhance our understanding of these dynamic processes during gastrulation. These technologies have enabled high-resolution analysis of gene expression patterns and cellular heterogeneity within developing embryos, identify that embryogenesis during gastrulation is orchestrated by precise schedules of gene expression within individual cells, regulated by complex intrinsic gene regulatory networks (GRNs) and influenced by intercellular signaling pathways3–6. These findings highlight the importance of both temporal and spatial regulation in cell fate and tissue patterning.

Despite these advances, a comprehensive characterization of the GRNs that govern the rapid cell fate commitment and spatial patterning during gastrulation remains elusive. Traditional single modality of single-cell omics often disrupts the spatial and temporal continuity of embryonic development, limiting our ability to fully understand the molecular mechanisms underlying embryo patterning. Spatiotemporal coordination of GRNs provides the developmental principles of lineage allocation and tissue patterning in the germ layers of the gastrulating embryo.

One illustrative example is the formation of the axial mesendoderm lineage in mice, a distinct structure emerging from the dorsal organizer during the late-streak stage. The anterior mesendoderm and primitive node cells form the prospective notochordal structure, with the primitive node acting as a reservoir for continuous production of caudal notochord cells7–9. The hereby forming notochord structure serves as an organizing center for proper dorsal-ventral and left-right body axes10–12. Disruptions in the notochord structure often result in severe embryonic defects, emphasizing its critical role in body axis formation13–17. However, the sequential cellular events and molecular GRNs with spatial-temporal information underlying axial mesendoderm formation are not well characterized.

To address these challenges, we generated a high-resolution, spatiotemporally resolved single-cell multi-omic reference for the mouse gastrulating embryo. Our dataset includes transcriptomics and chromatin accessibility profiles, collected at 6 h intervals across multiple stages of gastrulation. Using the developed Bi-Orientation Cis-Regulatory Elements Predictor (BioCRE) algorithm, we linked genes with their potential regulatory elements, establishing a comprehensive GRNs atlas underlying the rapid cell lineage commitment process. With reference to the three-dimensional spatial transcriptome coordinates, the dataset is rendered into the SpatioTemporal-Multi-omic Atlas of Gastrulating In-silico Cells (ST-MAGIC). By integrating transcription factor (TF)-target gene-target chromatin region cascades as well as signaling effector chromatin immunoprecipitation sequencing (ChIP-seq) data, we developed the ST-MAGIC(+) platform to explore the spatiotemporal dynamics of TF networks and determine intrinsic responsiveness to crucial developmental signals.

The learning gleaned from the ST-MAGIC and ST-MAGIC(+) atlas can inform the spatiotemporal coordination of GRNs for governing lineage allocation and tissue patterning in the germ layers of the gastrulating embryo. Here, from the ST-MAGIC and ST-MAGIC(+) atlas, we have determined the multi-omic basis for left-right symmetry breaking events in the gastrula, identified and experimentally validated the complex GRNs hierarchy involving TFs (EOMES, FOXA2, NOTO/POU6F1) and signaling pathways (such as NODAL and WNT signaling) that regulate mouse axial mesendoderm development, accompanied by a spatial developmental route from the anterior primitive streak to the distal region of the embryo.

Results
Single-cell multi-omics capture the GRNs underlying germ layer formation

To spatiotemporally dissect the intrinsic GRNs at single-cell resolution during mouse gastrulation, we have undertaken: (1) profiling the transcriptome and chromatin accessibility of individual cells of mouse gastrula at early-streak stage (E6.5) to the late-streak stage (E7.5) at 6 h (0.25 embryonic day) intervals; (2) constructing gene-peak linkages, that connect transcriptomic and epigenomic modalities; (3) integrating the inferred gene-peak linkages with the Geo-seq-based spatiotemporal information to reconstruct a real time-series three-dimensional spatiotemporal multiome atlas; and (4) analyzing the hierarchical logics of GRNs in parallel with signaling responsiveness to gain insights into the mechanistic attributes of the gene regulatory networks underpinning lineage development and embryonic patterning at gastrulation (Fig. 1a).

In total, we have acquired bi-modal omics data from 35,449 single cells across five embryonic stages (12 E6.5 embryos, 10 E6.75 embryos, 6 E7.0 embryos, 6 E7.25 embryos, and 5 E7.5 embryos) with high cell coverage and high data quality after stringent quality control measurements (Supplementary Fig. 1a–e). Uniform manifold approximation and projection (UMAP) was performed for both transcriptomic and chromatin accessibility (Fig. 1b and Supplementary Fig. 1f, g). To annotate the cell types, iterative transcriptomic clustering for cells from each stage was conducted, and specific gene expression pattern for each sub-cluster was refined. Cells with similar gene expression profiles were then grouped and annotated by inferring known gene expression signatures cross-referenced with the published mouse embryo transcriptome atlases3,6,18. The 31 annotated cell types exhibit consistent composition and high transcriptomic correlation with the published single-cell transcriptome gastrulation atlas3 and delineate the diversifying cell type composition during gastrulation with some cell types selectively enriched for earlier (e.g., anterior primitive streak) or later stages (e.g., cardiac mesoderm, etc.) (Supplementary Fig. 1h, 2a–d). Notably, in contrast to one homogeneous notochord cell cluster recognized in the previous report3,6, here we identified four previously unappreciated distinct subgroups within the notochord lineage, which was named as anterior mesendoderm precursors and node precursors at the E7.25 stage, and anterior mesendoderm and node at the E7.5 stage, respectively (Fig. 1b).

The annotated cell identities were then transferred to cells embedded in the snATAC-seq UMAP based on in silico cell-to-cell matches bridged by the common cellular barcodes incorporated during the preparation of the multi-omic sequencing library (Fig. 1b). Generally, the cells collated in this dataset cover the transitions from early pluripotent epiblast cells at E6.5 to fate-committed neuroectoderm, definitive endoderm and various diversified mesodermal subtypes at the E7.5 stage. We also applied the Trajectories Of Mammalian Embryogenesis (TOME) algorithm18 to systematically analyze the relationship among cell types from two adjacent stages and finally generated a clear relationship graph for all the cell types across gastrulation (Fig. 1c). The inferred lineage trajectories are broadly consistent with previous understandings of mouse early embryogenesis3,18,19, providing a developmental phylogeny to pinpoint the molecular dynamics and hierarchy of transcriptomic and epigenetic regulation during gastrulation.

To map the GRNs governing lineage development during gastrulation, we first determined cell type-specific gene expression, chromatin accessibility and motif enrichment. Intriguingly, despite distinct cell types exhibit distinctive transcriptomic features and epigenomic differences, the motif enrichment remained intermingled for cells within the same pedigree (Supplementary Fig. 3a–c). This observation may indicate that the formation of various cell types and the transcriptomic signatures may be regulated by a common set of TFs, but with frequent turnover of epigenomic landmarks.

Next, to capture the relationship between gene transcriptomic status and peak chromatin accessibility, we developed an algorithm, BioCRE, to capture the linkage between expressed genes and candidate regulatory elements, thereby gaining insight into the relationship between gene expression and chromatin peak accessibility. Distinct from existing tools, such as Signac and ArchR, BioCRE harnesses a bi-orientation regression model leveraging multi-omics data at the chromosome level to identify potential CREs (Fig. 1d). Jensen-Shannon Similarity (JSS) score (see “Methods”), which is calculated based on Jason-Shannon Divergence index, has been developed to measure the consistency between cell type transcriptomic features and linked chromatin accessibility diversities. As shown, the JSS score was higher for BioCRE than for Signac and ArchR (Fig. 1e). Moreover, we found that BioCRE-specific gene-peak linkages are enriched in regulating embryo development-related genes (Supplementary Fig. 4a–f), and show higher cell type-specificity as indicated by Silhouette score (Supplementary Fig. 4g, h). In-depth exploration of performance stability against variable sample sizes demonstrated that BioCRE is more robust with differential sampling sizes than Signac and ArchR, as marked by shorter runtime and higher consistency of linkage composition (Supplementary Fig. 4i, j). Cross-validation of predicted gene-peak linkages using publicly available gold standard promoter capture Hi-C dataset further demonstrated the performance superiority of BioCRE over Signac and ArchR (Supplementary Fig. 5a–j). This result indicates that BioCRE is more effective in predicting cell type-specific regulatory elements that modulate target gene expression.

Generally, BioCRE results show that one gene is linked to a median of 5 CREs with a median gene-to-peak distance of 127,175 bp (Supplementary Figs. 4b, 6a, b). Moreover, we observed that a considerable number of distal chromatin peaks (TSS ± 2-500 kb) exhibit higher JSS score (distal high group) than the regulated gene promoters, highlighting the distal regulatory regions may play prominent roles in regulating gene expression (Fig. 1f and Supplementary Fig. 6c). Co-variation of cell type-specific DEGs and BioCRE-linked peaks’ chromatin accessibility faithfully distinguished the identified cell types in both transcriptomic and epigenomic modalities (Fig. 1g and Supplementary Data 1). For example, ten distal peaks were linked to the expression of nascent mesoderm marker, Mesp1 (Fig. 1g). Amongst the distal linked peaks, d1 and d3 peaks are related to the previously known EME enhancer for Mesp120, while the potential regulatory logics for the remaining eight distal peaks were newly identified. To validate the regulatory force for the newly-identified Mesp1 linked distal peaks, we analyzed the chromatin accessibility attribute for one of the newly-identified distal peaks (named Neo_ME) and the EME element (Fig. 1h). Detailed exploration revealed that, apart from the co-accessible patterns in nascent mesoderm, lateral mesoderm, paraxial mesoderm 2, mesoderm progenitors and ExE mesoderm cells, the Neo-ME element was accessible in the caudal epiblast, caudal mesoderm cells as well as axial mesendoderm related cells (node and anterior mesendoderm). In contrast, the EME element was more accessible in paraxial mesoderm 1, cardiac mesoderm, and mesenchyme cells (Fig. 1h). Thus, while both Neo_ME and EME elements regulate Mesp1 expression, there may exhibit distinct usage preferences of regulatory elements among distinct cell types.

To further explore the role of the newly-identified element, we performed enhancer reporter assays, demonstrating that the specific distribution of the Neo-ME element was consistent with Mesp1 expression (Supplementary Fig. 6d). In addition, we established ESC cell lines with specific deletion of the Neo-ME element and EME element, respectively (Supplementary Fig. 6e). As expected, genetic removal of these two elements led to specific downregulation of Mesp1 expression during embryoid bodies (EBs) differentiation (Supplementary Fig. 6e). Thus, the newly-identified Neo_ME element is likely a critical regulatory element responsible for Mesp1 expression.

Together, through single-cell co-profiling of gene expression and chromatin accessibility in the mouse gastrula, combined with the gene-peak linkage capturing strategy (BioCRE), we established a comprehensive multi-omic atlas which captures the cell type-specific GRNs encompassing gene expression, chromatin accessibility of regulatory elements, as well as gene-peak linkages underlying mouse gastrulation.

Spatiotemporal multi-omic landscape of the in-silico gastrulating cells

Recent technological advances have enabled the measurement of gene expression across tissue sections using various spatial transcriptomic strategies21–25. However, most studies covered only the transcriptomic data module and provided limited real 3D spatial context. A comprehensive spatiotemporal multi-omic map that contains genuine spatiotemporal coordinates for embryo tissues remains elusive.

To construct a spatiotemporal multi-omic map that incorporates the spatial coordinates and timestamps of the cells of the gastrulating mouse embryo, we leveraged the published spatiotemporal registered transcriptome dataset of mouse embryonic tissues during gastrulation6,19 and the spatial transcriptome reconstruction algorithm-Tangram26. We refined a pipeline (see Methods) to reconstruct and visualize cell type distribution, single-cell based spatial transcriptome, and spatiotemporally resolved chromatin accessibility information from embryonic cells at the registered spatiotemporal resolution (Fig. 2a, Supplementary Fig. 7a and Supplementary Data 2) in a ST-MAGIC atlas (SpatioTemporal reconstructed Multi-omic Atlas of the Gastrulating In-silico Cells). Beyond spatial distribution, ST-MAGIC can be customized to investigate the spatial distribution of specific cell types, gene expression, and chromatin accessible domains (Fig. 2a). It is noteworthy that different from the previous gene activity score based chromatin accessibility spatial mapping5, the integrative usage of common cell barcodes present in both the transcriptome and chromatin accessibility modules would significantly enhance the spatial accuracy of the ST-MAGIC.

To assess the efficacy of ST-MAGIC, we found that cells were mapped evenly across spatial spots in the reconstructed atlas (Supplementary Fig. 7b), and the global gene distributions across five stages in ST-MAGIC were well correlated with Geo-seq results (Supplementary Fig. 7c, d), JSS scoring also revealed a consistent gene expression architecture between ST-MAGIC and the Geo-seq dataset (Supplementary Fig. 7e). Taking the expression of Sox2 as an example, the expression pattern of Sox2, gleaned from the ST-MAGIC, closely matched the Geo-seq resource and also the experimental results (Supplementary Fig. 7f–h). Histone modification H3K27ac has been usually used as an active histone marker frequently deposited at accessible chromatin regions27. Following this, we selected chromatin regions with region-specific H3K27ac distribution, as revealed in our previous study4, and checked their spatial distributions of chromatin accessibility from the ST-MAGIC. For example, the chromatin region (chr10:63347793-63348989), which locates upstream of the Sirt1 gene and is specially marked by H3K27ac in the anterior epiblast (A) region (Supplementary Fig. 7i), also exhibits anterior epiblast-specific accessible pattern as revealed by the ST-MAGIC (Supplementary Fig. 7i).

We further analyzed the spatial distributions of the embryonic cell types (Supplementary Fig. 8a). Mostly, cell types were assigned to their expected embryonic spatial positions (Supplementary Fig. 8a–c). For example, ectoderm precursors were mapped to the anterior region of the epiblast, while primitive streak cells were especially mapped to the expected posterior region (Supplementary Fig. 8a, c). Interestingly, we observed a clear regionalization of the mesoderm subtypes, particularly in the E7.5 embryos (Fig. 2b and Supplementary Fig. 8d–f). Specifically, the cardiac mesoderm and mesenchyme cells were mostly located in the proximal region, while the paraxial mesoderm 1 cells and the paraxial mesoderm 2 cells were situated in the medial-distal region (Fig. 2b). The spatial distributions for these mesoderm subtypes were validated by determining the related signatural gene expression in vivo and cross-validated in SEU-3D atlas28 (Fig. 2c and Supplementary Fig. 8d–f).

Apart from revealing the formation of cell type-specific spatial territories, ST-MAGIC also unveils the spatial-specific usage of distal linked peaks for gene regulation. For example, Otx2, which is broadly expressed in the epiblast (EPI) and visceral endoderm (VE)29 (Fig. 2d, e), is linked to three spatial types of chromatin peaks: EPI-specific, VE-specific, and universal peaks (Fig. 2f and Supplementary Fig. 8g). To specify, ST-MAGIC revealed that EPI-specific peaks were predominantly located in the epiblast region, while the VE-specific peaks were located in the visceral endoderm region (Fig. 2f and Supplementary Fig. 8g).

ST-MAGIC allows the exploration of specific cell type compositions, gene expression pattern, and chromatin accessibility at any spatial or temporal coordinate across mouse gastrulation (Fig. 2a). Previously, we reported the symmetry breaking event for the left-right body axis first emerges at the late gastrulation stage, manifesting as differential BMP signaling activity and target gene expression in the contralateral proximal mesoderm6. To trace the multi-omic basis for the initiation of left-right asymmetry, we extracted the molecular information from the proximal lateral region of the mesoderm layer using ST-MAGIC (Fig. 2g–l). Consistently, we found that genes related to left-right symmetry breaking also showed laterally biased expression (Supplementary Fig. 9a). Exploration of cell type composition revealed that the mesoderm progenitors, cardiac mesoderm, paraxial mesoderm 1 and mesenchyme cells were over-represented in the right proximal mesoderm region (RPM), while nascent mesoderm, lateral mesoderm, paraxial mesoderm 1 and paraxial mesoderm 2 cells were enriched in the left proximal mesoderm region (LPM) (Fig. 2g, h). Profiling of the top 1000 specific accessible peaks in both RPM and LPM regions revealed that asymmetric levels of chromatin accessibility was discernable between left and right cell type counterparts (Fig. 2i, k, Supplementary Fig. 9b and Supplementary Data 3). To identify the potential biological functions of these asymmetric peaks, we performed gene ontology analyses for these peaks linked genes. Interestingly, BMP signaling pathway-related genes were associated with RPM enriched peaks, while genes related to somitogenesis and heart development were regulated by peaks with higher accessibility in the LPM (Fig. 2j, l). We also traced the emergence of these asymmetric peaks during gastrulation (Supplementary Fig. 9c–n). Intriguingly, both left and right lateral asymmetric peaks became accessible from E6.75 onward (Supplementary Fig. 9c, l), coinciding with the emergence of mesoderm subtypes from the primitive streak cells. For example, two specific peaks, linked with Lefty2 expression, became accessible at the E6.75 nascent mesoderm cells when Lefty2 expression begins, and showed LPM higher distribution at the E7.5 stage when Lefty2 expression is higher on the left side (Supplementary Fig. 9d–g). One of the two peaks is related to a known ASE element30, while the other is a newly-identified Lefty2 regulatory element (Neo_LRE) (Supplementary Fig. 9e, f). Genetic deletion of Neo_LRE in mouse embryonic stem cells showed that Lefty2 expression was severely affected in the knockout (KO) cells, but no severe morphological changes compared to wild-type (WT) controls during gastruloid differentiation in vitro (Supplementary Fig. 9h–k).

Thus, the ST-MAGIC resource established by spatiotemporal reconstruction of the multi-omic atlas to temporal-matched spatial coordinates facilitates in-depth investigation of multi-dimensional molecular architectures of cell types in defined domains of the gastrulating embryo.

ST-MAGIC(+) infers the spatiotemporal turn-over of enhancer regulons

The hierarchical activation of gene regulatory networks by developmental signals and TFs is crucial for mammalian development31. To explore the logics of TF GRNs underpinning gastrulation, we applied the SCENIC+32, a method for the inference of enhancer-driven TF GRNs, to systematically profile the TF distribution, TF target gene expression as well as TF target peak chromatin accessibility across the annotated cell types. Major regulators for germ layer development were recovered (Fig. 3a and Supplementary Figs. 10, 11). Notably, while the distribution of TF expressions showed high cell type-specificity, the presence of TF targets, especially for the TF target peaks, often exhibited shared pattern among closely related neighbors in the same pedigree (Figs. 1c, 3a). For example, in the blood cell lineage, which encompasses hematoendothelial progenitors, blood progenitors 1 and blood progenitors 2, a clear TF expression hierarchy from Etv2, Gata2 to Tal1 was observed. However, the chromatin accessibility for these TF targets remains indistinguishable. This phenomenon supports that Etv2 could function as a priming factor responsible for enhancer opening prior to the hematoendothelial fate commitment (Fig. 3a)33. The sharing of accessible TF target peaks (such as Noto TF) was also detected in cells of axial mesendoderm lineage and definitive endoderm lineage (Fig. 3a), both of which were derived from the anterior primitive streak (Fig. 1c).

To explore the spatiotemporal distribution of these eRegulon imputed by SCENIC+, we developed a pipeline called ST-MAGIC(+), which projects the pre-defined TF target gene sets or target peak sets to the ST-MAGIC atlas through Area Under Curve (AUC) scoring (Fig. 3b). Theoretically, the so-constructed ST-MAGIC(+) atlas should be able to visualize the spatiotemporal turnover of TF GRNs during mouse gastrulation. To validate the fidelity of ST-MAGIC(+) profiling, we first checked the spatial distribution of region-specific H3K27ac modified peaks inferred from the epigenetic landscape of the mouse gastrula4. As shown, region-specific peak sets were accurately mapped to their sampling origins (Supplementary Fig. 12a). Next, we examined the spatial distribution of GRNs for the major regulators identified in Fig. 3a. Taken the transcription factor, T, as an example, Geo-seq and ST-MAGIC revealed that T expression is localized to the primitive streak and anterior mesendoderm in E7.5 embryo (Fig. 3c), and T’s target genes showed a similar distribution pattern (Fig. 3c). However, for the target peaks of T, we observed a broader spatial distribution pattern than the TF expression, manifesting as high level of chromatin accessibility in both the primitive streak and neighboring mesodermal regions (Fig. 3c). This observation was verified through determining the direct T binding peak around the gene locus of Tbx6, which is a well-known T target gene, using ST-MAGIC (Fig. 3d). Moreover, considering that T continues to express in the primitive streak throughout gastrulation and that the mesoderm cells are the immediate progeny of the primitive streak cells34, the broad accessible pattern of T binding chromatin regions suggests that T may act as a priming factor with the ability to open a broad spectrum of chromatin regions and instruct the subsequent fate commitment of mesoderm cells.

Finally, we dissected the dynamics of shared accessible patterns among trajectory neighborhoods for the axial mesendoderm lineage. As resolved by ST-MAGIC(+), we found that the target genes of NOTO exhibit a matched spatiotemporal distribution with the TF expression, but the chromatin accessibility of NOTO target peaks was widespread in the endoderm layer of the E7.5 embryo, where Noto is not expressed (Fig. 3e and Supplementary Fig. 12b, c). Further examination revealed that the chromatin accessibility of NOTO target peaks was elevated at early gastrulation stage, especially in the E7.0 anterior primitive streak region, well before Noto expression emerged (Fig. 3e; Supplementary Fig. 12b). Detailed exploration of the direct NOTO target chromatin region around the Noto locus supported the spatiotemporal kinetics of NOTO eRegulon by ST-MAGIC(+) (Fig. 3f). Global examination of the spatiotemporal concordance of TFs and their targets by Pearson correlation and SSIM index highlighted that the pervasive spatiotemporal asynchrony between TFs and their targets (Supplementary Fig. 12d, e).

These results indicate that crucial upstream players would be involved in setting up the chromatin level of TF GRNs, such as NOTO, during axial mesendoderm lineage development. It has been reported that the modulations of WNT and NODAL signaling are involved in the patterning of the anterior primitive streak35. We found that both WNT and NODAL signaling are highly enriched at these chromatin regions (Fig. 3g), suggesting crucial inductive signals may participate in establishing the pre-accessible chromatin pattern for NOTO targets.

ST-MAGIC(+) reflects the spatiotemporal responsiveness to developmental signaling

Inductive developmental signals such as WNT and NODAL signaling act as morphogen cues with localized production but imparting distant function, playing crucial roles in elaborating the proper embryo arrangement and related cell lineages. Multi-layered dimensions, including the morphogen dosage (concentration and duration), the completeness of the signaling cascade components, and the intrinsic competence of the receiving cells, collectively influence the interpretation of signals. Currently, a systematic evaluation of the intrinsic chromatin competence for signaling function remains unexplored36.

To characterize the spatial context-dependent mechanisms for interpreting developmental signals, we applied ST-MAGIC(+) to investigate the hierarchy of signal and effector distribution, associated chromatin responsiveness and gene-peak linkage inferred transcription output along the lineage trajectory (Fig. 4a). The spatial distributions of key signal components were characterized by determining the expression pattern of related genes, the chromatin responsiveness to signals was delineated by profiling the spatial accessibility of direct signal effector targets, and the spatial transcription output was measured by checking the expression pattern of signal effector chromatin binding regions linked genes (Fig. 4a).

We first analyzed the global pattern and the cellular responsiveness of WNT signaling across gastrulation (Supplementary Fig. 13a). Consistent with known biology37, the Wnt signal ligand encoded by Wnt3 was predominantly expressed in the primitive streak and adjacent mesoderm tissue, while the Wnt signal antagonist, Dkk1, was specifically expressed in the anterior visceral endoderm region (Fig. 4b and Supplementary Fig. 13b–e). However, the WNT signaling receptor, Lrp6, and the effector, Ctnnb1, (also known as β-Catenin) were widely distributed within the gastrula (Supplementary Fig. 13f).

To systematically characterize the embryo chromatin responsiveness to WNT signaling, we sourced two published ChIP-seq datasets for WNT signaling effector-β-Catenin from the in vitro ESCs38 and EpiLCs39. Through differential binding peaks analyses, the chromatin binding regions were classified into three groups of common, ESC-specific, and EpiLC-specific peaks (Supplementary Fig. 13g, h). Spatiotemporal reconstruction of these peaks though ST-MAGIC(+) revealed that ESC-specific peaks were accessible in the anterior epiblast, where cells remain in a pluripotent state40 (Supplementary Fig. 13i). Meanwhile, for the peaks with EpiLC-specific β-Catenin binding, ST-MAGIC(+) reported that these peaks were accessible in the embryonic region, largely overlapped with Wnt3 distribution (Fig. 4b, c). In-depth analyses of the ST-MAGIC(+) results for EpiLC-specific peaks across five stages of the gastrula showed that the PS region and neighboring posterior mesoderm region remained consistently accessible across the gastrulation (Fig. 4b, c and Supplementary Fig. 13j, k). Together, these results underscore that in response to WNT signaling input, cells first turnover the Wnt-responsive chromatin landscape during the transition from pluripotent state to lineage primed progenitors, and then modulate the chromatin accessibility to accommodate the specification of posterior mesoderm subtypes during gastrulation (Supplementary Fig. 13l).

We next analyzed the intrinsic cellular responsiveness to NODAL signaling (Supplementary Fig. 14a), a crucial signal for embryo patterning and germ layer formation41. Examination of the components of NODAL signaling cascade in the ST-MAGIC atlas revealed that the ligand-Nodal was initially expressed in the posterior epiblast at early gastrulation stage but abruptly shifted to the distal tip region at the E7.5 stage (Fig. 4d). In the meantime, the NODAL signaling antagonist-Cer1 was constantly expressed in the anterior visceral endoderm region (Fig. 4d and Supplementary Fig. 14b–f). However, the NODAL signaling receptors and the effectors were widely expressed in the gastrula (Supplementary Fig. 14g). The changes of expression domain of the NODAL ligand-antagonist pair may underpin the shift of the NODAL morphogen gradient during mouse gastrulation.

To chart the chromatin responsiveness along with the NODAL pattern shift, we incorporated the published ChIP-seq datasets for NODAL signaling effectors SMAD2 and SMAD342. We found that the chromatin binding ability of the Smad2/3 complex was largely de novo generated in EB cells, an in vitro counterpart of mesoderm cells (Supplementary Fig. 14h, i). GREAT analyses further indicated that these peaks are involved in the cell fate specification, anterior/posterior formation and primitive streak formation (Supplementary Fig. 14i). Then, we used ST-MAGIC(+) to profile the spatial distribution of these peaks during gastrulation. Interestingly, we found that the global chromatin accessibility of these EB-specific Smad2/3 binding peaks first appeared in both the endoderm region and PS regions (Fig. 4e). Subsequently, the accessibility of these peaks in the PS region gradually shifted and formed a distal-to-proximal accessibility gradient in the E7.5 gastrula (Fig. 4e and Supplementary Fig. 14j). Notably, the chromatin responsiveness to Nodal signaling started relocating at an earlier stage between E6.75 and E7.0, ahead of the re-alignment of the Nodal gradient in the germ layers occurring between E7.25 and E7.5 stage (Fig. 4f, g). These results highlight that the spatial turnover of NODAL signaling chromatin responsiveness occurs earlier than the repositioning of Nodal activity.

To determine the intrinsic features for the relocation of signaling chromatin responsiveness, we extracted the specific subset of signaling responsive chromatin peaks which show high accessibility at the E6.75 primitive streak region (E6.75_PS) or at the E7.0 anterior primitive streak and adjacent mesoderm region (E7.0_APS_M) (Fig. 4h). Generally, the signaling responsive peaks for NODAL and WNT signaling were largely independent between E6.75_PS region and E7.0_APS_M region (Fig. 4h and Supplementary Data 4, 5). As shown, the E7.0_APS_M specific WNT signaling responsive peaks were gradually extended to the whole posterior region (Fig. 4i), and enriched with typical WNT signaling co-effector LEF1 and TCF7L1 motifs (Fig. 4j). The posterior enriched WNT signaling responsiveness (Fig. 4k) support the notion that WNT signaling is instrumental for regulating posterior embryo development36. Interestingly, for NODAL signaling, profiling of the E7.0_APS_M specific signal responsive peaks exhibited a gradual shift pattern of these peaks towards the distal region of the E7.5 embryo, where Noto and Nodal-expressing cells reside (Figs. 3e, 4d, l). Functional characterization of E7.0_APS_M NODAL signaling responsive peaks showed enrichment in chordate embryonic development and embryo pattern specification process, and the knockout of related genes can frequently lead to abnormal germ layer morphology, abnormal rostral-caudal axis patterning and absent floor plate (Supplementary Fig. 14k). Motif enrichment analyses of these peaks revealed significant enrichment of major TF regulators especially for axial mesendoderm development such as Zfp281, Foxa2, and Foxj1 (Fig. 4m). These results indicate that the spatial re-distribution of embryonic responsiveness to Nodal signaling may play specific roles in the forthcoming development of axial mesendoderm lineage (Fig. 4n).

Molecular hierarchy underlying axial mesendoderm lineage development

The axial mesendoderm has been reported to be the direct developmental antecedent for the midline notochord cells, which instruct the following somitogenesis and neural patterning7. To better understand axial mesendoderm lineage development, it is essential to document the spatiotemporal context of the stepwise appearance of molecular features and cellular states during the lineage formation process, with accurate spatiotemporal information.

Here, based on the inferred lineage trajectory for the gastrula (Figs. 1c, 5a) and the ST-MAGIC reconstructed cell distribution for the mouse gastrula (Supplementary Fig. 8a), we found that cells residing in the E7.5 distal tip, where the prospective notochord cells first emerge43, are composed of the node cells and anterior mesendoderm cells (Fig. 5a and Supplementary Fig. 15a). Expression of known notochord markers, Shh, Noto, and Foxj1, revealed that the presence of two distinct cell subtypes (Shh+Noto-Foxj1- and Shh+Noto+Foxj1+) at the ventral distal surface for the E7.5 embryo, and the two distinct cell subtypes can persist through organogenesis (Fig. 5b and Supplementary Fig. 15b–d).

The inferred lineage trajectory indicates that the node and the anterior mesendoderm cells are spatially derived from the E7.0 APS cells (Fig. 5a). Concurrently, chromatin responsiveness to NODAL and WNT signaling gradually increased in the E7.0 anterior primitive streak cells and is maintained at high accessibility level till the E7.5 stage (Fig. 5c). To chart the molecular dynamics along the developmental trajectory, we applied CellRank44 to assign the cell fate probabilities for all the related cells at single-cell resolution. RNA velocity and pseudotime ordering revealed that the cells at an early stage (E6.5 primitive streak, E6.75 anterior primitive streak) are unspecified for cell fate directions, while cells at more advanced stages show directed flow towards the E7.5 node and E7.5 anterior mesendoderm (Fig. 5d). In-depth analyses revealed that the two cell subtypes showed distinct transcriptomic and epigenomic features (Supplementary Fig. 15e). To specify, the node cells express high levels of cilia and dynein-related genes, such as Foxj1 and Dnah11, while the anterior mesendoderm cells express high level of collagen-related genes, such as Col2a1 and Srd5a2 (Fig. 5e and Supplementary Fig. 15f–i). Molecularly, smoothed trends of gene expression along the pseudotime order showed that cilium related genes, such as Dynlrb2 and Rsph9, were gradually upregulated down the road to node, while the collagen gene Col2a1 was expressed till the terminal stage of the anterior mesendoderm (Fig. 5f, g).

Pseudotime-based tracking of gene expression cascades and associated chromatin peaks revealed three successive stages (G1, G2, G3) enriched in the early, intermediate and terminal populations, respectively (Fig. 5h–k). Intriguingly, for both lineages, a subset of the terminal stage expressed genes (G3a) was regulated by chromatin peaks which get accessible at an earlier stage (Fig. 5h–k). For example, in node cells, Foxj1 expression was evidently upregulated at the E7.25 stage, but the chromatin accessibility for the linked distal peak-Foxj1-DRE was accessible in the E7.0 anterior primitive streak cells (Fig. 5i). Similarly, for anterior mesendoderm cells, genes like Sox9 showed the same pattern (Fig. 5k). The presence of these pre-accessible chromatin peaks points to early epigenomic priming for the development of axial mesendoderm lineage.

To systematically determine the molecular hierarchy of the enhancer-driven TF GRNs for the axial mesendoderm lineage, we used SCENIC+ to infer the potential candidate TFs and the dynamics of downstream target genes and linked chromatin peaks (Fig. 5l). A clear hierarchy of TF usage turnover was captured. To specify, pluripotency-related TFs, such as Nanog, were enriched at the E6.5 stage; early lineage factors, such as Eomes, were enriched at the E6.75 stage; subsequently, for the intermediate stage factors, such as Mixl1, Foxa2 and Lhx1, were enriched in the E7.0 anterior primitive streak cells, E7.25 node precursors, and E7.25 anterior mesendoderm precursors, respectively; and terminal stage factors, such as Noto and Sox9, were abundant in E7.5 node and anterior mesendoderm cells, respectively (Fig. 5l). Detailed analyses of these eRegulons revealed that the early and intermediate TFs first upregulate TFs gene expression, followed by increasing the target gene expression as well as target peak chromatin accessibility (Fig. 5m and Supplementary Fig. 16a). In contrast, for the terminal stage enriched TFs, the chromatin of TF target peaks becomes accessible at an earlier stage, followed by TF expression and TF target gene expression (Fig. 5n). This temporal turnover of TF GRNs may reflect a molecular relay underlying the sequential cascade of cell fate commitment during axial mesendoderm development. Notably, among multiple known TF GRNs, we identified a novel TF, POU6F1, which exhibits highly specific expression of both the TF itself (adjusted P-value, 8.4e-29) and its target genes (adjusted P-value, 1.3e-55) in the node (Fig. 5l and Supplementary Fig. 16b). However, the chromatin regions targeted by POU6F1 get accessible prior to the expression of the TF and its target genes (Fig. 5n).

Distinct roles of stage-related TFs in regulating gene expression and setting up chromatin accessibility

To determine the relevance of TF hierarchy and signal responsiveness during axial mesendoderm lineage development, we first checked the distribution of the TF GRNs for early stage TF (EOMES), intermediate stage TFs (MIXL1, FOXA2, LHX1), and terminal stage TFs (NOTO and POU6F1) during gastrulation using ST-MAGIC(+). Consistent with the SCENIC+ results (Fig. 5m, n), we found that for the early and intermediate stage TFs, expression of these TFs commences prior to or concurrently with the target gene expression and target peak accessibility (Supplementary Fig. 16c–f). In contrast, for the terminal stage TFs-NOTO and POU6F1, we observed that their target chromatin peaks first become accessible at the expected APS region and endoderm region at the E7.0 stage, when TFs and their target genes only show minimal expression (Figs. 4e, 6a, b and Supplementary Fig. 17a–g). With the progressing of embryo development, once Noto and Pou6f1 reached the expression summit at the E7.5 node region, the target genes show abundant enrichment at the same region (Figs. 4e, 6a, b and Supplementary Fig. 17a–g). These results suggest that the early and intermediate stage TFs may function as priming factors by opening a broad spectrum of chromatin regions, thereby setting up a permissive chromatin environment for the following signaling and the terminal stage TFs-NOTO and POU6F1. The expression of TFs-NOTO and POU6F1 at the terminal stage can then interact with the pre-accessible chromatin regions, and enhance the expression of target genes to establish the requisite transcriptomic state for further lineage differentiation.

To demonstrate this molecular framework, we integrated datasets of both EOMES ChIP-seq (from embryoid bodies in vitro) and FOXA2 ChIP-seq (from mesendoderm cells differentiated in vitro) and chromatin accessibility data of Eomes-KO and Foxa2-KO cells. By profiling the enrichment around the pre-opening Foxj1-DRE element (Fig. 5i), we found that the chromatin accessibility for Foxj1-DRE element was markedly downregulated in Eomes-KO and Foxa2-KO cells (Supplementary Fig. 18a), which strongly suggests the roles of EOMES and FOXA2 in opening Foxj1-DRE chromatin region. We also systematically analyzed the enrichment of EOMES and FOXA2 around the E7.0_APS_M signaling responsive peaks, the NOTO target peaks, and also the POU6F1 target peaks, we found that both EOMES and FOXA2 were enriched around these pre-defined genomic regions (Fig. 6c and Supplementary Fig. 18b, c). Moreover, the knockout of either Eomes or Foxa2 led to the reduction of chromatin accessibility around these loci (Fig. 6c and Supplementary Fig. 18b–d). Therefore, EOMES and FOXA2 may play roles in pre-opening signaling responsive elements and setting up chromatin accessibility.

Next, we investigated whether the TFs-NOTO and POU6F1 play roles in establishing the terminal transcriptomic state but not the chromatin setup. Cross-referencing published atlas3 confirmed the enrichment of these two genes in the notochord cells (Supplementary Fig. 19a). We then generated two mouse mutants with the genetic deletions of the Noto or Pou6f1 genes (Fig. 6d and Supplementary Fig. 19b). No visible morphological phenotypes were detected in either mutants at E7.5 stage (Supplementary Fig. 19c). To probe the potential molecular abnormalities, we collected E7.5 WT control, Noto KO, and Pou6f1 KO embryos and performed 10x scMultiome sequencing (snRNA-seq+snATAC-seq). Cells were annotated by computationally mapping their transcriptome onto our E7.5 WT atlas (Fig. 6e). The cell type compositions in the KO embryos remained comparable to WT, except for the axial mesendoderm-related cells, in which Noto and Pou6f1 are expressed and now successfully removed (Fig. 6f and Supplementary Fig. 19d, e). Moreover, in both Noto KO and Pou6f1 KO embryos, the expression of node cell-related genes and ciliogenesis was severely disrupted (Fig. 6g–j, Supplementary Fig. 19f, g and Supplementary Data 6). Phenotypically, by E11.5, directional axis turning, which is related to proper notochord function17, was randomized (Supplementary Fig. 19h, i). Importantly, the chromatin accessibility of the TF target chromatin peaks remained unchanged (Fig. 6k and Supplementary Fig. 19j). These results demonstrate that NOTO and POU6F1 function primarily as the final-step transcriptional regulators in the TF relay during the sequential development of the mouse axial mesendoderm.

Together, through systematic exploration of the ST-MAGIC (+) resource, we unraveled the cellular events and associated highly-organized regulatory cascades of TF and signal GRNs, and demonstrated the differential roles for TFs from different stages in shaping the transcriptomic architecture and chromatin accessibility landscape at different milestones of the lineage trajectory (Fig. 7).

Methods
Ethics statement
Mouse used in this study were housed and bred in the SPF facilities of Guangzhou National Laboratory. All animal experiments were performed in compliance with the guidelines of the Animal Core Facility.

Mouse embryo collection and multi-omic profiling
For embryo sampling, C57BL/6J embryos were harvested from pregnant mice at day 6.5, 6.75, 7.0, 7.25, and 7.5 of gestation (day of vaginal plug detection = Day E0.5). Plugged female mice were picked after mating and marked as embryonic day 0.5 (E0.5). Female mice were sacrificed for embryos collection at specific gestational stages. Embryos were isolated from the uterus and carefully transferred into pre-cool PBS in petri dishes, and surrounding decidua and parietal endoderm tissues were removed using needles under an Olympus stereoscope. Sex of the embryos was not considered in this study.

Careful morphological staging of the acquired embryos was performed before single-cell multi-omic data preparation. The developmental timepoints of embryos were staged by the proximal-distal span of the PS and the anterior-posterior span of the mesoderm layer. To meet the loading requirements for 10x genomics, embryos from the same stage were pooled together and subjected to TrypLE Express enzyme (Gibco, 12604) incubation at 37 °C for 7–10 min. The acquired single-cell suspension was carefully washed and filtered to ensure proper integrity and avoid cell clumping of each single cell. Cell counts were then assessed with a haemocytometer counted under a microscope. Nuclei isolation and multiome library preparation were performed by following the manufacturer’s instructions (https://www.10xgenomics.com/cn/support/single-cell-multiome-atac-plus-gene-expression/documentation).

For the single-cell multi-omic profiling of the Noto KO and Pou6f1 KO embryos, the parents of the heterozygotic mutation were crossed and checked for the vaginal plug. E7.5 embryos were freshly harvested, and a tiny portion of extra embryonic tissues were collected and genotyped. During the genotyping process, the embryos were freshly frozen in liquid nitrogen. Embryos with the same genotypes (WT, Noto KO homozygote, Pou6f1 KO homozygote) were grouped and dissociated into single nuclei. Approximately 15,000 nuclei for each group were collected and loaded for further single-cell multiome library preparation, following the manufacturer’s instructions.

Enhancer activity reporter assay
Generally, the enhancer reporter assay was performed as previously reported4. In brief, DNA fragments for potential regulatory elements were cloned from the C57BL/6J mouse genome and then ligated into the plasmid construct containing the minimal Hsp68 promoter and LacZ. The acquired purified plasmids were then linearized and used for pronuclear injections of PN4 stage zygotes with a FemtoJet Microinjection System (Eppendorf). The injected embryos were cultured to the 2-cell stage in KSOM medium with amino acids at 37 °C under 5% CO2, and then transferred to the oviduct of pseudo-pregnant ICR females and marked as 0.5 dpc. Embryos were collected at the corresponding stage for LacZ staining. LacZ staining was performed using a commercialized β-gal staining kit (Beyotime, RG0039).

RNAscope staining for whole mount or sections of embryo
RNAscope probes including mm-Lefty2-C2 (436291-C2), mm-Hand1-C1 (429651), mm-Irx5-C2 (513871-C2), mm-Isl1-C3 (451931-C3), mm-Sox2-C1 (401041-C1), mm-Otx2-C3 (444381-C3), mm-T-C3 (423511-C3), mm-Eomes-C2 (429641-C2), mm-Foxa2-C4 (409111-C4), mm-Shh-C2 (314361-C2), mm-Noto-C3 (1253281-C3), mm-Foxj1-C3 (317091-C3), mm-Dynlrb2-C3 (1243011-C3), mm-Col2a1-C4 (407221-C4), mm-Rsph9-C2 (430201-C2), and mm-Pou6f1-C1 (801931-C1) were bought from the Advanced Cell Diagnostics.

For whole mount RNAscope staining, embryos were fixed in 4% PFA overnight. After serial dehydration and rehydration of embryos using gradient methanol-PBS solution, embryos were subject to the RNAscope protocol following the manufacturer’s instruction (https://acdbio.com/sites/default/files/MK%2050016%20Zebrafish_WISH_Tech%20Note_12042017.pdf). Images of stained embryos were acquired using the LiTone XL system (Light Innovation Technology Limited).

For RNAscope staining with embryo sections, dissected embryos were first fixed in 4% PFA, then dehydrated in 20% Sucrose-PBS and 30% Surcrose-PBS solution, respectively. The dehydrated embryos were embedded in OCT matrix and then cryo-sectioned using Leica CM1950. RNAscope workflow was performed by following the manufacturer’s instructions (https://acdbio.com/ebook/introduction/materialsmethod). Images were acquired by using the Carl Zeiss LSM980 system.

Wholemount in situ hybridization
Wholemount in situ hybridization of RNA transcripts was performed by following the published protocol61. Briefly, DNA fragments encoding the probes of Lety2 were firstly PCR amplified by using the oligos in Supplementary Data 7 using a mouse embryo cDNA library. Embryos at relevant stages were collected in DMEM media and then fixed in 4% PFA at 4 °C overnight. Fixed embryos were then washed in PBS solution at room temperature to remove residual PFA, followed by serial dehydration and rehydration in PBS, 25% Methanol/PBS, 50% Methanol/PBS, 75% Methanol/PBS, and 100% Methanol. Afterward, the embryos were treated with 10 μg/mL proteinase K and incubated with DIG-labeled RNA probes at 70 °C overnight. To remove the remaining RNA probes, embryos were washed in TBST buffer with frequent buffer changing for at least 2–4 h. Then the embryos were incubated with anti-DIG-AP antibody at 4 °C overnight. The embryos were subjected to sufficient TBST wash before final staining with NBT and BCIP solution. The final images of stained embryos were collected using an Olympus SZX16 microscope.

The generation of genome edited mouse embryonic stem cell or mouse model
In this study, we have prepared the genome-edited mouse embryonic stem cell lines with Mesp1 Neo_ME element, Mesp1 EME element, and Lefty2 Neo_LRE element knockout, respectively. The CRISPR-Cas9 system was used to generate the genome-edited mouse embryonic stem cell lines. Guided RNAs were designed by using the online tool Chop-chop (http://chopchop.cbu.uib.no/). The synthesized sgRNA DNA fragments were then ligated into the px330 plasmid with a Cas9 protein expression cassette. The acquired purified plasmids were then transfected into the WT cells using Lipofectamine 3000. The genome-edited clones were then selected and genotyped as previously reported62. The potential off-target editing sites were tested based on the website’s indication.

To determine the roles of Pou6f1, we also generated a mouse model with Pou6f1 deletion. The specific sequences for sgRNA pairs were included in Supplementary Data 7. DNA fragments for the sgRNA and T7-Cas9 were transcribed and purified in vitro using MMESSAGE MMACHINE T7 Ultra Kit (Invitrogen, AM1345) and MEGAclear kit (Invitrogen, AM1908). To prepare a sufficient number of fertilized zygotes, C57BL/6J female mice (4 weeks old) were superovulated and mated with the male C57BL/6 mice. Twenty-four hours later, fertilized embryos were collected from the oviducts. RNA for Cas9 protein (100 ng/µl) and corresponding sgRNA (100 ng/µl) were mixed in HEPES-CZB medium containing 5 μg/ml cytochalasin B (CB) and injected into the cytoplasm of fertilized eggs using a FemtoJet microinjector (Eppendorf) with constant flow settings. The injected embryos were cultured in KSOM with amino acids at 37 °C under 5% CO2 in the air to reach the 2-cell stage after 24 h in vitro. Two-cell embryos were transferred into pseudo-pregnant ICR female mice. The acquired mouse individuals were then subjected to genotyping for successful Pou6f1 deletion. Mice with Pou6f1 deletion were retained and recorded as F0. To prepare the stable Pou6f1 deletion mouse line, F0 mice were crossed, and the acquired F1 mice were checked. Pou6f1 KO heterozygotes F1 mice were maintained, and the population was expanded for the following experiments.

The Noto KO mouse line was purchased from GemPharmatech (Strain ID: T017401).

Scanning electron microscopy
Embryo samples were freshly collected and fixed in 2.5% glutaraldehyde overnight at 4 °C. Subsequently, the samples were thoroughly washed with PBS four times at 10-minute intervals to ensure the complete removal of the fixative. Fixed embryos were treated with osmium tetroxide for 1 h to enhance the contrast and stability of cellular structures under the electron beam. The samples were then washed with ddH2O five times at 10 min intervals to remove any residual osmium tetroxide. Dehydrate the embryos in an ethanol series for 5–10 min each: 50%, 70%, 85%, and three times in 100%. Transferred the embryos in ethanol to baskets for critical point drying (CPD) in a critical point dryer machine (Quorum K850 critical point dryer). After that, the samples were mounted on sample stubs and sputter-coated with gold using a Quorum Q150RS sputter coater. Finally, the samples were imaged using the Zeiss GeminiSEM 300 scanning electron microscope.

The ST-MAGIC atlas
Pre-processing of Geo-seq and 10x scMultiome: The Geo-seq data utilized in this study were obtained from our previously published datasets, specifically including E6.5, E7.0, and E7.5 from GSE120963, and E6.75, E7.25, and E7.5 (left-right resolved) from GSE171588. The five stages of 10x scMultiome data were newly generated in this study. For the Geo-seq data, the quantification of the FASTQ files was performed using RSEM (v1.3.3), leveraging annotations from the Cell Ranger GTF file. The normalized TPM values were processed using the ‘NormalizeData’ function in the Seurat (v4.1.1) package and were then used for further analytical procedures. For the scMultiome data, the normalized snRNA-seq and snATAC-seq data were processed using Seurat and Signac (v1.6.0) and were subsequently employed for further analytical procedures.

Alignment of 10x scMultiome to Geo-seq: The alignment of snRNA-seq data to Geo-spatial spots was achieved using Tangram (v 1.0.4), with the precise spatial location of each cell determined from the mapping matrix with the highest probability. The gene expression of each reconstructed spot was obtained by averaging the gene expression levels of the cells mapped to the indicated spot, and the gene expression percentage of a spot represents the proportion of cells within the indicated spot that express the gene (gene reads > 0). For snATAC-seq data mapping, the spatial location of each cell from snRNA-seq is shared. The peak accessibility of each reconstructed spot was obtained by averaging the peak accessibility level of the cells mapped to the indicated spot, and the peak accessibility percentage of a spot represents the ratio of cells within the indicated spot with an open peak (peak reads > 0). The cell type distribution was obtained using the tg.project_cell_annotations function within the Tangram toolkit, and the cell type percentage of a spot represents the ratio of cells of that type within the spot.

In the application of Tangram, the training genes are specified as the zipcode genes relevant to each stage, and the mapping strategy is customized for individual cells. For the E6.5 and E6.75 embryos, the top 100 cell-type marker genes identified by the COSG (v 0.9.0) package63 were designated as zipcode genes. For the E7.0, E7.25, and E7.5 embryos, the PC loading genes calculated by the Seurat package (E7.0, top 500 highest and 500 lowest genes of PC1-10; E7.25, top 50 highest and 50 lowest genes of PC1-10; E7.5, top 500 highest and 500 lowest genes of PC1-5) were designated as zipcode genes. The zipcode list has been provided in Supplementary Data 2.

Since Geo-seq primarily captures the proximal to distal embryonic regions of the gastrula, cell types from the scMultiome data that is not present within these zones were excluded from the analysis. Specifically, the following cell types were filtered out: ExE endoderm, ExE ectoderm, Haematoendothelial progenitors, Blood progenitors 1, and Blood progenitors 2.

Smooth of corn plot: The smoothing of mapped features in ST-MAGIC is achieved through Gaussian weighting of neighboring spots, including the mapped feature smooth denoted as F and the percent smooth denoted as P, defined as:
F_smooth = Σ n j=1 w i,j (θ) × F j, w i,j (θ) = e - ||Fi - Fj||2 / 2θ2
P_smooth = Σ n j=1 w i,j (θ) × P j, w i,j (θ) = e - ||Pi - Pj||2 / 2θ2

For the selection of neighboring spots, ST-MAGIC provides three types of input: column neighboring spots along the proximal-distal axis by default, row neighboring spots along the anterior-posterior axis, and customized neighboring spots. The query is designed to recognize user-defined consistent patterns, whether they are genes, peaks, or cell types. The spatial pattern of interest was defined in binary, with the value of the corresponding spot set to 1 and the value of all other spots set to 0. Similar omics data matching this binary pattern were selected using cosine similarity, as implemented in the COSG package. For example, ST-MAGIC detects the specified peaks in customized RPM and LPM regions (distal peaks in the left and right proximal mesoderm regions shown in Fig. 2g, h). This process will facilitate the identification of spatiotemporally regulated patterns in individual cells.

The ST-MAGIC(+) atlas
Calculation of enrichment score: Building on the single gene or peak analysis capabilities of ST-MAGIC, ST-MAGIC (+) extends these features by enabling the spatial monitoring of gene sets and peak sets from SCEINIC+ or public resources through enrichment analysis by the AUCell (v1.18.1) package. The enrichment of gene sets and peak sets of each spot was obtained by averaging the AUC score of the cells mapped to that spot, and the enrichment percentage of a spot represents the average ratio of expressed genes or accessible peaks within that spot.

Spatial signaling responsiveness pattern: The binding of signaling molecules to multiple genomic regions highlights the complexity and interconnectedness of gene regulatory networks. It highlights how single factors can coordinate with diverse downstream targets to orchestrate cellular responses and developmental programs. ST-MAGIC(+) enables the depiction of the spatial activities of signaling responsiveness, such as WNT and Nodal (Fig. 4), using data from public ChIP-seq datasets. Before calculating the AUC score, the peaks from public datasets must first be aligned to those identified in our own dataset using the GenomicRanges (v1.48.0) package to ensure a consistent genomic context for comparative analysis.

BioCRE model
The construction of BioCRE: Unlike published tools that detect CREs using only a single peak or a single gene, BioCRE employs a systematic approach to calculate these linkages. BioCRE harnesses a bi-orientation regression model leveraging multi-omics data decoding GRN at the chromosome level to identify potential CREs.

To formulate the linear optimization procedure using paired snRNA-seq and snATAC-seq data, let Gθ ∈ Rg×n denotes snRNA-seq data profiling g genes across n cells on a given chromosome θ, let Pθ ∈ Rp×n denotes snATAC-seq data with p peaks for the same set of n cells on the chromosome θ. We assumed the linkage matrix Lθ ∈ Rg×p is the cis-regulation between genes and peaks on chromosome θ. Our goal is to learn the linkage matrix L that captures cis-regulation under two conditions: (1) the predicted snRNA-seq matrix G^θ, generated using the linkage matrix Lθ and snATAC-seq matrix Pθ, should be similar to the real snRNA matrix Gθ; (2) the predicted snATAC-seq matrix P^θ, generated using the linkage matrix Lθ and snRNA matrix Gθ, should be similar to the real snATAC-seq matrix Pθ. We minimize the following cost function as:
arg min Lθ ∈ {chr1, chr2, ..., chrY} (E ||Gθ - G^θ||2 + E ||Pθ - P^θ||2 + λθ × REG loss)

The cost function consists of three parts:
1) The expected loss between the real snRNA-seq matrix Gθ and the predicted snRNA-seq matrix G^θ:
E ||Gθ - G^θ||2 = Σg i=1 Σn k=1 (Gθ i,k - G^θ i,k)2 = Σg i=1 Σn k=1 (Gθ i,k - Σp j=1 Lθ i,j × Pθ j,k)2
Where Gθ i,k is the expression level of gene i in cell k on chromosome θ, Pθ j,k is the chromatin accessibility of peak j in cell k on chromosome θ, the Lθ i,j is the strength of linkage between gene i and peak j on chromosome θ.
2) The expected loss between the real snATAC-seq matrix Pθ and the predicted snATAC-seq matrix P^θ:
E ||Pθ - P^θ||2 = Σp j=1 Σn k=1 (Pθ j,k - P^θ j,k)2 = Σp j=1 Σn k=1 (Pθ j,k - Σg i=1 Lθ i,j × Gθ i,k)2
3) The L2 regularization loss of the linkage matrix Lθ with hyperparameter λθ:
REG loss = ||Lθ||2 = sqrt(Σg i=1 Σp j=1 (Lθ i,j)2)

Under this model, we systematically considered all combinations of genes and peaks within the same chromosome, effectively simulating the complexity of GRN.

To mitigate the risk of gradient vanishing in a model with a large number of parameters, ‘He’ initialization was employed for the weight matrices:
W ∼ N(0, 2/p)
Where p is the number of peaks. The L-BFGS algorithm is employed to achieve the minimization of the cost function, the training process is configured for 500 epochs with a learning rate of 0.5.

After the optimization of L, the extraction of significant linkages is required. Firstly, we computed bi-orientation p-values, which involved assessing the significance of a given CRE in relation to all CREs linked to the same gene by z-score, followed by evaluating its significance among all CREs associated with the same peak, the p-values were adjusted using the Bonferroni correction method applied separately to each set of comparisons. Then, the two p-values are combined into a single value through Pearson’s method in the ‘SciPy’ package:
-2 Σn i=1 log(1 - pi)
Finally, the CREs are constrained to those within a 1 megabase range, and only those with a combined p-value less than 10-7 are reserved.

Benchmark analysis of BioCRE
Performance on gastrula multi-omic dataset: In the absence of a genome-wide gold standard for gene-to-peak linkages in the mouse gastrula, we used the Jensen-Shannon Similarity score, which is calculated based on the Jensen-Shannon Divergence index, to evaluate the concordance between genes and peaks in linkage for BioCRE, compared to linkages established using default parameters in ArchR and Signac. The JSS score of each linkage was calculated using the normalized gene expression Gi from snRNA data and the peak accessibility Pj from snATAC data across n cells.
JSS = 1 - E(G, P) = 1 - (sqrt(E(G) + E(P)) - sqrt(E(G + P))) / 2
E(P) = - Σn k=1 pk × log(pk)

The JSS score was also used to evaluate the robustness of the ST-MAGIC atlas, which maps snRNA-seq data to Geo-seq data. Specifically, it assessed the concordance between the reconstructed gene expression in spot from the snRNA-seq and the raw gene expression in spot from the Geo-seq data.

In addition, we computed the silhouette score to evaluate the effectiveness of cell type-specific gene-linked peaks in clustering cells in UMAP visualizations. Using pre-annotated cell types as clustering labels, the silhouette score was calculated based on the 2nd to 10th latent semantic indexing components via the ‘cluster’ package. A higher silhouette score indicates that cells of the same type are more closely grouped.

To systematically evaluate the stability and robustness of BioCRE in predicting CREs across varying sequencing depths, comprehensive down-sampling experiments were conducted. Initially, paired single-cell transcriptomic (snRNA-seq) and single-cell chromatin accessibility (snATAC-seq) datasets were jointly down-sampled to 20, 40, 60, 80, and 100% of their original sequencing depth. This allowed for an assessment of prediction consistency across methods under increasingly sparse data conditions. Subsequently, two complementary experimental paradigms were implemented: first, snATAC-seq data were maintained at full depth while the sequencing depth of snRNA-seq was progressively reduced; second, snRNA-seq data were preserved in their entirety while the depth of snATAC-seq data was incrementally decreased. In all experimental conditions, the stability of predicted regulatory linkages was quantitatively evaluated based on the proportion of shared linkages identified across different sequencing depths.

Performance on public multi-omic datasets with corresponding promoter capture Hi-C data: To systematically compare the predictive performance of BioCRE, we conducted additional evaluations using publicly available data. We analyzed published multi-omics datasets (PBMC: https://www.10xgenomics.com/cn/datasets/pbmc-from-a-healthy-donor-granulocytes-removed-through-cell-sorting-10-k-1-standard-1-0-0 and GM12878: GSE166797), with corresponding promoter capture Hi-C data (PBMC: EGAS00001001911 and GM12878: https://zenodo.org/records/3255048) as reference datasets for validating predicted linkages.

Standardized genomic scope of predicted linkages across all three tools were restricted to ±500kb upstream and downstream of the gene body. Notably, although the total numbers of predicted linkages are theoretically consistent across methods, minor discrepancies may arise due to each tool’s internal preprocessing steps—such as gene and peak filtering—which may exclude certain elements prior to linkage inference. The filtered geneset and peakset shared by all three tools were used for further linkage prediction.

For comparison between the public multi-omic datasets and the reference dataset, any predicted linkages overlapping the reference set was considered as positives, while non-overlapping linkages were treated as negatives. Due to the distinct features for different data modules between the predicted linkages and the reference dataset (PCHiC), the overlap is inherently limited, leading to a sparse set of positive linkages. Here, to ensure balanced and robust evaluation, we randomly sampled negative linkages ten times at each positive-to-negative ratio. Specifically, for the PBMC dataset, we evaluated the AUPRC through ‘sklearn.metrics’ package under varying positive-to-negative ratios: 1:1, 1:2, 1:3, 1:4, and 1:5. For GM12878, we used ratios of 1:1 and 1:2, as the larger number of positive linkages—relative to PBMC—split the accessibility of AUPRC. Besides, we used the p-value for BioCRE, the p-value for Signac, and the FDR for ArchR.

Multi-omics data pre-processing
Raw FASTQ files of snRNA-seq and snATAC-seq were processed by Cell Ranger ARC (v2.0.1) with default mapping parameters. Reads were mapped to the mouse reference genome (mm10) and quantified with gene annotation (GRCm38.98). The quality control of snRNA-seq and snATAC-seq data was performed separately, retaining only high-quality cells that passed both criteria for subsequent analyses. For snRNA-seq data, the following quality control criteria were applied: (1) cells that had fewer than 1000 expressed genes, or over 10% of UMIs derived from mitochondrial genome were removed; (2) cells identified as doublet (DoubletFinder, v2.0.3) were removed; (3) clusters generated by Seurat using default parameters were removed if AY036118 and Gm42418 genes were highly expressed. For snATAC-seq data, cells were required to have a minimum of 3000 fragments, a minimum fraction of fragments in peaks of 0.15, a minimum of TSS enrichment of 9, a maximum of blacklist ratio of 0.05 and a maximum of nucleosomal signal strength of 4.

snRNA-seq processing
After converting the filtered snRNA-seq count matrix to a Seurat object, the standard procedures were performed using the Seurat package. For each stage (E6.5, E6.75, E7.0, E7.25, and E7.5), the following analysis steps were applied: (1) the count matrix was log-normalized for each cell with a scaling factor of 10,000; (2) 2000 variable genes were selected by ‘vst’ method; (3) the ‘nFeature_RNA’ was regressed out during the scaling process. For the combined dataset, the ‘stage’, ‘nFeature_RNA’, and ‘percent.mt’ were regressed out; (4) The scaled variable genes were projected into a low-dimensional space using principal component analysis; (5) The shared nearest neighbor (SNN) graph was built based on the distance of PC space (nPCs=50); (6) The Louvain algorithm was used to identify clusters based on SNN graph (resolution=1); (7) The UMAP embedding was calculated using ‘uwot’ with the top 10 PCs.

snATAC-seq processing
As the precise boundaries of peaks differed across embryonic stages, a common set of peaks (20 bp < length of peak < 10,000 bp) across all stages was generated using the ‘reduce’ function in GenomicRanges (v1.48.0) package. The fragment file was then remapped to these common peaks to get the chromatin accessibility matrix. Then, the remapped chromatin accessibility matrix and correspondent fragment file were processed using the Signac package. The matrix was normalized by Term Frequency Inverse Document Frequency, and all peaks were used for dimension reduction using Singular Value Decomposition. The UMAP embedding was computed utilizing dimensions from 2 to 10, with min.dist set to 0.2.

Clustering and cell type annotation
The annotation of multi-omics is based on snRNA-seq data, while snATAC-seq annotation is facilitated via the consistent cellular barcodes between the transcriptome and chromatin accessibility data, ensuring harmonized identification across datasets. Given that cell type diversification during gastrulation, the cell types of snRNA-seq data were annotated in a stage-by-stage manner. Initially, the snRNA-seq data is clustered through the standard Seurat pipeline (vars.to.regress=‘nFeature_RNA’, k.param=20, dims=1:10, resolution=1) with a preliminary annotation to the three primary groups: Embryonic, ExE endoderm, and ExE ectoderm, then cells within the Embryonic group were subjected to an iterative clustering process to achieve more finely detailed cell types. The iterative clustering also employs the Seurat pipeline (vars.to.regress=c(‘nFeature_RNA’, ‘percent.mt’), k.param=30, dims=1:30, resolution=0.9, these parameters remain constant throughout the iterative process) to conduct a clustering process independently on each cluster, with the stipulation that clusters cannot be further subdivided or any cluster containing fewer than 50 cells will not be subjected to additional clustering. Following several iterations, the definitive clusters are established. The integration of cluster-specific genes (defined by the COSG algorithm63) and marker genes from the published gastrulation transcriptome atlas3,6 were mainly used to annotate the clusters. Ambiguous clusters were clarified by correlating them with the closest matching cell type found in existing public datasets (GSE87038). The integration of the public dataset with our snRNA-seq data is facilitated through the use of Canonical Correlation Analysis (CCA) within the Seurat package. The Pearson correlation between the pseudo-bulk transcriptomic profiles of our iterative clusters and those of cell types in the public dataset is employed to refine annotations.

The similarity between snRNA-seq data and published datasets
To verify the reliability of the annotation, we compared it with the public dataset (GSE87038) at three levels: the proportion of cell types, the correlation of marker genes, and the correlation of the whole transcriptome. First, the shared stages and corresponding cell type data were extracted from the public dataset. Then, the CCA method in Seurat was used to integrate these two datasets. The proportion of each cell type and the expression of marker genes (Sox2, Mesp1, Sox17, and Hnf4a) were statistically compared between the two datasets, along with the pseudo-bulk expression of 2000 high-variance genes across cell types.

Cellular trajectory analysis with TOME
The cellular trajectory spanning from E3.5 to E7.5 was reconstructed based on scRNA/snRNA-seq data using the TOME. The processed dataset covering the developmental stages from E3.5 to E6.25 was made available by the TOME. For the subsequent stages from E6.5 to E7.5, we employed our snRNA-seq data. The identification of the nearest ancestral node for each subsequent stage node is repeated for 500 times based on 50 nearest neighbors. Only those edges with a probability greater than 0.2 will be preserved. Given that the annotated cell types do not align perfectly with those in the pre-stage dataset, we have designated “Embryonic visceral endoderm” as “visceral endoderm,” and “Extraembryonic visceral endoderm” as “ExE endoderm” within our dataset for consistency.

WNN embedding for embryonic cells
The RNA and ATAC modality from embryonic cells, mapping to Geo-seq data, across five development stages, were integrated to construct a weighted nearest neighbor graph by the ‘FindMultiModalNeighbors’ method in the Seurat package. In the snRNA-seq data process, the ‘stage’, ‘nFeature_RNA’, and ‘percent.mt’ were regressed out during normalizing. In snATAC-seq data process, the ‘min.cutoff’ was set to ‘q0’ in finding top features.

Inference of the TF regulatory network with SCENIC+
The pre-processed multi-omic RNA and ATAC matrices, encompassing the full dataset and axial mesendoderm lineages, are utilized as inputs for SCENIC+ (v 0.1.dev452), facilitating the construction of the TF regulatory network. SCENIC+ was executed with the majority of its settings at their default values, with adjustments made only to parameters specific to the dataset. For the full dataset, the peak matrix is regenerated utilizing stage-specific cell type pseudobulk data through MACS2 (v2.2.7.1). During the cisTopic (v1.0.2) module, an optimized set of 16 topics is selected for further analysis. Peaks within a 2kb~500kb range upstream or downstream of genes are incorporated into the SCENIC+ modeling process. eRegulons featuring positive region-to-gene associations, along with a correlation exceeding 0.15 between the TF cistrome and region-based activity, were retained, resulting in a set of 302 high-quality eRegulons. For axial mesendoderm lineages, stage-specific peaks for each cell type are also generated. In cisTopic, the optimized number of topics is set at 16. The range for peaks is consistent with that of the full datasets. eRegulons are retained based on a correlation threshold of 0.15 between the TF cistrome and region-based activity, yielding a total of 118 high-quality eRegulons.

Differential ChIP-seq peaks of the NODAL and WNT signaling
For each ChIP-seq dataset, raw sequencing data was downloaded from the relative resource (NODAL signaling: GSE70486; WNT signaling: GSE43565, GSE162774) and then mapped to the mouse genome (mm10). Peak calling was performed against the corresponding input sample using MACS with the parameters “--shiftsize=100 --nomodel --keep-dup=all”. Subsequently, we merged all resulting peaks of samples into a consensus list of genomic regions and counted the reads within those regions using MAnorm2_utils (v1.0.0) with the parameters “--min-peak-gap=150 --typical-bin-size=2000 --shiftsize=100 --filter=blacklist”. The blacklist regions of mm10 were obtained from Amemiya et al.64. Then, MAnorm2 (v 1.2.2) were used to perform differential peak analysis between samples with default parameters (refer to https://cran.r-project.org/web/packages/MAnorm2/vignettes/MAnorm2_vignette.html). The differential peaks were identified using the cutoffs adjusted p-value < 0.01 and fold change > 2.

Axial mesendoderm lineage cells
The selection of axial mesendoderm lineage cells was guided by trajectory analysis depicted in Fig. 1c. This lineage originates from the primitive streak cells at E6.5, and subsequently progresses through the Anterior primitive streak at E6.75 and E7.0. It then developed into Node precursors and Anterior mesendoderm precursors at E7.25. The lineage culminates at E7.5 with Node and Anterior mesendoderm. To ensure more reliable lineage-related cells, we retain only those cells from ancestral stages that are found among the nearest neighbors as calculated by TOME, except for cells in the E6.5 Primitive streak.

Driver genes of axial mesendoderm lineage with CellRank
The sorted BAM files derived from CellRanger were processed using Velocyto (v0.17.17) to compute RNA velocities for individual cells. The mask GTF file was obtained from the UCSC database, and the genome annotation file originates from CellRanger. This procedure enables detailed analysis of transcriptional dynamics at the single-cell level by leveraging the spliced and unspliced transcript information contained within the BAM files. The velocity matrix specific to the axial mesendoderm lineage cells has been extracted from the loom file utilizing the SeuratWrappers (v 0.3.0) R package. Subsequently, the velocity matrix was employed for velocity estimation through the application of the ‘dynamical’ model in scVelo (v0.2.5). The cell-cell transition matrix was calculated employing the ‘cellrank.tl.transition_matrix’ function in CellRank (v 1.5.1). Terminal states were set as Node and Anterior mesendoderm manually. Following this, the GPCCA estimator was utilized to compute both the fate probabilities and driver genes. The top 1000 smoothed driver genes, which exhibited the highest correlation with pseudotime, were selected as the key drivers. The driver genes identified in the BioCRE linkage analysis of the axial mesendoderm lineage and their linked peaks were selected to plot the heatmap. The pseudotime for the lineage was computed utilizing the ‘scanpy.tl.dpt’ function (scanpy, v 1.9.1), incorporating the eRegulon RNA and ATAC AUC matrices obtained from the SCENIC+.

KO-omics data processing
The preprocessing for KO-omics data (WT, Noto KO, Pou6f1 KO embryos at the E7.5 stage) is consistent with those previously mentioned. The annotation method for WT snRNA-seq data remains consistent with the previous approach, whereas for Noto KO and Pou6f1 KO samples, snRNA-seq annotation is performed using the ‘MapQuery’ in Seurat, with annotated WT snRNA-seq data serving as a reference. The annotation for snATAC-seq data was aligned with that of snRNA-seq by leveraging shared cell barcodes.

Web server implementation
The 10x scMultiome data and reconstructed Geo-seq data are available through the ST-MAGIC web server, which was implemented using Shinyapps from R. This web service provides three main functionalities: (1) Transcriptomics in snRNA-seq data, which illustrates gene expression on UMAP and cell type grouping; (2) Epigenomics in snATAC-seq data, which shows peak accessibility, linkage, and motif information; (3) Spatial omics of multiome data, which includes gene, peak, cell type, gene set, and peak set enrichment visualized on a corn plot. ST-MAGIC enables users to explore spatial and temporal patterns in a user-friendly manner. For transcriptomics, profiled genes can be searched for their expression patterns over space and time. For epigenomics, profiled peaks and motifs can be used to investigate spatial and temporal accessibility patterns. In addition, ST-MAGIC provides a selection list of gene sets or peak sets from SCENIC+ to TF enrichment analysis.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In “Spatiotemporal Regulatory Logics of Mouse Gastrulation”, Yang et al describe two resources, ST-MAGIC and ST-MAGIC+, that project new single cell multimodal transcriptome and ATAC-seq data collected for 4 timepoints between E6.5-E7.5 during mouse gastrulation onto the published Geo-seq spatial profiling atlas. To identify gene regulatory networks, the authors introduce a new software package, BioCRE, to analyze multimodal transcriptome information with ATAC-seq data to link genes with accessibility peaks. The authors describe several vignettes including dissection of the GRNs and enhancers that guide mesendoderm development to showcase how these resources can be used. Through those vignettes, they identify novel enhancer elements and transcription factors that may affect differentiation.\n\nThe manuscript is a tour-de-force work integrating many published data sets and tools. Indeed, most figures both in the main text and supplementary had a very large number of panels. The manuscript was reasonably well written so that I could follow the logic through analysis of this complex data set and I commend the authors for providing many data panels. The schematics as the first panel for several figures was helpful. Moreover, I appreciated the generation of new data to follow up new hypotheses (Pou6f1 KO and Noto KO multi-seq data) and the RNA-scope experiments for validation.\n\nNonetheless, it was difficult to assess many of the claims due to a lack of quantification, statistical tests, and useful visualizations. While the development of new analysis platforms like BioCRE is useful, BioCRE was not adequately benchmarked against other methods and so it is difficult to evaluate its contribution to this manuscript.\n\nGeneral comments:\nI recommend a different title for this manuscript as the current one suggests more comprehensive profiling across time.\nThe legend captions should be more descriptive and explain nomenclature used in the figure. For example, I think peaks are often shortened to “p” and genes are often shortened to “g”.\n\nPoint-by-point comments sectioned by the main figure referenced:\nFigure 1\nED Fig 1a: The quality of the E6.75 data set looks poor compared to the other data sets (notably fewer genes and UMIs detected, more mitochondrial RNA). While I believe the majority of the vignettes presented in the manuscript rely on the later timepoints, the utility of this data set as a resource would require that the data is of high quality. Is it possible that the lower quality is reflective of biology or was there noticeable difference in sample quality or preparation for the E6.75 data set?\nLines 114-116: Please include the number of embryos pooled for each time point.\nED Fig 1g: What method was used to count approx. cell numbers of embryos? A table highlighting the number of cells per embryo per timepoint would be useful to give a sense of embryo staging.\nED Fig 3a: Please add details to the legend. How are the values computed? Is it the cluster compared against the rest of the cells? Does the colorbar refer to log2 fold change?\nLines 154-165: There is not enough analysis benchmarking BioCRE against Signac and ArchR. What other metrics can be used to measure the consistency of gene expression and peaks? What are the set characteristics between the linkage predictions of the three tools (does BioCRE identify a super set of linkage? Are many predictions non-overlapping?). What are the run-time characteristics for BioCRE compared to the other tools? How does sampling affect BioCRE performance? Example: downsample the RNA-seq and ATAC-seq together, down-sample one while holding the other constant, etc. How does the median gene-to-peak distance compare between methods?\nFigure 2\nED Fig 4b: There is a general trend of increasing JSD score with timepoint. Can the authors speculate why? Is it due to the number of cell types?\nED Fig 4d: The claim is that the expression of Sox2 closely matches between the GEO-seq data and ST-MAGIC gene expression. This is hard to evaluate from the figure panel. A scatter plot would aid evaluation of the claim. Moreover, the proposed visualization can be done for all genes and provide a better measurement to understand the similarity between GEO-seq and ST-MAGIC.\nFig 2g and h: The schematic is not useful.\nLines 272-274: Are the differences statistically significant?\nED Fig 6g: How many gastruloids where generated and measured in this experiment? Did the authors observe morphological differences between the cell lines or different generation efficiencies?\nFigure 3\nLines 307-308: How were Etv2, Gata2 and Tal1 selected from the SCENIC+ candidates? Do the other TF candidates that fit the TF expression hierarchy display the same paradigm of the upstream factor functioning as a priming factor responsible for opening chromatin for use by later factors?\nLines 342-346: How many target genes and target peaks are there for NOTO?\nFig 3e: It is difficult to evaluate the claim that target genes of NOTO and NOTO expression have matched spatiotemporal distribution with the provided visualization. Moreover, the authors should leverage their comprehensive integrated atlas, to evaluate how many genes and their targets nominated by SCENIC+ display similar behavior. It is not clear why Noto is selected from the candidates SCENIC+ nominates for mesendoderm development. It does not appear to be highlighted in ED Fig 7.\nFigure 4\nLines 374-380: There are claims about the specific and global spatial expression patterns for different factors. While these assessments can be made by eye using the overlay onto the spatial coordinates, a quantitative assessment is more valuable. For example, if the authors plot a histogram of the values for each spot, the factors that show specific expression patterns and those that are global will be obvious. To be clear, I think both the plots the authors provide and the quantifications I recommend are needed.\nFig 4b,c: The claim is that Wnt3 and Dkk1 have overlapping distributions. This is difficult to assess by the provided panels. Moreover, can any statistics be applied to this analysis? Are the peaks enriched?\nFigure 5\nLines 463-466: There is not a figure panel that shows the co-expression pattern for Shh, Noto, and Foxj1. I believe Fig 5e might be the panel but it does not have Shh or Noto highlighted. It is difficult to assess the claim that these two distinct subtypes exist.\nLine 514: What does a “highly specific distribution” mean for POU6F1? Can the authors support this claim with a statistic?\nFigure 6\nLine 531: “The target genes start to show abundant enrichment at the same region”. From inspection of Fig 6a, it appears to me that the POU6F1 target gene set expression appears before the expression of Pou6f1 with the former appearing at E7.25. Can the authors explain the differential timing of these events?\nLines 532-535: The authors suggest that “early and intermediate stage TFs may function as priming factors by opening a broad range of chromatin regions”. This is a reasonable hypothetical model. Can the authors support or refute this claim with additional analysis? For example, are the target peaks from early and intermediate stage TFs enriched for terminal TF targets? From Fig 5m,n, the number of peaks and genes targeted by Eomes (42p, 35g respectively) is far less than the number of peaks and genes targeted by Pou6f1 (127p, 99g respectively). It's not clear to me how to incorporate these numbers into the proposed model.\nLines 547-551, Fig 6c: The authors described that EOMES and FOXA2 ChIP-seq data shows enrichment of these TFs around the peaks of NOTO and POU6F1 targets. First, the text does not describe the ChIP-seq data. The citations are in the legend and only by looking at the papers could I discern what the ChIP-seq data is, which is from EBs and in vitro differentiation, respectively I believe. Second, there is no statistical test for the enrichment but I agree that by eye, there is a peak for both Eomes and FoxA2. However, it does not appear from the data presented in Fig 6c that the Eomes KO has an effect on these peaks implying that the chromatin regulation may be controlled by other factors. I would agree that F2-KO does have an effect but does not entirely explain the regulation at those target loci. Will the authors please explain or modify the following statement in the text: “EOMES and FOXA2 are likely to be responsible for pre-opening… target chromatin regions”.\nMinor points:\nFig 1e: I believe JSD index is incorrect. JSD measures the similarity between two distributions. The Jensen-Shannon diversity index is a special case of JSD.\nED Fig 1g: needs a color legend for cell types\nED Fig 2b: The text describes cell types that are enriched at earlier or later time points. To aid this visual assessment, it would be helpful if the x-axis had the same scale for each time point.\nED Fig 3h: units on y axis. How many replicates per condition?\nED Fig 10e: Needs a summary plot at the top of the heatmap such as in Fig 2i.\nLine 486: It is not clear what Combinatorial pseudo-time based tracking means.\nED Fig 12f: Caption missing explanation of acronyms in panel (EB – embryoid body, MEN - ?)\nLine 569, ED Fig 13h: I realize the numbers are small, but it would nonetheless be useful for the authors to derive a likelihood (p-value) for the axis turning randomization.", "reviewer_2": "In this manuscript, Yang and colleagues present a spatiotemporal dataset capturing mouse gastrulation, generated using 10x Genomics Multiome data. The study involves a significant amount of data integration and projection, incorporating external datasets such as Geo-seq and ChIP-seq. To enable gene regulatory network (GRN) analysis, the authors developed a new computational tool called BioCREs, which they claim performs better than existing pipelines. They also adapted the Tangram algorithm—originally designed for high-resolution spatial transcriptomics—for use with their data to align it with Geo-seq references. Through this combined approach, they aim to map spatial and temporal regulatory changes during gastrulation.\n\nThe scope of this work is impressive, and the integration of multiple data types is interesting. The authors also acknowledge the limitations of current spatial resolution and mismatched cell type references, which are important caveats. I appreciate the validation experiments included, which lend support to some of their conclusions. That said, I do have several concerns, especially regarding the BioCREs pipeline and the modified Tangram tool.\n\nMajor Concerns:\n1. The claim that BioCREs outperforms Signac and ArchR is based on a single performance index. This is insufficient to justify the claim. Benchmarking against established regulatory element datasets (e.g., ENCODE) is necessary to validate accuracy and specificity.\n2. Distal vs. Proximal regulatory region comparison: The comparison between \"Distal High\" and \"Proximal\" groups is biased by excluding the \"Distal Low\" group. A more appropriate comparison would be between all distal regions and proximal ones, with proper statistical testing (e.g., p-values) to support conclusions.\n3. Insufficient benchmarking for BioCREs and modified Tangram: The manuscript lacks technical detail and performance benchmarking for both pipelines. More rigorous testing is needed to demonstrate their robustness and reliability.\n4. Incorporating data from high-resolution platforms such as Visium HD or Stereo-seq would significantly strengthen the spatial claims. Even one matching experiment would add credibility.\n\nMinor Comments:\n• The manuscript is dense, especially in later sections discussing biological processes. Streamlining these parts would improve accessibility.\n• Some extended data figures are poorly organized and should be revised for clarity.\n• Clarity in Extended Data Fig. 2a: The plot is overcrowded due to overlapping cells, making it difficult to interpret. A side-by-side comparison with reference data in separate panels would improve readability.\n• Missing data in Extended Fig. 4: The claim that ST-MAGIC has higher JSD values than Geo-seq is not supported by the figure. This should be clarified or removed."}, "author_response": "Responses to Reviewers\nSummary of major new data\nWe would like to thank the reviewers for their valuable input. We provide a revised version\nof the manuscript with a substantial body of new analyses and experiment data that further\nconfirm and strengthen the conclusions. Please note that we used Fig. 1, 2, 3, etc. to refer\nto figures in the revised manuscript and Response Fig. R1, R2, R3, etc. to refer to figures\nin this letter.\n1. We systematically benchmarked the performance of BioCRE through evaluating the\nbiological relevance, cell type specificity of gene-to-peak linkages, the computational\nefficiency, and consistency of linkage composition under various sample sizes. Moreover,\nthrough using publicly available single-cell multiome and matched promoter-capture HiC\ndatasets, we identified the reliability and accuracy of linkage prediction. These analyses\ncollectively demonstrated that BioCRE outperforms Signac and BioCRE in predicting\nreliable and accurate gene-to-peak linkages.\n2. We compared the ST-MAGIC atlas with a high-resolution Stereo-seq dataset (SEU-3D)\nfor E7.5 embryo. As demonstrated by analyzing cell type-specific transcriptomic\narchitectures, signature gene set distributions, spatial locations of cell types and gene\nexpression patterns, ST-MAGIC shows high spatial consistency with SEU-3D. Furthermore,\nin-depth evaluation highlighted the effective attributes of ST-MAGIC in detecting low-\nabundance genes.\n3. We incorporated additional metrics (such as Pearson correlation, structural similarity\nindex, P-value) and versatile plots (histogram, scatter plot, density plot) to quantitatively\nevaluate the global similarities between ST-MAGIC and Geo-seq atlas, and\ncomprehensively assessed the spatiotemporal relationships for distinct signal cascades\nand eRegulons. Overall, the original spatial claims have been substantially strengthened.\n4. We conducted additional gastruloid differentiation experiments to determine the roles of\nNeo_LRE element and collected more mutant embryos to summarize the phenotypic\nsignificances for Noto mutant and Pou6f1 mutant embryos. Consequently, we\ndemonstrated the roles of Neo_LRE in regulating Lefty2 expression, and confirmed the\naxis turning defects of Noto mutant and Pou6f1 mutant embryos.\nPoint-by-point response letter\nReviewer #1 (Remarks to the Author):\nIn “Spatiotemporal Regulatory Logics of Mouse Gastrulation”, Yang et al describe two\nresources, ST-MAGIC and ST-MAGIC+, that project new single cell multimodal transcriptome\nand ATAC-seq data collected for 4 timepoints between E6.5-E7.5 during mouse gastrulation\nonto the published Geo-seq spatial profiling atlas. To identify gene regulatory networks, the\nauthors introduce a new software package, BioCRE, to analyze multimodal transcriptome\ninformation with ATAC-seq data to link genes with accessibility peaks. The authors describe\nseveral vignettes including dissection of the GRNs and enhancers that guide mesendoderm\ndevelopment to showcase how these resources can be used. Through those vignettes, they\nidentify novel enhancer elements and transcription factors that may affect differentiation.\nThe manuscript is a tour-de-force work integrating many published data sets and tools. Indeed,\nmost figures both in the main text and supplementary had a very large number of panels. The\nmanuscript was reasonably well written so that I could follow the logic through analysis of this\ncomplex data set and I commend the authors for providing many data panels. The schematics\nas the first panel for several figures was helpful. Moreover, I appreciated the generation of new\ndata to follow up new hypotheses (Pou6f1 KO and Noto KO multi-seq data) and the RNA-scope\nexperiments for validation.\nNonetheless, it was difficult to assess many of the claims due to a lack of quantification,\nstatistical tests, and useful visualizations. While the development of new analysis platforms\nlike BioCRE is useful, BioCRE was not adequately benchmarked against other methods and\nso it is difficult to evaluate its contribution to this manuscript.\nResponse: We highly appreciate the reviewer’s positive assessment and constructive\nsuggestions on our manuscript. In this revision, we have carefully addressed the reviewer’s\ninsightful concerns, and the original conclusions have been substantially strengthened.\nGeneral comments:\nI recommend a different title for this manuscript as the current one suggests more\ncomprehensive profiling across time.\nResponse: Thanks for the suggestion. We have revised the title of the manuscript into\n‘Integrated multi-omic atlas reveals the hierarchy of spatiotemporal regulatory networks of\nmouse gastrulation’, which would more accurately capture the central theme of our work,\nnamely, the systematic reconstruction of multi-omic data to elucidate the hierarchical\norganization of regulatory networks in space and time during gastrulation.\nThe legend captions should be more descriptive and explain nomenclature used in the figure.\nFor example, I think peaks are often shortened to “p” and genes are often shortened to “g”.\nResponse: Thanks for the reminder. In this revision, we have carefully revised the figure legend\nto provide more detailed descriptions and clarified the nomenclature in the figures. As the\nreviewer correctly mentioned, ‘p’ is short for ‘peaks’ and ‘g’ is short for ‘genes’.\nPoint-by-point comments sectioned by the main figure referenced:\nFigure 1\nED Fig 1a: The quality of the E6.75 data set looks poor compared to the other data sets (notably\nfewer genes and UMIs detected, more mitochondrial RNA). While I believe the majority of the\nvignettes presented in the manuscript rely on the later timepoints, the utility of this data set as\na resource would require that the data is of high quality. Is it possible that the lower quality is\nreflective of biology or was there noticeable difference in sample quality or preparation for the\nE6.75 data set?\nResponse: Thanks for your insightful reminder. We would like to clarify that the E6.75 samples\nwere processed using the same experimental pipeline and quality control procedures as all\nother developmental stages. In this revision, to further assess the biological reliability of the\nE6.75 dataset, we have conducted several validation analyses as follows.\nFirst, we determined the transcriptome similarities between our E6.75 dataset and the\npublished dataset1 (Response Fig. R1a). As expected, we found that the corresponding cell\ntype from both datasets is well correlated.\nNext, we performed Pearson correlation analyses of the E6.75 transcriptome with that of\nneighboring stages (E6.5, E7.0). As shown, all relevant cell types exhibit high correlation with\ntheir counterparts from the other two stages (Response Fig. R1b).\nFinally, we examined the association of the cell type-specific genes and cell type-specific\naccessible peaks for each cell type at E6.75. As shown in Response Fig. R1c, expression of\ncell type-specific marker genes is coupled with enrichment of accessible peaks. Taking the\nhaematoendothelial progenitors as an example, we found that they show both specific\nexpression of haematoendothelial marker-Fli13 (Response Fig. R1d) and specific enrichment\nof accessible peaks around the Fli1 locus (Response Fig. R1e).\nTogether, although the E6.75 dataset exhibits slightly lower sequencing quality metrics, it\nremains biologically informative and robust for the analyses performed in this study.\nResponse Figure R1. Molecular characterization of E6.75 multi-omic data. a. Heatmap showing the Pearson\ncorrelation between E6.75 dataset in this study and the published E6.75 dataset.1 Color legend represents the\ncorrelation. b. Heatmap showing the Pearson correlation for cell counterparts between neighboring stages (comparison\nbetween E6.5 and E6.75 dataset at the left, comparison between E6.75 and E7.0 at the right). Color legends represent\nthe correlation. c. Heatmap showing the differentially expressed genes (left) and differential accessible peaks (right)\nacross all E6.75 cell types. Color legends represent the scaled level of gene expression (left) and chromatin\naccessibility (right). d. UMAP plot showing the expression of Fli1 in haematoendothelial progenitors. e. Coverage plot\nshowing the distribution of accessible peaks around the Fli1 locus. Accessible peaks were highlighted in red.\nLines 114-116: Please include the number of embryos pooled for each time point.\nResponse: Thanks for the kind suggestion. We have included the exact number of embryos\npooled for each time point in the revised manuscript (Line 116-117). Specifically, a total of 12\nembryos for E6.5, 10 embryos for E6.75, 6 embryos for E7.0, 6 embryos for E7.25, and 5\nembryos for E7.5.\nED Fig 1g: What method was used to count approx. cell numbers of embryos? A table\nhighlighting the number of cells per embryo per timepoint would be useful to give a sense of\nembryo staging.\nResponse: Thanks for the insightful reminder. The approximate cell numbers and\ncorresponding developmental stages for each embryo were determined based on a\ncombination of well-established literatures1,2 (as shown in Response Fig. R2a, b) and\nmorphological assessment from our team’s extensive experience in early mouse\nembryogenesis. The representative images for embryo per timepoint has been shown\n(Response Fig. R2c; Supplementary Fig. 1a).\nAs suggested, to deliver a clear sense of cell numbers, we have revised and attached a table\nhighlighting the numbers of cells per timepoint in the revised manuscript (Response Fig. R2d;\nSupplementary Fig. 1e).\nResponse Figure R2. Estimated cell numbers for the mouse gastrula. a. Table summarizing the cell numbers of\nthe mouse embryo. Table was from a published study.2 b. Bar plot summarizing the estimated cell numbers of the\nmouse embryo. Figure was from a published study.1 Estimated cell numbers for correspondent embryonic stages\ncovered in our study were highlighted with red triangles. c. Representative images of embryos collected at the time\npoints indicated. d. Detailed information for estimated cell numbers and multi-omic profiled cell numbers in this study.\nED Fig 3a: Please add details to the legend. How are the values computed? Is it the cluster\ncompared against the rest of the cells? Does the colorbar refer to log2 fold change?\nResponse: Thank you for these comments. We apologize for the lack of clarity regarding the\ndetails of this figure. To clarify, the values are scaled transformations of the original raw data\nand do not represent cluster-versus-rest comparisons or log2 fold changes. Specifically, for the\ngene expression heatmap, the values are scaled results derived from gene count data, adjusted\nfirst for library size and then subjected to log transformation. In the peak chromatin accessibility\nheatmap, the values are scaled results derived from peak count data, corrected for library size\nand further processed using TF-IDF (Term Frequency-Inverse Document Frequency)\nnormalization. Lastly, the motif enrichment heatmap displays scaled motif activity scores, which\nwere computed from the peak count matrix using the chromVAR method to quantify the\nenrichment of specific transcription factor motifs within open chromatin regions.\nThe details of this figure have been included in the revised figure legend of Supplementary\nFig. 3a-c.\nLines 154-165: There is not enough analysis benchmarking BioCRE against Signac and ArchR.\nWhat other metrics can be used to measure the consistency of gene expression and peaks?\nWhat are the set characteristics between the linkage predictions of the three tools (does\nBioCRE identify a super set of linkage? Are many predictions non-overlapping?). What are the\nrun-time characteristics for BioCRE compared to the other tools? How does sampling affect\nBioCRE performance? Example: downsample the RNA-seq and ATAC-seq together, down-\nsample one while holding the other constant, etc. How does the median gene-to-peak distance\ncompare between methods?\nResponse: Thanks for this valuable suggestion. We fully agree that a more comprehensive\nbenchmarking strategy is essential to substantiate our claim that BioCRE outperforms existing\ntools.\nIn this revision, to benchmark the performance of BioCRE, we have conducted a series of\nsystematic evaluations to rigorously compare the performance of BioCRE against the state-of-\nthe-art methods: Signac and ArchR.\n1. The characteristics of predicted gene-to-peak linkages:\nUsing our multi-omic dataset from the mouse gastrula, we evaluated the number, biological\nrelevance, and cell type specificity of predicted gene-to-peak linkages. As shown, BioCRE\nidentified 44,644 gene-to-peak linkages, with a median gene-to-peak distance of 127,175 bp\n(Response Fig. R3a, b). Comparative analyses revealed that 28,472 linkages (64%) were\nspecifically identified in BioCRE (Response Fig. R3c, d). Gene ontology analyses of genes\nassociated with these BioCRE-specific linkages (28,472 linkages) revealed significant\nenrichment in embryo development, including heart development, liver development, and\nnervous system development (Response Fig. R3e). This result suggests that these BioCRE-\nspecific linkages could be responsible for the regulation of embryo development related genes.\nTo assess the cell type specificity of predicted linkages, we analyzed the accessibility of peaks\nlinked to cell type-specific genes identified by each tool (Response Fig. R3f, g). Using UMAP\nvisualization and quantifying clustering coherence with the Silhouette score, we found that\nlinkages predicted by BioCRE gets a higher score compared to those from Signac and ArchR,\nindicating a more effective capability for BioCRE to capture cell type-specific linkages\n(Response Fig. R3f, g).\nTogether, these results demonstrate that BioCRE could be adept at identify cell type-specific\nlinkages, which are likely to play key roles in orchestrating embryo development.\n2. The computational efficiency and robustness for predicting linkages:\nWe assessed the computational efficiency and robustness of different methods under varying\nsampling conditions. Overall, we found that the built-in deep-learning bi-orientation regression\nmodel for BioCRE does not incur excessive run-time costs, showing comparable performance\nto ArchR (Response Fig. R3h). To evaluate the robustness of linkage identification, we\nchecked the consistency of predicted linkages across three distinct down-sampling scenarios:\n In the first scenario, where both RNA and peak reads were down-sampled concurrently,\nBioCRE demonstrated the highest recovery rate of linkage identification (Response Fig.\nR3i). To specify, 40.9% linkages (31,911 out of 78,109) are constantly identified in BioCRE,\nwhile only 8.2% linkages (3,677 out of 44,705) for Signac and 5.7% linkages (4,479 out of\n79,067) for ArchR show consistent results among variable sample sizes (Response Fig.\nR3i-1).\n In the second scenario, where only RNA reads were down-sampled while peak reads\nremain unchanged, we found that 44.6% linkages (21,406 out of 48,043) are consistently\ncaptured for BioCRE (Response Fig. R3i-2). In the meantime, 42.9% linkages (18,335\nout of 42,937) for Signac and 61.9% linkages (38,510 out of 62,226) for ArchR are stably\nidentified (Response Fig. R3i-2).\n In the last scenario, we specifically down-sampled peak reads and then performed linkage\nidentification (Response Fig. R3i-3). Consequently, we found that 52.1% linkages\n(49,476/95,042) for BioCRE are captured across all down-sampling groups (Response\nFig. R3i-3). In contrast, only 24.1% linkages (9,122/37,894) for Signac and 46.2% linkages\n(28,388/61,422) for ArchR show the same pattern (Response Fig. R3i-3).\nCollectively, these analyses demonstrate that BioCRE achieves a favorable balance between\ncomputational efficiency and robustness, exhibiting consistency and high reproducibility in\nlinkage prediction across diverse and reduced sequencing depths—particularly under\nconditions of limited chromatin accessibility data (Scenario 1 and 3)—highlighting its reliability\nin practical, resource-variable settings.\n3. The reliability and accuracy of predicted linkages:\nGiven the absence of a genome-wide gold standard for gene-to-peak linkages in the mouse\ngastrula, we turned to external single-cell multiome datasets with matched, high-confidence\npromoter-enhancer interaction data to benchmark performance. Specifically, we applied\nBioCRE, Signac and ArchR to publicly available single-cell multiome data from human PBMC\ncells and Gm12878 cells, using promoter-capture Hi-C (PCHiC) data as an orthogonal\nreference for true regulatory interactions (Response Fig. R3j).\nTo enable a comprehensive and balanced evaluation, we defined both positive and negative\ncontrols: positive linkages were inferred from PCHiC-confirmed interactions, while negative\nlinkages were established as peak-gene pairs showing no chromatin interactions within a 1 Mb\nrange in the PCHiC data. We evaluated the performance using multiple following metrics:\n The metrics of Macro-F1 score across all linkages showed that BioCRE consistently\nachieved higher scores than both Signac and ArchR across all tested cell types, with strong\nreproducibility across replicates (Response Fig. R3j-1). This indicates that BioCRE’s\npredictions more closely align with the experimentally determined chromatin architecture.\n The Macro-F1 score for linkages of all three tools shared genes revealed that BioCRE\noutperformed both Signac and ArchR in most cell types, where its advantage was less\npronounced (Response Fig. R3j-2).\n The True Positive Rate (TPR) for linkages unique to each method was also assessed.\nBioCRE-specific linkages show comparable or higher level of TPR than those of Signac\nand ArchR (Response Fig. R3j-3), suggesting that BioCRE-specific predictions are not\nonly robust but also biologically credible.\nThus, through incorporating external dataset and gold standard chromatin interaction dataset,\nwe demonstrated that BioCRE exhibits consistent superiority in linkage identification than\nSignac and ArchR.\nThese results have been included in the revised manuscript (Line 162-176; Supplementary\nFig. 4), the detailed method has been provided in Line 982-1038, and we also acknowledged\nthe significance of systematically benchmark of BioCRE performance and discussed about the\ncurrent limitations in Line 623-632, Line 703-705.\nResponse Figure R3. Benchmark of BioCRE. a. Bar plot showing the number of gene-peak linkages identified by\nBioCRE, Signac, and ArchR, respectively. b. Histogram showing the gene-peak distance for linkages identified by\nBioCRE, Signac, and ArchR, respectively. c. Venn diagram showing the intersection of predicted gene-peak linkages.\nd. Bar plot highlighting the ratio of non-overlapping linkage ratios for each tool. e. Gene ontology analyses for BioCRE-\nspecific captured linkages associated genes. f. UMAP plot based on the chromatin accessibility of cell type-specific\ngenes associated peaks. g. Silhouette score evaluating the cell type distribution specificity. h. Line plot depicting the\nruntime for each tool with variable proportion of reads. BioCRE was running on a single GPU (A100), while ArchR and\nSignac were executed using 20 CPU cores for parallel computing, respectively. i. The effects of sequencing reads\ndown-sampling on linkage composition consistency for each tool. The constant linkage numbers for each scenario were\nhighlighted in red. j. The prediction and validation of gene-peak linkages using external multiome dataset and gold\nstandard promoter-enhancer interaction dataset.\nFigure 2\nED Fig 4b: There is a general trend of increasing JSD score with timepoint. Can the authors\nspeculate why? Is it due to the number of cell types?\nResponse: Thanks for the insightful observation. The increasing trend in JSD scores (revised\nas JSS score, Line 162-164) across developmental timepoints (from E6.5 to E7.5) reflects key\nbiological dynamics during mouse gastrulation. As the reviewer correctly pointed out, the\nnumber of cell types is increasing during mouse gastrulation, which coincides with the gradual\ncell fate specification of pluripotent epiblast cells from E6.5 to E7.5 (Response Fig. R4a;\nSupplementary Fig. 1h).\nConcurrently, the embryo develops pronounced spatial organization, with cells adopting\nposition-specific identities along anterior-posterior, proximal-distal, left-right axes (Response\nFig. R4b; Supplementary Fig. 7a). This emerging regionalization leads to sharper, more\nspatially restricted gene expression patterns. As JSS scores measure the specificity and\nconfidence of spatial gene-to-spot associations, they naturally increase as gene expression\nbecomes more localized and cell type composition more structured.\nResponse Figure R4. Cell type diversification during mouse gastrulation. a. Dynamics of cell type frequency per\nembryonic stage. b. Regionalized cell type distribution across the five embryonic stages as visualized by ST-MAGIC.\nED Fig 4d: The claim is that the expression of Sox2 closely matches between the GEO-seq\ndata and ST-MAGIC gene expression. This is hard to evaluate from the figure panel. A scatter\nplot would aid evaluation of the claim. Moreover, the proposed visualization can be done for all\ngenes and provide a better measurement to understand the similarity between GEO-seq and\nST-MAGIC.\nResponse: Thanks for your insightful reminder. As suggested, we have added scatter plots to\nassess the spot-by-spot similarities between Geo-seq and ST-MAGIC (Response Fig. R5a)\nand quantitatively measure the spatial similarities using Pearson correlation and Structural\nSimilarities index (SSIM), a metric used for pattern comparison of spatial transcriptome4-6.\nAs shown, scatter plots and quantitative analyses reveal a high degree of concordance in Sox2\nexpression across five developmental stages (Pearson correlation index > 0.7, SSIM > 0.75)\n(Response Fig. R5a). Moreover, we also determined the global similarities of all zipcode genes,\nthose exhibiting spatially restricted expression during gastrulation, between Geo-seq and ST-\nMAGIC (Response Fig. R5b, c; Supplementary Table 2). Density plots and similarity metrics\n(Response Fig. R5b, c) demonstrate a strong overall agreement in gene expression\narchitecture, with high correlation and SSIM scores across the gene set. This global similarity,\ncombined with the JSS analysis (Supplementary Fig. 6e), provides convergent evidence that\nST-MAGIC reliably recapitulates the spatial transcriptomic landscape captured by Geo-seq.\nThese results have been included in the revised manuscript (Supplementary Fig. 6c, d, h).\nResponse Figure R5. Quantitative evaluation of gene expression pattern between Geo-seq and ST-MAGIC. a.\nScatter plot revealing the correlation of Sox2 expression between Geo-seq and ST-MAGIC. b. c. Density plots showing\nthe genome-wide spot-by-spot transcriptomic similarities between Geo-seq and ST-MAGIC based on metrics of\nPearson correlation (b) and SSIM index (c).\nFig 2g and h: The schematic is not useful.\nResponse: Thanks for the comment. We have revised the schematic to make a clear\nillustration for digital sampling from E7.5 right and left proximal mesoderm cells, respectively\n(Response Fig. R6a, b; Fig. 2g, h).\nResponse Figure R6. Schematic diagram showing the digital sampling from the right (a) and left (b) proximal\nmesoderm spots.\nLines 272-274: Are the differences statistically significant?\nResponse: Thanks for your insightful suggestion. We have calculated the average signal\ndensities for chromatin accessibility around the top 1000 right proximal mesoderm (RPM)-\nspecific (top) and left proximal mesoderm (LPM)-specific (bottom) peaks (log transformed\nsignal intensities around peaks center ± 500 bp) (Response Fig. R7a, b). Wilcox signed-rank\ntest analyses were applied to calculate the significance and the resulting P-values were labeled\nin the figure. As shown, the differences of chromatin accessibility between RPM and LPM are\nstatistically significant.\nThis figure has been included in the revised Supplementary Fig. 8b.\nResponse Figure R7. Statistical analyses for chromatin accessibility of RPM-specific peaks (top) and LPM-\nspecific peaks between left and right cell type counterpart. Wilcox signed-rank test analyses were used for\nsignificance calculation.\nED Fig 6g: How many gastruloids where generated and measured in this experiment? Did the\nauthors observe morphological differences between the cell lines or different generation\nefficiencies?\nResponse: Thanks for pointing this out. For the original Supplementary Fig. 6g, we utilized\nthree independent Neo_LRE KO mESC lines along with wild-type counterparts. At each\ntimepoint during gastruloid differentiation, eight gastruloids were collected and pooled together\nfor each genotype to determine Lefty2 expression (original ED Fig 6g). The detailed information\nregarding the number of replicates and pooling strategy has now been explicitly stated in the\nfigure legend of revised Supplementary Fig. 8k.\nAs to the morphological differences and differentiation efficiencies, we conducted additional\nindependent gastruloid differentiation experiments with triple biological replicates. Careful\nexamination of both undifferentiated mESCs and differentiating gastruloids revealed no\ndiscernible morphological differences between wild-type and Neo_LRE KO lines (Response\nFig. R8a, b). We also quantified the geometry (gastruloid length and width) of the gastruloids\nand the efficiencies of forming well-structured gastruloids. Statistical analysis showed no\nsignificant changes between the groups, with measurements from 16 wild-type and 27\nNeo_LRE KO gastruloids (Response Fig. R8c, d).\nThese results have been included in the revised Supplementary Fig. 8h-k.\nResponse Figure R8. Morphological characterization of Neo_LRE KO gastruloids. a. Bright field image showing\nthe morphologies of stem cell clone for both wildtype and Neo_LRE KO ESCs. b. c. Representative images (b) and\ngeometrical summary (c) of D5 wild-type and Neo_LRE KO gastruloids. Student’s t-test was used to calculate the\nsignificance. 16 wildtype gastruloids and 27 Neo_LRE gastruloids were used. Source data has been provided. d. Bar\nplot showing the differentiation efficiency for wild-type and Neo_LRE KO cells. Differentiation efficiencies were\ncalculated by summarizing the formation of well-structured gastruloids. Three replicates were performed for each\ngroup. Student’s t-test was used to calculate the significance.\nFigure 3\nLines 307-308: How were Etv2, Gata2 and Tal1 selected from the SCENIC+ candidates? Do\nthe other TF candidates that fit the TF expression hierarchy display the same paradigm of the\nupstream factor functioning as a priming factor responsible for opening chromatin for use by\nlater factors?\nResponse: Thanks for pointing this out. For your information, among the TF candidates in the\nSCENIC+ results (Fig. 3a), the molecular hierarchy from Etv2, Gata2, and Tal1 during\nhematoendothelial differentiation with ETV2 functioning as a priming factor responsible for\nopening hematoendothelial enhancer has been well studied (Response Fig. R9)7,8. Thus, we\nselected these factors to exemplify the phenomenon that the expression of TF shows high cell\ntype-specificity, while the distributions of TF targets (peaks and genes) often display shared\npattern among closely related cell types (hematoendothelial progenitors, blood progenitors 1,\nblood progenitors 2) (Fig. 3a; Line 324-329).\nWe thought it would be a reasonable hypothesis that ‘TF candidates that fit the TF expression\nhierarchy display the same paradigm of the upstream factor functioning as a priming factor\nresponsible for opening chromatin for use by later factors’. We acknowledged that a definitive\nconclusion cannot be made without sufficient experimental validation. In this manuscript, we\nhave experimentally validated the TF hierarchy of Eomes-Foxa2-Noto/Pou6f1, and found that\nearly and intermediated stage TFs (EOMES and FOXA2) could be responsible for opening the\nchromatin and then bound by later TFs (NOTO and POU6F1) during axial mesendoderm\nlineage development (Fig. 5 and Fig. 6).\nResponse Figure R9. Published study highlighting the TF hierarchy during hematoendothelial differentiation.\nLines 342-346: How many target genes and target peaks are there for NOTO?\nResponse: Thanks for the reminder. The number of NOTO target genes is 125 genes, and the\nnumber of NOTO target peaks is 207 peaks. These numbers have been included in the revised\nFig. 3e.\nFig 3e: It is difficult to evaluate the claim that target genes of NOTO and NOTO expression\nhave matched spatiotemporal distribution with the provided visualization. Moreover, the\nauthors should leverage their comprehensive integrated atlas, to evaluate how many\ngenes and their targets nominated by SCENIC+ display similar behavior. It is not clear why\nNoto is selected from the candidates SCENIC+ nominates for mesendoderm development.\nIt does not appear to be highlighted in ED Fig 7.\nResponse: Thanks for these insightful suggestions. Regarding the claim about the\nspatiotemporal concordance between Noto expression and its target genes, we have now\nquantitatively assessed the spatiotemporal alignment using both Pearson correlation and SSIM\nmetrics. As shown in the scatter plot (Response Fig. R10a; Supplementary Fig. 11c), the\nsimilarities between Noto expression and the spatial distribution of its predicted target genes\ngradually increase along with the progression of gastrulation. Specifically, at late gastrulation\nstages (from E7.0 to E7.5), when Noto becomes specifically expressed at the distal embryonic\ntip, its target genes exhibit a striking co-enrichment in the same region (Response Fig. R10a;\nSupplementary Fig. 11c). This dynamic alignment supports the biological relevance of the\npredicted NOTO regulon within the emerging axial structure.\nIn response to the suggestion to evaluate the broader landscape of TF–target coordination\nacross the atlas, we globally analyzed all transcription factors identified by SCENIC+. As shown,\nonly 15%-30% TFs exhibit high spatiotemporal concordance (Pearson correlation > 0.7) with\nthe expression patterns of their target genes, and about 5%-15% TFs exhibit strong spatial\nalignment (Pearson correlation > 0.7) with the accessibility patterns of their target peaks\n(Response Fig. R10b, c; Supplementary Fig. 11d, e).\nResponse Figure R10. Quantitative measurement of spatiotemporal relationship between TF and their targets.\na. Scatter plot showing the relationship between Noto expression and NOTO target genes expression. The Pearson\ncorrelation and SSIM index were also indicated at the top of each plot. b. c. Density plots showing the global correlation\n(b) and SSIM (c) between TF and TF targets.\nThis widespread asynchrony likely reflects the complexity of TF co-regulatory networks and the\nhierarchical kinetics of TF regulatory network during embryogenesis. These results underscore\na key strength of our ST-MAGIC (+) atlas: by integrating spatial, temporal, and multi-omic layers,\nwe are able to distinguish coherent regulatory relationships from those are temporally or\nspatially decoupled, thus providing a more nuanced view of developmental gene control.\nWe also acknowledge the confusion regarding the selection of Noto among SCENIC+\ncandidates, we sincerely apologize for the lack of clarity in our previous statement. Noto was\ninitially selected based on Fig. 3a, which lists cell type-specific TFs and Noto is specifically\nexpressed in Node cells. For the original Extended Data Fig. 7 (now Supplementary Fig. 10),\nit displays cell type-specific distribution of eRegulon based on target peaks specificity. We have\nnow supplemented a plot of eRegulon based on target genes specificity (Response Fig. R11;\nSupplementary Fig. 9), from which we can clearly detect that Noto target genes ranks the top.\nTogether, for Noto eRegulon, the expressions of Noto (Fig. 3a) and target genes (Fig. 3a;\nSupplementary Fig. 9) show specific distribution in node cells, while Noto target peaks are\npervasive accessible in the axial mesendoderm lineage (Fig. 3a; Supplementary Fig. 10).\nWe are sorry for the unclarified description here, and have now revised this in this revised\nversion (Line 320-324).\nResponse Figure R11. Identification of specific eRegulon based on specificity score of target genes. The\neRegulon rank in node cells was shown. The full result of eRegulon identified based on specificity score of target\ngenes has been shown in Supplementary Fig. 9.\nFigure 4\nLines 374-380: There are claims about the specific and global spatial expression patterns\nfor different factors. While these assessments can be made by eye using the overlay onto\nthe spatial coordinates, a quantitative assessment is more valuable. For example, if the\nauthors plot a histogram of the values for each spot, the factors that show specific\nexpression patterns and those that are global will be obvious. To be clear, I think both the\nplots the authors provide and the quantifications I recommend are needed.\nResponse: Thanks for the insightful recommendation. Here, we have profiled the scaled\nexpression of Wnt3 and Dkk1 using histograms. In consistent with the ST-MAGIC corn plot\nresults (Fig. 4b), we found that Wnt3 is predominantly expressed in the primitive streak and\nadjacent posterior mesoderm regions, while Dkk1 is specifically expressed in the anterior\nendoderm region (Response Fig. R12a, b). Moreover, we also profiled the histograms for\nNodal and Cer1 (Response Fig. R12c, d), which show consistent results as Fig. 4d.\nThese results have been incorporated in Supplementary Fig. 12b-e and Supplementary Fig.\n13b-e.\nResponse Figure 12. Spatial expression patterns of signaling ligand and antagonists. a. b. Histogram showing\nthe spatial expression patterns of Dkk1 and Wnt3 in E6.5 embryos (a). The comparison of Dkk1 and Wnt3 scaled\nexpression level at indicated embryonic regions was shown as box plot (b). Student’s t-test was used to calculate the\nsignificance. c. d. Histogram showing the spatial expression patterns of Dkk1 and Wnt3 in E7.5 embryos (c). The\ncomparison of Dkk1 and Wnt3 scaled expression level at indicated embryonic regions was shown as box plot (d).\nStudent’s t-test was used to calculate the significance. e. f. Histogram showing the spatial expression patterns of Cer1\nand Nodal in E6.5 embryos (e). The comparison of Cer1 and Nodal scaled expression level at indicated embryonic\nregions was shown as box plot (f). Student’s t-test was used to calculate the significance. g. h. Histogram showing the\nspatial expression patterns of Cer1 and Nodal in E7.5 embryos (g). The comparison of Cer1 and Nodal scaled expression\nlevel at indicated embryonic regions was shown as box plot (h). Student’s t-test was used to calculate the significance.\nSource data was provided.\nFig 4b,c: The claim is that Wnt3 and Dkk1 have overlapping distributions. This is difficult to\nassess by the provided panels. Moreover, can any statistics be applied to this analysis?\nAre the peaks enriched?\nResponse: Thanks for your kind suggestion. As mentioned above, Wnt3 and Dkk1 exhibit\nspecific and distinct distribution pattern (Response Fig. R12a, b; Fig. 4b).\nFor Fig 4b,c, we found that Wnt signaling target peaks were accessible in the embryonic region\nthat largely overlapped with Wnt3 distribution. Here, we incorporated the Pearson correlation\nand SSIM to quantitatively measure the similarities between these two plots. As shown in\nResponse Fig. R13a, b, the distribution of Wnt3 expression in ST-MAGIC exhibit high similarity\nwith that of targeting peaks in ST-MAGIC (+), and also with the linked genes expression in ST-\nMAGIC (+). This result has been incorporated in revised Supplementary Fig. 12k.\nAs to the peaks, we firstly compared published -Catenin ChIPseq datasets from in vitro ESCs\nand EpiLCs, an in vitro counterpart to cells in the mouse gastrula, from the published dataset\n(Supplementary Fig. 12g, h). Then, peaks enriched with specific -Catenin binding in EpiLCs\nwere selected for the ST-MAGIC (+) analyses here (Fig. 4c).\nResponse Figure R13. Spatial distribution similarities between Wnt3 expression and signaling targets. a.\ncorrelation between Wnt3 expression and Catenin target peaks. b. Correlation between Wnt3 expression and\nCatenin target peaks linked genes. Both structural similarity index and Pearson correlation analyses were applied.\nFigure 5\nLines 463-466: There is not a figure panel that shows the co-expression pattern for Shh,\nNoto, and Foxj1. I believe Fig 5e might be the panel but it does not have Shh or Noto\nhighlighted. It is difficult to assess the claim that these two distinct subtypes exist.\nResponse: Thanks for the kind reminder. As requested, we have profiled the co-expression\npattern of Shh, Noto, Foxj1 by ST-MAGIC (+) and also included the experimental results for\neach individual gene. As shown in Response Fig.14a, the expressed domains for these three\ngenes are regionalized to the distal tip and anterior part of the embryo. Detailed analyses for\neach gene revealed that these three genes were co-expressed in the distal tip (node), while\nonly Shh was also expressed in the anterior 2EA region, where the anterior mesendoderm cells\nlocate (Response Fig.14b). This was further validated by wholemount RNAscope staining\nexperiments (Response Fig.14c).\nTogether, these results revealed the presence of two distinct subtypes, anterior mesendoderm\ncells (Shh+Noto-Foxj1-) and node cells (Shh+Noto+Foxj1+) at the E7.5 mouse embryo. The co-\nexpression profile has been included in the revised Fig. 5b and Supplementary Fig. 14c.\nResponse Figure R14. Co-expression of typical axial mesendoderm markers demarcates the presence of\ntwo subtypes at the E7.5 mouse embryo. a. ST-MAGIC (+) visualization of the co-expression pattern of Shh, Noto,\nFoxj1 at the E7.5 mouse embryo. b. c. ST-MAGIC showing the expression of known markers for node and anterior\nmesendoderm cells (Shh, Noto, Foxj1) (b), which were validated by whole-mount RNAscope (c). AME is short for\nanterior mesendoderm.\nLine 514: What does a “highly specific distribution” mean for POU6F1? Can the authors\nsupport this claim with a statistic?\nResponse: Thanks for the reminder. We are sorry for the confusion. As shown in Response\nFig. R15a (Fig. 5l), we found that Pou6f1 and its target genesets show specific distribution in\nE7.5 node cells. Statistical analyses revealed that both Pou6f1 and its target genes show\nspecific distribution in Node cells (Response Fig. R15b, c).\nThese results have been included in the revised manuscript as\n‘Notably, among multiple known TFs GRNs, we identified a novel TF, POU6F1, which exhibits\nspecific distribution of both the TF itself (adjusted P-value equal to 8.4e-29) and its target genes\nin the node (adjusted P-value equal to 1.3e-55) (Fig. 5l).’\nin Line 538-541 and Supplementary Fig. 15b.\nResponse Figure R15. SCENIC+ analyses revealing the transcription factor hierarchy during axial\nmesendoderm development. a. Heatmap showing the TF eRegulon hierarchy during axial mesendoderm\ndevelopment. Pou6f1 regulon was highlighted in red. The numbers for selected TF-target genes and TF-target peaks\nwere labelled, where g is short for TF-target genes, p is short for TF-target peaks. Wilcox tests were used to calculate\nthe significance. b. c. Bar plots showing the specific enrichment of Pou6f1 expression (b) and its target genes (c) in\nnode cells. Wilcox tests were used to calculate the significance.\nFigure 6\nLine 531: “The target genes start to show abundant enrichment at the same region”.\nFrom inspection of Fig 6a, it appears to me that the POU6F1 target gene set expression\nappears before the expression of Pou6f1 with the former appearing at E7.25. Can the\nauthors explain the differential timing of these events?\nResponse: Thanks for the valuable reminder. We also noticed that POU6F1 target gene set\nexpression seems to appear at E7.25. After carefully checking the specific gene set of POU6F1\ntarget gene, we found that a certain number of target genes are co-regulated by multiple\nupstream transcription factors. For example, as shown in Response Fig. R16a, 31 out of 99\nPOU6F1 target genes are co-regulated by SP5. In comparison with the POU6F1 own regulated\ngenes, we found that the overlapping co-regulated genes and SP5 show earlier expression\npatterns (Response Fig. R16b, c). We proposed that these genes may influence the final\nvisualization results.\nNevertheless, to make a more precise description of the results, we have revised the original\ndescription into ‘the target genes show abundant enrichment at the same region.’ in Line 557-\n558.\nResponse Figure R16. Detailed investigation of the spatiotemporal distribution of Pou6f1 regulon. a. Venn\ndiagram showing the overlapping of target genes between Pou6f1 regulon and Sp5 regulon. b. c. ST-MAGIC\nvisualization of Pou6f1 own regulating genes (b) and overlapping genes (c).\nLines 532-535: The authors suggest that “early and intermediate stage TFs may function\nas priming factors by opening a broad range of chromatin regions”. This is a reasonable\nhypothetical model. Can the authors support or refute this claim with additional analysis?\nFor example, are the target peaks from early and intermediate stage TFs enriched for\nterminal TF targets? From Fig 5m,n, the number of peaks and genes targeted by Eomes\n(42p, 35g respectively) is far less than the number of peaks and genes targeted by Pou6f1\n(127p, 99g respectively). It's not clear to me how to incorporate these numbers into the\nproposed model.\nResponse: Thanks for pointing out this important issue. We highly appreciate the significance\nof seeing whether the target peaks from early and intermediate stage TFs are enriched for\nterminal TF targets.\nThe current SCENIC+ pipeline starts with identifying candidate enhancers from 10xscMultiome-\nsnATAC-seq dataset. The scATAC-seq results can only capture a certain portion of\nrepresentative targets due to data sparsity. Given this, towards comprehensive analyses of the\nrelationship between early/intermediate and late TFs, we used the ChIP-seq dataset to explore\nthe relationship between early/intermediate and late stage TFs. In this study, using the\npublished dataset for early/intermediate (EOMES/FOXA2) factor from in vitro differentiation\ncounterparts, we demonstrated that they are enriched at the late stage NOTO and POU6F1\ntarget sites (Fig. 6c).\nHowever, for the late stage TF (Noto and Pou6f1), neither published datasets for stage-\nmatched NOTO ChIP-seq and POU6F1 ChIP-seq nor ChIP-grade antibodies for these two TFs\nare available currently. Thus, we were unable to analyze the enrichment of late stage TFs on\nthe target peaks of early and intermediate stage TFs. Nonetheless, we analyzed the\nintersection between intermediate factor-Foxa2 and late stage TFs (Noto and Pou6f1). As\nshown, we found that the overlapped target peaks are largely related to the prospective function\nof axial mesendoderm (GREAT-phenotype analysis, Response Fig. R17a-d). These results\nfurther support our notion that early and intermediate stage TFs may function as priming factors\nby opening a broad range of chromatin regions.\nFor the number of targets for early/intermediate stage TFs and late stage TFs, as referenced\nin the current work model, the opening of chromatin targets is regulated by intricate TF\ncascades, and pioneer factors tend to initialize the opening of a small number of compacted\nchromatin and then expand the effects by enabling the action of nucleosome remodeler and\nother co-factors, thus facilitating the terminal binding and functioning of terminal stage TFs9,10.\nFor example, pioneer factor-PU.1 can locally open compacted chromatins, and then engage\nother proteins, including ATP-dependent chromatin remodeler-SWI/SNF, to enable further\nchromatin opening by late stage TFs (Response Fig. R17e)11,12. Thus, the target peaks of\npioneer factors (early/intermediate) are not necessarily larger than those of late-acting factors.\nResponse Figure R17. TF cascade. a. Venn diagram highlighting the overlapping of Foxa2 target peaks and Noto\ntarget peaks. b. GREAT phenotype analyses of the overlapping peaks in a. c. Venn diagram highlighting the overlapping\nof Foxa2 target peaks and Pou6f1 target peaks. d. GREAT phenotype analyses of the overlapping peaks in c. e.\nPublished model highlighting the TF cascade driving cell differentiation.\nLines 547-551, Fig 6c: The authors described that EOMES and FOXA2 ChIP-seq data\nshows enrichment of these TFs around the peaks of NOTO and POU6F1 targets. First, the\ntext does not describe the ChIP-seq data. The citations are in the legend and only by\nlooking at the papers could I discern what the ChIP-seq data is, which is from EBs and in\nvitro differentiation, respectively I believe. Second, there is no statistical test for the\nenrichment but I agree that by eye, there is a peak for both Eomes and FoxA2. However,\nit does not appear from the data presented in Fig 6c that the Eomes KO has an effect on\nthese peaks implying that the chromatin regulation may be controlled by other factors. I\nwould agree that F2-KO does have an effect but does not entirely explain the regulation at\nthose target loci. Will the authors please explain or modify the following statement in the\ntext: “EOMES and FOXA2 are likely to be responsible for pre-opening… target chromatin\nregions”.\nResponse: Thanks for the kind reminder. We deeply apologize for not explaining the detailed\ninformation of EOMES ChIP-seq and FOXA2 ChIP-seq in the text. For your information,\nEOMES ChIP-seq was from differentiated embryoid bodies (EB)13 and FOXA2 ChIP-seq was\nfrom differentiated mesendoderm cells (MEN)14, respectively. We have now included this\ninformation in the revised manuscript (Line 566-568).\nFor the statistical analyses, we qualified the signal intensities at the indicated peaks (log\ntransformed signal intensities around peaks center ± 500 bp). We also included the biological\nreplicates for WT ATAC-seq data as the control group. As shown, both Eomes knockout and\nFoxa2 knockout affect the accessibility of relevant chromatin regions (NOTO target peaks,\nPOU6F1 target peaks, WNT signaling responsive peaks, NODAL signaling responsive peaks),\nwhereas the two replicates in the control group exhibit consistently comparable levels\n(Response Fig. R18a-d). This has been included in the revised Supplementary Fig. 15j.\nNevertheless, for a better description of the results, we revised the statement into ‘EOMES and\nFOXA2 may play roles in pre-opening signaling responsive elements and NOTO and POU6F1\ntarget chromatin regions.’ in Line 578-579.\nResponse Figure R18. Statistical significance of chromatin accessibility change upon Eomes or Foxa2\nknocking-out. Wilcox test was used to measure the significance of signal intensities change.\nMinor points:\nFig 1e: I believe JSD index is incorrect. JSD measures the similarity between two\ndistributions. The Jensen-Shannon diversity index is a special case of JSD.\nResponse: Thanks for the insightful reminder. We are sorry for the inaccurate description. It is,\nas said, Jensen-Shannon Divergence index measures the divergence of two distributions. Here,\nwe used Jensen-Shannon Similarity (JSS) score, which was calculated based on JSD index,\nto measure the consistency between gene expression and linked peaks’ chromatin accessibility\n(Response Fig. R19). The detailed method for JSS score has been provided in Line 980-992.\nResponse Figure R19. The formula for JSS score calculation.\nBesides, we have also revised the description in the manuscript (Line162-165) as\n‘Jensen-Shannon Similarity (JSS) score (see Methods), which is calculated based on Jason\nShannon Divergence index, has been developed to measure the consistency between cell type\ntranscriptomic features and linked chromatin accessibility diversities. As shown, JSS score was\nhigher for BioCRE than for Signac and ArchR (Fig. 1e)’.\nED Fig 1g: needs a color legend for cell types.\nResponse: Thanks for the reminder. We are sorry for missing the color legend. We have now\nadded the color legend in the revised Supplementary Fig. 1h.\nED Fig 2b: The text describes cell types that are enriched at earlier or later time points. To\naid this visual assessment, it would be helpful if the x-axis had the same scale for each\ntime point.\nResponse: Thanks for the insightful suggestion. We have uniformed the x-axis for each time\npoint as Response Fig. R20. The revised figure has been included in the revised\nSupplementary Fig. 2b.\nResponse Figure R20. Bar plots showing the annotated cell type compositions in this study (top) and the\npublished dataset1 (bottom).\nED Fig 3h: units on y axis. How many replicates per condition?\nResponse: Thanks for the reminder. We have added the units on y axis, and the information\nof replicate number (n=4) for each genotype were included in the figure legend\n(Supplementary Fig. 5e).\nED Fig 10e: Needs a summary plot at the top of the heatmap such as in Fig 2i.\nResponse: Thanks for the insightful suggestion. We have added the summary plot as\nResponse Fig. R21.\nThe revised figure has been included in the revised Supplementary Fig. 13i.\nResponse Figure R21. Smad2/3 chromatin binding profile in ESC and EB cells. Left: heatmap plot of Smad2/3\nChIP-seq data; middle: the summary plot of signal intensities; right: GREAT annotation for each group of peaks.\nLine 486: It is not clear what Combinatorial pseudo-time based tracking means.\nResponse: Thanks for the comment. We apologize for this unclear wording. For the original\nline 486 (Fig. 5h-k), we have tracked the gene expression dynamics and the chromatin\naccessibility of associated peaks following the pseudotime ordering.\nTo make a clear description, we have revised the text as\n‘Pseudotime-based tracking of gene expression cascades and associated chromatin peaks\nrevealed three successive stages (G1, G2, G3) enriched in the early, intermediate and terminal\npopulations, respectively (Fig. 5h-k).’\nThis has been included in the revised manuscript (Line 511).\nED Fig 12f: Caption missing explanation of acronyms in panel (EB – embryoid body, MEN\n- ?)\nResponse: Thanks for the kind reminder. We are deeply sorry for missing explanations for\nthese abbreviations. The acronym for EB represents embryoid bodies, and the acronym for\nMEN represents mesendoderm. We have now added these in the revised figure legend (Line\n566-568, figure legend for the revised Supplementary Fig.15g).\nLine 569, ED Fig 13h: I realize the numbers are small, but it would nonetheless be useful\nfor the authors to derive a likelihood (p-value) for the axis turning randomization.\nResponse: Thanks for the insightful suggestion. In this revision, we have collected more\nembryos and summarized the phenotypes for each genotype (WT, Noto KO, and Pou6f1 KO).\nThe chi-squared test was performed to determine the statistical significance. As shown in\nResponse Figure R22, both Noto KO and Pou6f1 KO embryos show significant differences in\naxis turning direction in comparison with the WT control embryos. This result has been included\nin the revised Supplementary Fig. 16h, i.\nResponse Figure R22. Both Noto KO and Pou6f1 KO embryos show significant differences in axis turning\nin comparison with the WT control. a. Table showing the specific numbers and percentages for each phenotype.\nb. Stacked bar plot highlighting the axis turning differences among three groups.\nReviewer #2 (Remarks to the Author):\nIn this manuscript, Yang and colleagues present a spatiotemporal dataset capturing mouse\ngastrulation, generated using 10x Genomics Multiome data. The study involves a\nsignificant amount of data integration and projection, incorporating external datasets such\nas Geo-seq and ChIP-seq. To enable gene regulatory network (GRN) analysis, the authors\ndeveloped a new computational tool called BioCREs, which they claim performs better\nthan existing pipelines. They also adapted the Tangram algorithm—originally designed for\nhigh-resolution spatial transcriptomics—for use with their data to align it with Geo-seq\nreferences. Through this combined approach, they aim to map spatial and temporal\nregulatory changes during gastrulation.\nThe scope of this work is impressive, and the integration of multiple data types is interesting.\nThe authors also acknowledge the limitations of current spatial resolution and mismatched\ncell type references, which are important caveats. I appreciate the validation experiments\nincluded, which lend support to some of their conclusions. That said, I do have several\nconcerns, especially regarding the BioCREs pipeline and the modified Tangram tool.\nResponse: We highly appreciate the reviewer’s positive assessment and constructive\nsuggestions on our manuscript. In this revision, we have carefully revised the manuscript,\nespecially for the benchmarking of BioCRE pipeline and comparing the ST-MAGIC results with\nStereo-seq results. Generally, these efforts further strengthen our results.\nMajor Concerns:\n1. The claim that BioCREs outperforms Signac and ArchR is based on a single\nperformance index. This is insufficient to justify the claim. Benchmarking against\nestablished regulatory element datasets (e.g., ENCODE) is necessary to validate accuracy\nand specificity.\nResponse: Thanks for this valuable suggestion. We fully agree that a more comprehensive\nbenchmarking strategy is essential to substantiate our claim that BioCRE outperforms existing\ntools.\nIn this revision, to systematically benchmark the performance of BioCRE, we have conducted\na series of systematic evaluations to rigorously compare the performance of BioCRE against\nthe state-of-the-art methods: Signac and ArchR.\n1. The set characteristics of predicted gene-to-peak linkages:\nUsing our multi-omic dataset from the mouse gastrula, we evaluated the number, biological\nrelevance, and cell type specificity of predicted gene-to-peak linkages. As shown, BioCRE\nidentified 44,644 gene-to-peak linkages, with a median gene-to-peak distance of 127,175 bp\n(Response Fig. R23a, b). Comparative analyses revealed that 28,472 linkages (64%) were\nspecifically identified in BioCRE (Response Fig. R23c, d). Gene ontology analyses of genes\nassociated with these BioCRE-specific linkages (28,472 linkages) revealed significant\nenrichment in embryo development, including heart development, liver development, and\nnervous system development (Response Fig. R23e). This result suggests that these BioCRE-\nspecific linkages could be responsible for the regulation of embryo development related genes.\nTo assess the cell type specificity of predicted linkages, we analyzed the accessibility of peaks\nlinked to cell type-specific genes identified by each tool (Response Fig. R23f, g). Using UMAP\nvisualization and quantifying clustering coherence with the Silhouette score, we found that\nlinkages predicted by BioCRE gets a higher score compared to those from Signac and ArchR,\nindicating a more effective capability for BioCRE to capture cell type-specific linkages\n(Response Fig. R23f, g).\nTogether, these results demonstrate that BioCRE could be adept at identify cell-type-specific\nlinkages, which are likely to play key roles in orchestrating embryo development.\n2. The computational efficiency and robustness of predicting linkages:\nWe assessed the computational efficiency and robustness of different methods under varying\nsampling conditions. Overall, we found that the built-in deep-learning bi-orientation regression\nmodel for BioCRE does not incur excessive run-time costs, showing comparable performance\nto ArchR (Response Fig. R23h). To evaluate the robustness of linkage identification, we\nchecked the consistency of predicted linkages across three distinct down-sampling scenarios:\n In the first scenario, where both RNA and peak reads were down-sampled concurrently,\nBioCRE demonstrated the highest recovery rate of linkage identification (Response Fig.\nR23i). To specify, 40.9% linkages (31,911 out of 78,109) are constantly identified in\nBioCRE, while only 8.2% linkages (3,677 out of 44,705) for Signac and 5.7% linkages\n(4,479 out of 79,067) for ArchR show consistent results among variable sample sizes\n(Response Fig. R23i-1).\n In the second scenario, where only RNA reads were down-sampled while peak reads\nremain unchanged, we found that 44.6% linkages (21,406 out of 48,043) are consistently\ncaptured for BioCRE (Response Fig. R23i-2). In the meantime, 42.9% linkages (18,335\nout of 42,937) for Signac and 61.9% linkages (38,510 out of 62,226) for ArchR are stably\nidentified (Response Fig. R23i-2).\n In the last scenario, we specifically down-sampled peak reads and then performed linkage\nidentification (Response Fig. R23i-3). Consequently, we found that 52.1% linkages\n(49,476/95,042) for BioCRE are captured across all down-sampling groups (Response\nFig. R23i-3). In contrast, only 24.1% linkages (9,122/37,894) for Signac and 46.2%\nlinkages (28,388/61,422) for ArchR show the same pattern (Response Fig. R23i-3).\nCollectively, these analyses demonstrate that BioCRE achieves a favorable balance between\ncomputational efficiency and robustness, exhibiting consistently high reproducibility in linkage\nprediction across diverse and reduced sequencing depths, particularly under conditions of\nlimited chromatin accessibility data (Scenario 1 and 3), highlighting its reliability in practical,\nresource-variable setting.\n3. The reliability and accuracy of gene-to-peak linkages:\nGiven the absence of a genome-wide gold standard for gene-to-peak linkages in the mouse\ngastrula, we turned to external single-cell multiome datasets with matched, high-confidence\npromoter-enhancer interaction data to benchmark performance. Specifically, we applied\nBioCRE, Signac and ArchR to publicly available single-cell multiome data from human PBMC\ncells and Gm12878 cells, using promoter-capture Hi-C (PCHiC) data as an orthogonal\nreference for true regulatory interactions (Response Fig. R23j).\nTo enable a comprehensive and balanced evaluation, we defined both positive and negative\ncontrols: positive linkages were inferred from PCHiC-confirmed interactions, while negative\nlinkages were established as peak-gene pairs showing no chromatin interactions within a 1 Mb\nrange in the PCHiC data. We evaluated the performance using multiple following metrics:\n The metrics of Macro-F1 score across all linkages showed that BioCRE consistently\nachieved higher scores than both Signac and ArchR across all tested cell types, with strong\nreproducibility across replicates (Response Fig. R23j-1). This indicates that BioCRE’s\npredictions more closely align with the experimentally determined chromatin architecture.\n The Macro-F1 score for linkages of all three tools shared genes revealed that BioCRE\noutperformed both Signac and ArchR in most cell types except in CD16 monocyte, , where\nits advantage was less pronounced (Response Fig. R23j-2).\n The True Positive Rate (TPR) for linkages unique to each method was also assessed.\nBioCRE-specific linkages show comparable or higher level of TPR than those of Signac\nand ArchR (Response Fig. R23j-3), suggesting that BioCRE-specific predictions are not\nonly robust but also biologically credible.\nThus, through incorporating external dataset and gold standard chromatin interaction dataset,\nwe demonstrated that BioCRE exhibits consistent superiority in linkage identification than\nSignac and ArchR.\nThese results have been included in the revised manuscript (Line 162-176; Supplementary\nFig. 4), the detailed method has been provided in Line 982-1038, and we also acknowledged\nthe significance of systematically benchmark of BioCRE performance and discussed about the\ncurrent limitations in Line 623-632, Line 703-705.\nResponse Figure R23. Benchmark of BioCRE. a. Bar plot showing the number of gene-peak linkages identified by\nBioCRE, Signac, and ArchR. b. Histogram showing the median gene-peak distance for linkages identified by BioCRE,\nSignac, and ArchR. c. Venn diagram showing the intersection of predicted gene-peak linkages. d. Bar plot highlighting\nthe ratio of non-overlapping linkage ratios for each tool. e. Gene ontology analyses for BioCRE-specific captured\nlinkages associated genes. f. UMAP visualization based on the chromatin accessibility of cell type-specific genes\nassociated peaks. g. Silhouette score evaluating the cell type distribution specificity. h. Line plot depicting the runtime\nfor each tool with variable proportion of reads. BioCRE was running on a single GPU (A100), while ArchR and Signac\nwere executed using 20 CPU cores for parallel computing, respectively. i. The effects of sequencing reads down-\nsampling on linkage composition consistency for each tool. The constant linkage numbers for each scenario are\nhighlighted in red. j. The prediction and validation of gene-peak linkages using external multiome dataset and gold\nstandard promoter-enhancer interaction dataset.\n2. Distal vs. Proximal regulatory region comparison: The comparison between \"Distal High\"\nand \"Proximal\" groups is biased by excluding the \"Distal Low\" group. A more appropriate\ncomparison would be between all distal regions and proximal ones, with proper statistical\ntesting (e.g., p-values) to support conclusions.\nResponse: Thanks for your constructive suggestion. As recommended, we have directly\ncompared the JSS score between proximal and distal linked peaks—including both 'distal high'\nand 'distal low' groups—thereby providing a more balanced and representative assessment.\nStatistical testing (Wilcoxon rank-sum test) confirms that, on average, proximal peaks exhibit\nsignificantly higher JSS scores than distal peaks (Response Fig. R24a), supporting their\nstronger association with gene expression variation.\nAdditionally, we retained the comparison with the \"distal high\" group as a secondary analysis\nto emphasize the particularly strong regulatory potential of this subset, clearly distinguishing it\nfrom the broader distal category. We observed that a considerable number of distal linked\npeaks exhibit higher JSS score than those of the proximal peaks (Response Fig. R24b),\nhighlighting the functional importance of certain distal regulatory elements during gastrulation.\nThese updated analyses have been incorporated into the Fig. 1f and Supplementary Fig. 5c.\nResponse Figure R24. Bar plots highlighting the JSS score measuring the consistency between gene\nexpression and linked peaks. a. The comparison between proximal and distal linked peaks. b. The distal peaks\nwere subdivided into distal high and distal low group based on the comparison with proximal peaks. Distal peaks were\nannotated as peaks located outside TSS ± 2 kb, proximal peaks were annotated as peaks located within TSS ± 2 kb.\nWilcox tests were used to calculate the significances.\n3. Insufficient benchmarking for BioCREs and modified Tangram: The manuscript lacks\ntechnical detail and performance benchmarking for both pipelines. More rigorous testing is\nneeded to demonstrate their robustness and reliability.\nResponse: Thanks for raising this important issue. We are deeply sorry that the initial\nsubmission lacked sufficient technical details and comprehensive benchmarking.\nRegarding BioCRE, we have now conducted an extensive performance evaluation, as detailed\nin our response to comment #1 and illustrated in Response Figure R23. This includes\ncomparisons with established methods—Signac and ArchR—across multiple metrics such as\ncell type specificity (Silhouette score), computational efficiency, robustness to various sampling\ndepth, and accuracy against orthogonal chromatin interaction data (Macro-F1 and TPR scores).\nThe results consistently demonstrate that BioCRE, built on a bidirectional regression framework,\noutperforms the other tools in key aspects relevant to gene-to-peak linkage prediction in\ndevelopmental contexts. These analyses have been integrated into the revised manuscript\n(Line 162-176; Supplementary Fig. 4) to enhance methodological transparency and rigor.\nWith respect to the use of Tangram for spatiotemporal mapping, we clarify that our approach\ndoes not involve fundamental modification to the Tangram algorithm itself. Instead, we have\nadopted a more refined and biologically informed input strategy, as recommended in the official\nTangram documentation (https://github.com/broadinstitute/Tangram). The refinement strategy\nof gene selection for Tangram has been benchmarked in a recent study15. Specifically, rather\nthan using all genes shared between the Geo-seq and 10x scMultiome datasets, we selected\nstage-specific \"zipcode\" genes — those showing spatially restricted expression patterns during\ngastrulation — to improve mapping accuracy. This refinement aligns with best practices and\nenhances the biological relevance of the spatial inference (Supplementary Fig. 6c-i). We are\nsorry that this aspect was not clearly described in the original submission. We have now\nprovided the complete list of genes used (Supplementary Table 2) and updated the method\nwith a detailed and accurate description of using tangram (Line 859-881).\n4. Incorporating data from high-resolution platforms such as Visium HD or Stereo-seq\nwould significantly strengthen the spatial claims.\nResponse: Thanks for the valuable suggestion. We fully agree that incorporating high-\nresolution spatial data, such as Visium HD or Stereo-seq, would strengthen the spatial claims\nof the current study.\nIn this revision, we have incorporated a recently published Stereo-seq dataset (SEU-3D)16\ncovering the mouse early organogenesis embryos, to which we contributed, for direct\ncomparison with our ST-MAGIC atlas. As shown, to check the consistency of spatial results\nbetween these two datasets (SEU-3D and ST-MAGIC), we first spatial-temporally aligned the\ntwo datasets, then determined the transcriptional correlation of cell types and spatial\nconcordance of cell type marker geneset (CTGs) expression, and finally checked the\nreproducibility of key spatial findings reported in the current study (Response Fig. R25a).\nGenerally, spatial transcriptomic data from the E7.5 stage were used, the spatial coordinate\nalong z-axis (proximal-distal axis) and the cell types annotated in both datasets were aligned\naccording to their respective sampling strategies (Response Fig. R25b). Pearson correlation\nanalyses of cell type transcriptome reported that most cell types were well correlated\n(Response Fig. R25c). Moreover, the spatial distribution of CTGs were also highly consistent\n(Response Fig. R25d). We first present the AUC scores of the top 500 marker gene sets for\neach cell types (illustrated for VE, SE, and CE) in SEU-3D, where the spatial distribution of\nthese scores aligns well with the known anatomical localization of each cell type (first column).\nApplying the identical analytical to ST-MAGIC reveals a highly similar spatial pattern (second\ncolumn), demonstrating the consistency of marker gene enrichment across reconstruction\nmethods. To enable a more granular comparison, we then performed a layer-by-layer analysis\nrelative to Stereo-seq and GEO-seq. For each cell type, we extracted corresponding single-cell\nprofiles from Stereo-seq and matched them with spot-level single-cell compositions in each\nlayer of GEO-seq. The similarity between these matched cell populations was assessed using\nthe marker gene sets. As shown in the results, high similarity is observed at spatially\ncorresponding locations (third column), while gray regions indicate positions where data were\nnot jointly present in both datasets. For example, the expression of marker geneset for visceral\nendoderm cells were highly enriched in outer layer cells of SEU-3D map and in EA and EP\nregions of the ST-MAGIC atlas (Response Fig. R25d). Spatial correlation analysis confirmed\nthat this enrichment is anatomically coherent, as expression of these markers correlates\nobviously only in regions where visceral endoderm is known to reside (Response Fig. R25d).\nThen, we checked the spatial distribution of specific cell types and genes between SEU-3D and\nST-MAGIC atlas. As shown in Response Fig. R25e (Fig. 2b in the manuscript), mesenchyme\ncells were located at the proximal posterior region and paraxial mesoderm 1 cells were\ndistributed at the anterior medial mesoderm region according to both ST-MAGIC atlas and\nexperimental results (Fig. 2b, c). Here, we found that mesenchyme and paraxial mesoderm\ncells exhibit similar spatial patterns in SEU-3D atlas. Furthermore, the expression patterns of\ncanonical markers such as Hand1 (mesenchyme) and Irx5 (paraxial mesoderm) align closely\nbetween the two platforms, further corroborating the accuracy of our ST-MAGIC atlas.\nWe acknowledge that direct integration is challenged by inherent technical differences: SEU-\n3D achieves near single cell bin resolution but with lower transcriptomic depth (~400–600\ngenes per bin), whereas ST-MAGIC spots capture richer molecular profiles (~3,500–4,500\ngenes per cell) but at lower spatial resolution (~5 - 10 cells per spot) (Response Fig. R25f).\nDue to this disparity in molecular sensitivity, some spatially-restricted, low-abundance\ntranscripts—such as Wnt3, Nodal and Noto, which are clearly resolved in ST-MAGIC and\nsupported by prior experimental data—do not exhibit evident expression patterns in SEU-3D\n(Response Fig. R25g).\nDespite these technical constraints, the overall concordance in cell type localization and gene\nexpression architecture between SEU-3D and ST-MAGIC strongly supports the robustness of\nour spatial reconstructions. We have included the results of cell type spatial allocation in SEU-\n3D to demonstrate the credibility of ST-MAGIC in revised Supplementary Fig. 7f and Line\n263-266, and also discussed this important issue as\n‘Transcriptome coordinate profiling using Stereo-seq17 or MERFISH18, incorporation of\nadditional data types, such as histone modifications, proteomics, and metabolomics, and\noptimizations of bioinformatic toolkits would further enhance the molecular understanding of\nearly embryogenesis.’\nin the revised manuscript (Line 696-700).\nResponse Figure R25. Incorporation of high resolution spatial transcriptomic data. a. Schematic overview of\nthe work flow to check the spatial consistency between Stereo-seq dataset (SEU-3D) and ST-MAGIC. b. Temporal,\nspatial, cell type alignment between these two datasets. c. Pearson correlation for cell types annotated in both\ndatasets. d. Spatial concordance of cell type distribution between these two datasets. Three representative cell types\nwere shown. top: Visceral endoderm; middle: Surface ectoderm; bottom: Caudal epiblast. e. The consistency of\nspecific spatial claims between SEU-3D and ST-MAGIC atlas. f. Bar plot highlighting the detected nFeature_RNA\nper cell in SEU-3D and ST-MAGIC. g. The expression patterns of Wnt3 (left), Nodal (middle), and Noto (right) as\nvisualized in experimental results, ST-MAGIC, and SEU-3D.\nMinor Comments:\n• The manuscript is dense, especially in later sections discussing biological processes.\nStreamlining these parts would improve accessibility.\nResponse: Thanks for the comment. In this revision, we have revised the original manuscript\nto make more streamlined descriptions.\n• Some extended data figures are poorly organized and should be revised for clarity.\nResponse: Thanks for the kind reminder. In this revision, we have revised the organization of\nthe original extended figures to improve the clarity.\n• Clarity in Extended Data Fig. 2a: The plot is overcrowded due to overlapping cells, making\nit difficult to interpret. A side-by-side comparison with reference data in separate panels\nwould improve readability.\nResponse: Thanks for the insightful suggestion. We have now revised the presentation as in\nResponse Fig. R26. This result has been included in the revised Supplementary Fig. 2a.\nResponse Figure R26. Consistent transcriptome architectures between our dataset and publicly available\ndataset. a. UMAP co-embedding of both our dataset and the published transcriptomic atlas. b. A side-by-side comparison\nbetween the transcriptomic dataset (top) and the published single-cell data (bottom).\n• Missing data in Extended Fig. 4: The claim that ST-MAGIC has higher JSD values than\nGeo-seq is not supported by the figure. This should be clarified or removed.\nResponse: Thanks for the reminder. We are sorry for the unclarified description in the original\ndescription regarding Extended Fig. 4 (now Supplementary Fig. 6e). The JSD scores (now\nrevised as JSS score, Line 162-164) were in fact used to assess the consistency of gene\nexpression architecture between ST-MAGIC and the Geo-seq dataset.\nTo correct this, we have revised the main text to provide a more accurate and precise\ninterpretation. The updated description now reads:\n‘To assess the efficacy of the ST-MAGIC, we found that cells were mapped evenly across\nspatial spots in the constructed atlas (Supplementary Fig. 6b), the global gene distributions\nacross five stages in ST-MAGIC were well correlated with Geo-seq results (Supplementary Fig.\n6c, d), JSS scoring also revealed a consistent gene expression architecture between ST-\nMAGIC and the Geo-seq dataset (Supplementary Fig. 6e).’\nThis has been incorporated in the revised manuscript (Line 239-243).\nReferences\norganogenesis. Nature 566, 490-495 (2019).\nDevelopment 42, 293-303 (1977).\nnetwork driving blood and endothelial development. Curr Biol 18, 1234-1240 (2008).\nfrom error visibility to structural similarity. IEEE Trans Image Process 13, 600-612\n(2004).\nexpression via fusing spot-based spatial transcriptomics, location, and histology using\nGCN. Brief Bioinform 26 (2024).\n6533 (2025).\nhematoendothelial fate commitment. Cell Rep 42, 112665 (2023).\nrecruiting the repressor REST to restrict alternative lineage commitment. bioRxiv,\n595971 (2024).\nchromatin opening. Nat Commun 10, 3807 (2019).\nregulators of chromatin and gene expression. Nat Rev Genet 24, 809-815 (2023).\ntargeted ATP-dependent nucleosome remodeling. Nat Struct Mol Biol 30, 31-37 (2023).\nTrends Genet 40, 134-148 (2024).\nsegregation by changing the chromatin state. Nature Cell Biology 21, 1518-1531\n(2019).\nthe pioneering activity of Foxa2. Nucleic Acids Res 47, 9069-9086 (2019).\nmapping. Bioinformatics 41, i552-i560 (2025).\nCell (2025).\nnanoball-patterned arrays. Cell 185, 1777-1792.e1721 (2022).\nMERFISH. Nature 598, 137-143 (2021)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have improved the manuscript and addressed many of my comments. I commend the authors again for including many data panels that make it possible to evaluate the data and claims. Nonetheless, I have some remaining concerns.\n\nFigure 1:\nPlease include the number of embryos into the table in ED Fig 1e. I’m surprised that the number of cells in each embryo included in this study is not available and a literature based reference is given, especially because 10x genomics recommends loading numbers as part of the protocol. I appreciate that a reference number is included now as I believe it is important to assess sampling.\n\nED Fig 3a: From my understanding of the authors’ response, I believe that ED Fig 3a is mislabeled as it does not represent differential gene expression (DEG). Differential gene expression is an analysis where two or more groups are compared against each other. From the description, I believe that ED Fig 3a is Cell type-specific genes and ED Fig 3b is Cell type-specific accessible peaks. It is not clear from the manuscript how these genes are selected, which could be from DEG analysis. The Methods state “marker genes specific to each cell type were mainly used to annotate clusters…” but it is not stated how the marker genes were identified. My apologies if I missed it. Will the authors please explain how the genes in ED Fig 3 are selected?\n\nBioCRE:\nThe authors use GO analysis to show that BioCRE-specific linkages are enriched in embryo development. Since the purpose of the analysis is to benchmark against the other methods, would the authors please also show the GO analysis results for Signac and ArchR?\n\nThe downsampling analysis revealed interesting behavior from BioCRE. In Scenario 1, BioCRE returns approximately the same number of linkages for all of the downsampled sets, with a trend where more data leads to a smaller set. This is in contrast to Signac and ArchR which show an increased set size with increased data input. Scenarios 2 and 3 show that this difference in behavior is driven by the ATAC-seq data. The behavior, to me, implies that Signac and ArchR are prioritizing precision while BioCRE prioritizes recall. The desire to prioritize precision is sometimes preferred due to the need to prioritize targets for experimental validation.\n\nThe next analysis seeks to assess the tools against a gold standard. Will the authors please show the results in terms of precision and recall rather than the combined Macro-F1 score? Moreover, it is not clear to me why results are shown relative to BioCRE’s scores rather than, for example, showing each tool’s precision score. Given my comment above, I am curious how precision and recall behave for the downsampled data sets.\n\nFigure 5.\nI appreciate that the authors included a panel showing the enrichment of cells that are Shh+, Noto+, Foxj1+ (ED F14c). Indeed, it does look enriched in the position of the node. Nonetheless, I do not believe the new plots are sufficient to support their claim that there are two distinct subtypes (Shh+Noto-Foxj1- and Shh+Noto+Foxj1+). From the scRNA-seq data, what proportion of the cells fall into these two categories. Presumably they form two subclusters if the authors are confident in defining the subtypes. Is there a plot they can provide to support this claim?\n\nFig 5l. I believe the authors mean “specific expression” and not “specific distribution”\n\nFigure 6.\nThank you for changing the wording in the manuscript to better match the figure panel. Given the complexity about the dynamics of the system, I believe it will be worthwhile to include a quantitative figure panel better illustrating the temporal dynamics of the regulation. Can the authors please include bar plots or line plots?", "reviewer_2": "I am satisfied with the revised version."}, "author_response": "We highly appreciate the reviewer’s recognition of our efforts in addressing the previous comments. Here, we have addressed the remaining concerns in a point-by-point manner below.\n\nFigure 1:\nPlease include the number of embryos into the table in ED Fig 1e. I’m surprised that the number of cells in each embryo included in this study is not available and a literature based reference is given, especially because 10x genomics recommends loading numbers as part of the protocol. I appreciate that a reference number is included now as I believe it is important to assess sampling.\n\nResponse: Thank you for this insightful comment. We have included the number of embryos into the revised ED Fig. 1e (Response Fig. R1). We are sorry for the imprecise description in our manuscript and appreciate for the opportunity to clarify this issue. As the reviewer correctly pointed out, here to meet the loading requirements for 10x genomics, embryos at each developmental stage were pooled and dissociated into single-cell suspensions together. Therefore, while we can estimate the total number of cells per embryonic stage, it is technically challenging to determine the exact number of cells originating from each individual embryo. Here, we used the theoretical cell numbers for embryos at each embryonic stage based on published literatures1,2 and the cell numbers captured for each embryonic stage to reflect the sampling coverage for our multi-omic dataset. We have added the specific details for sample preparation in the revised method (Line 737-739) and revised the figure legend for Supplementary Fig. 1e to specify the definition for each number.\n\nResponse Figure R1. Top: Bar plots showing the number of cells captured in this study compared to the theoretical number of cells in the corresponding embryonic stage. Bottom: table summarizing the embryo numbers and number of single cells sampled from the embryos at each stage. Theoretical cell numbers represent the theoretical cell numbers for embryos at each embryonic stage inferred from published literatures; Numbers of embryos pooled represent the numbers of embryos pooled together for each stage before single cell dissociation; Cell numbers captured represent the final cell numbers captured with multi-omic data for each embryonic stage.\n\nED Fig 3a: From my understanding of the authors’ response, I believe that ED Fig 3a is mislabeled as it does not represent differential gene expression (DEG). Differential gene expression is an analysis where two or more groups are compared against each other. From the description, I believe that ED Fig 3a is Cell type-specific genes and ED Fig 3b is Cell type-specific accessible peaks. It is not clear from the manuscript how these genes are selected, which could be from DEG analysis. The Methods state “marker genes specific to each cell type were mainly used to annotate clusters…” but it is not stated how the marker genes were identified. My apologies if I missed it. Will the authors please explain how the genes in ED Fig 3 are selected?\n\nResponse: Thank you for pointing out this mistake. The genes and peaks displayed in the heatmaps were selected using the COSG algorithm3, which identifies genes or peaks that are specifically expressed or accessible in one cell type based on cosine similarity. As suggested, we have revised the figure labels and related captions as ‘cell type-specific genes’ and ‘cell type-specific accessible peaks’ for ED Fig 3 (Response Fig. R2) to accurately reflect the content. For the marker genes of each cell type, we have referenced the marker gene sets from the published gastrulation transcriptome atlas2,4. We have now clarified the details of this analysis in the figure legend of Supplementary Fig. 3a, b in the Methods section (Line 1101-1103) as follows:\n“The integration of cluster-specific genes (defined by COSG algorithm) and marker genes from the published gastrulation transcriptome atlas were mainly used to annotate the clusters.”\n\nResponse Figure R2. The scMultiome atlas delineates the molecular features of the annotated gastrulating cells. a. Heatmap showing the cell type-specific genes across the 31 annotated cell types. Representative genes were highlighted at the right. The defined genes for each cell type were visualized by profiling the scaled gene expression level. Color legend represents the scaled gene expression level. Cell type-specific genes were identified through COSG algorithm. b. Heatmap showing the cell type-specific accessible peaks across the 31 annotated cell types. Color legend represents scaled values of peak count data. Cell type-specific peaks were identified through COSG algorithm. c. Heatmap showing the TF-motifs enrichment across the 31 annotated cell types. The representative motifs were highlighted at the right. Color legend represents scaled values of motif activity scores generated by chromVAR.\n\nBioCRE:\nThe authors use GO analysis to show that BioCRE-specific linkages are enriched in embryo development. Since the purpose of the analysis is to benchmark against the other methods, would the authors please also show the GO analysis results for Signac and ArchR?\n\nResponse: Thank you for this insightful reminder. As suggested, to benchmark the performance of development-related gene ontology terms enrichment, we have included the gene ontology results for Signac and ArchR (Response Fig. R3a, b). For this analysis, gene ontology terms with ‘development’, ‘stem cell’ and ‘signaling pathway’ are selected as development-related terms. As shown, most of the development-related GO terms show high significant levels for BioCRE predicted linkages (Response Fig. R3a, b). These results indicate that BioCRE may more accurately anchor specific linkages governing developmental gene expression compared with ArchR and Signac. This has been included in the revised Supplementary Fig. 4e, f.\n\nResponse Figure R3. The gene ontology analyses for tool-specific linkages related genes. a. Heatmap showing the enrichment of development-related gene ontology terms for specific linkages predicted by BioCRE, ArchR, and Signac, respectively. b. Bar-plot showing the significance levels of development-related gene ontology terms for each method.\n\nThe downsampling analysis revealed interesting behavior from BioCRE. In Scenario 1, BioCRE returns approximately the same number of linkages for all of the downsampled sets, with a trend where more data leads to a smaller set. This is in contrast to Signac and ArchR which show an increased set size with increased data input. Scenarios 2 and 3 show that this difference in behavior is driven by the ATAC-seq data. The behavior, to me, implies that Signac and ArchR are prioritizing precision while BioCRE prioritizes recall. The desire to prioritize precision is sometimes preferred due to the need to prioritize targets for experimental validation.\n\nResponse: Thank you for this insightful observation that BioCRE returns consistent results for linkage prediction across all down-sampled sets, whereas Signac and ArchR yield progressively larger sets of predicted linkages as data volume grows (Scenario 1). It seems that the linkages identified by Signac and ArchR are strongly dependent on sequencing depth, while the majority of linkages detected by BioCRE remain robust across varying sequencing depths—indicating that BioCRE would be less sensitive to sequencing depth bias (consistent with our findings on Scenario 2 and 3). In-depth analysis reveals that Signac and ArchR rely on correlation-based approaches, this behavior aligns with previous studies demonstrating that such methods might be susceptible to technical artifacts, particularly sequencing depth effects in sparse and compositional data5,6. In contrast, BioCRE employs a regression-based framework, which has been shown to be relatively more robust to these technical variations due to its ability to model biological dependencies while explicitly accounting for potential confounding factors5.\nWe highly appreciate your thought-provoking interpretation regarding a potential trade-off between precision and recall. However, given the absence of a reference-standard orthogonal dataset for the mouse gastrula, we are unable to definitively assess the precision and recall of the predictions here. Alternatively, we have examined and specified the precision and recall results using published datasets in the following section (Response Fig. R4, 5).\n\nThe next analysis seeks to assess the tools against a gold standard. Will the authors please show the results in terms of precision and recall rather than the combined Macro-F1 score? Moreover, it is not clear to me why results are shown relative to BioCRE’s scores rather than, for example, showing each tool’s precision score. Given my comment above, I am curious how precision and recall behave for the downsampled data sets.\n\nResponse: Thank you for this insightful comment. We agree that precision and recall provide more granular information about model performance. To better leverage these metrics, we have now incorporated the Area Under the Precision-Recall Curve (AUPRC)7, which systematically evaluates the trade-off between precision and recall across all possible score thresholds. While precision and recall provide intuitive, threshold-specific insights, the AUPRC offers a comprehensive assessment of model performance and reflects the robustness of a method in ranking true positives above false linkages.\nTo specify, external multi-omic datasets for both PBMC cells (CD4 cells, CD14 cells, and CD16 cells) and Gm12878 cells (rep1, rep2) were used to determine the gene-peak linkages predicted by all three tools (BioCRE, Signac, and ArchR) (Response Figure R4a). Promoter capture Hi-C (PCHiC) datasets for matched cells were used as reference dataset (Response Figure R4a). Standardized genomic scope of predicted linkages across all three tools were restricted to ±500 kb upstream and downstream of gene body. Notably, although the total numbers of predicted linkages are theoretically consistent across methods, minor discrepancies may arise due to each tool’s internal preprocessing steps—such as gene and peak filtering—which may exclude certain elements prior to linkage inference. The filtered geneset and peakset shared by all three tools were used for further linkage prediction.\nFor comparison between the external multi-omic datasets and reference datasets, any predicted linkages overlapping the reference set was considered as positives, while non-overlapping linkages were treated as negatives. Due to the distinct features for different data modules between the predicted linkages and the reference dataset (PCHiC), the overlap is inherently limited, leading to a sparse set of positive linkages. Here, to ensure balanced and robust evaluation, we randomly sampled negative linkages ten times at each positive to negative ratios. Specifically, for the PBMC dataset, we evaluated the AUPRC under varying positive-to-negative ratios: 1:1, 1:2, 1:3, 1:4, and 1:5. For GM12878, we used ratios of 1:1 and 1:2, as the larger number of positive linkages—relative to PBMC—limits the feasible extent of negative down-sampling. Besides, tool-specific metrics to predict the linkages were applied: the combined p-value (as detailed in the Methods, Line 1018-1045) for BioCRE, the p-value for Signac, and the FDR for ArchR.\nAs shown, in consistent with our previous evaluation, we found that BioCRE achieves higher AUPRC scores than Signac and ArchR in most comparison groups (CD4 cells, CD16 cells, and GM12878 cells rep1 and rep2), except for CD14 cells (Response Figure R4b-j). Moreover, BioCRE maintained relatively superior performance across different positive-to-negative linkage ratios (Response Figure R4b-j). These results have been included in the revised Supplementary Figure. 5.\n\nResponse Figure R4. Benchmark of BioCRE, Signac, and ArchR using external datasets. a. The scheme illustrating the prediction and validation of gene-peak linkages using external multiome dataset and referenced promoter-enhancer interaction dataset. b-j. Benchmark of BioCRE, Signac, and ArchR using AUPRC metric. Random sampling for positive-to-negative linkages were performed for 10 replicates at each ratio level, with the summarized bar plot showing at the left and the representative AUPRC curve showing at the right. Comparisons between CD4 cells multi-omic and PCHiC dataset (b), CD14 cells multi-omic and PCHiC dataset (c), CD16 cells multi-omic and PCHiC dataset (d), Gm12878 multi-omic rep1 and PCHiC-1 dataset (e), Gm12878 multi-omic rep1 and PCHiC-2 dataset (f), Gm12878 multi-omic rep1 and PCHiC-3 dataset (g), Gm12878 multi-omic rep2 and PCHiC-1 dataset (h), Gm12878 multi-omic rep2 and PCHiC-2 dataset (i), and Gm12878 multi-omic rep2 and PCHiC-3 dataset (j) were shown. Source data was provided.\n\nAs kindly reminded, here we also checked the performance for BioCRE, Signac, and ArchR for the down-sampled datasets (Response Figure R5). As shown, 20%, 40%, 60%, and 80% of the external multiome dataset were down-sampled for each comparison. Generally, the AUPRC metrics for each tool gradually increases along with the increment of sampled sizes. Importantly, BioCRE exhibits consistent higher AUPRC score than Signac and ArchR, especially for CD4 cells, Gm12878 cells rep1 and rep2 (Response Figure R5a, d, e, f, g, h, i). Meanwhile, for CD14 cells, we found that the BioCRE scores at an intermediate level, performing higher than Signac but lower than ArchR (Response Figure R5b), and the superiority of BioCRE emerges with the increment of sampled data sizes for CD16 cells (Response Figure R5c). These results indicate that BioCRE could achieved stable and relative higher performance in most settings.\n\nResponse Figure R5. Benchmark of BioCRE, Signac, and ArchR through dataset down-sampling. Bar plots showing the summarized AUPRC scores for each group. Four down-sampled groups (20%, 40%, 60%, and 80%) for each comparison were performed. Comparisons between CD4 cells multi-omic and PCHiC dataset (a), CD14 cells multi-omic and PCHiC dataset (b), CD16 cells multi-omic and PCHiC dataset (c), Gm12878 multi-omic rep1 and PCHiC-1 dataset (d), Gm12878 multi-omic rep1 and PCHiC-2 dataset (e), Gm12878 multi-omic rep1 and PCHiC-3 dataset (f), Gm12878 multi-omic rep2 and PCHiC-1 dataset (g), Gm12878 multi-omic rep2 and PCHiC-2 dataset (h), and Gm12878 multi-omic rep2 and PCHiC-3 dataset (i) were shown.\n\nThus, through incorporating various external datasets and systematic measurement of the trade-off between precision and recall with various sampled depths, we observed that BioCRE exhibits a consistent better performance in linkage identification than Signac and ArchR.\nDespite these, we acknowledged more work may be required for a firm benchmark of BioCRE with other tools, such as incorporating more external multiome and reference datasets. We have now discussed about this important issue as\n‘Although standard datasets are not yet available to benchmark the efficacy of BioCRE, preliminary results using external datasets indicate its ability to identify key regulatory elements and their interactions at high precision. The incorporation of more external multiome data and matched gold standard dataset will foster a more systematic evaluation of BioCRE.’ in Line 629-633.\n\nFigure 5.\nI appreciate that the authors included a panel showing the enrichment of cells that are Shh+, Noto+, Foxj1+ (ED F14c). Indeed, it does look enriched in the position of the node. Nonetheless, I do not believe the new plots are sufficient to support their claim that there are two distinct subtypes (Shh+Noto-Foxj1- and Shh+Noto+Foxj1+). From the scRNA-seq data, what proportion of the cells fall into these two categories. Presumably they form two subclusters if the authors are confident in defining the subtypes. Is there a plot they can provide to support this claim?\n\nResponse: Thank you for this thoughtful comment. In the single cell data, we have identified 159 node cells and 65 anterior mesendoderm cells, representing approximately 71% and 29% within the E7.5 axial mesendoderm lineage, respectively. As suggested, we re-analyzed the molecular features for E7.5 axial mesendoderm cells based on RNA and ATAC data, respectively. The resulting UMAP visualization shows that the two subpopulations show clear separation in both transcriptome and epigenome modalities (Response Fig. R6a). Through cell type specific gene expression analyses, we also revealed that the two subtypes exhibit distinct gene expression signatures (Fig. 5e; Response Fig. R6b).\nMoreover, we also checked the presence of known marker gene (Shh) for axial mesendoderm lineage8 and the subtype-specific genes identified in this study. As shown, the axial mesendoderm cells (Shh positive) can be clearly divided into Shh+Noto-Foxj1-Col2a1+Srd5a2+ (Anterior mesendoderm) and Shh+Noto+Foxj1+Col2a1-Srd5a2- (Node) subtypes (Response Fig. R6c). RNAscope experiments further validated the distribution of subtype-specific genes at the E7.5 embryo (Response Fig. R6d, e).\nThese results support our interpretation that there are two transcriptionally and spatially distinct subtypes (anterior mesendoderm and node) within the axial mesendoderm population at the E7.5 embryo. These results have been included in the revised Supplementary Fig. 15e, g, h, i and Line 502-506.\n\nResponse Figure R6. The presence of two subtypes within the axial mesendoderm lineage at the E7.5 embryo. a. UMAP clustering of axial mesendoderm cells at the E7.5 mouse embryo. Cells were colored based on the annotation of node and anterior mesendoderm cells, respectively. b. Heatmap showing the direct comparison of differentially expressed genes for the two axial mesendoderm subtypes (anterior mesendoderm and node). c. Violin plot showing the expression of Shh, Col2a1, Srd5a2, Noto, and Foxj1 in anterior mesendoderm and node cells. d. e. RNAscope validating the expression of axial mesendoderm lineage marker-Shh (d), and co-staining of Col2a1 and Foxj1 (e, top panel), Srd5a2 and Noto (e, bottom panel) at the E7.5 embryos. AME is short for anterior mesendoderm.\n\nFig 5l. I believe the authors mean “specific expression” and not “specific distribution”\n\nResponse: Thank you for pointing this out. We have revised it in the manuscript (Line 540), accordingly.\n\nFigure 6.\nThank you for changing the wording in the manuscript to better match the figure panel. Given the complexity about the dynamics of the system, I believe it will be worthwhile to include a quantitative figure panel better illustrating the temporal dynamics of the regulation. Can the authors please include bar plots or line plots?\n\nResponse: Thank you for your valuable suggestion. To better illustrate the temporal dynamics of the regulation and provide a clearer quantitative assessment, we have added the bar plots as recommended (Response Fig. R7a-e). As shown, the bar plots clearly depict the temporal dynamics of Pou6f1 expression together with its target geneset and target peakset. Specifically, the target peakset for POU6F1 (blue bars) becomes accessible at the expected P, M, EA, and EP regions (highlighted by dashed box) at the E7.0 stage, when both Pou6f1 (red bars) and the target geneset (yellow bars) remain minimally expressed (Response Fig. R7a-c). By the E7.5 stage, both Pou6f1 and its target geneset show abundant expression at the same regions (1EAP, 1A, and 2EA), while with the target peakset show a more broad distribution (Response Fig. R7e).\nFurthermore, we quantitatively assessed the spatial similarities among different modules of POU6F1 TF regulon using Pearson correlation and Structural Similarities index (SSIM) (Response Fig. R7f, g). As shown in the scatter plots, we found that both Pou6f1 expression and its target geneset show limited expression before E7.5 stage (Response Fig. R7f). Notably, only 1EAP, 1A, and 2EA regions of the E7.5 embryo exhibit high similarity and co-expression for both modules (Response Fig. R7f). In contrast, for Pou6f1 expression and its target peakset, although the target peakset get accessible as early as E7.0, Pou6f1 expression becomes prominent only at the E7.5 stage at the 1EAP, 1A, and 2EA regions (Response Fig. R7g).\nTogether, these quantitative results provide additional consistent evidences for the temporal dynamics of the POU6F1 regulon. We have now included in the revised Supplementary Fig. 17a-g.\n\nResponse Figure R7. Quantitative measurement of the temporal dynamics for POU6F1 regulon. a-e. Bar plots showing the relative enrichment of Pou6f1 expression (red bars), its target geneset expression (yellow bars), target peakset accessibility (blue bars) at the E6.5 (a), E6.75 (b), E7.0 (c), E7.25 (d), and E7.5 (e) embryos. f. g. Scatter plots showing the relationship among different modules of POU6F1 regulon (f. Pou6f1 expression and its target geneset expression; g. Pou6f1 expression and its target peak accessibility."}]

Decision letter