Spatial patterning of the epigenome during vertebrate gastrulation

Domain: biology
Published: 2025-12-15
Editorial outcome: major_revision

Abstract
A central question in developmental biology is how embryonic cells acquire and store positional information during pattern formation. In vertebrates, this process begins with the localized activation of signaling systems that mediate axial specification. How these spatial cues are recorded within the regulatory landscape of cells has remained unclear. Here, we report that the chromatin landscape of embryonic cells is spatially patterned during gastrulation. Using spatially resolved genomic analysis in avian embryos, we observed that the epigenome becomes organized in gradients of accessibility and activity along the embryonic axes. These chromatin gradients are established within the loci of developmental genes at the onset of gastrulation and can be used to infer the position of cells in the embryo. Our results show that axial specification involves the spatial organization of the epigenome, linking patterns of chromatin activity to the emergence of the embryonic body plan.

Introduction
Embryonic patterning depends on molecular systems that convey spatial cues to developing cells. This requires the generation of positional information, a set of molecular coordinates that define where cells reside within the embryo1–3. Positional information has been a fundamental concept in developmental biology, offering a framework to explain how progenitor cells are patterned in space. It was outlined by Lewis Wolpert1, based on his work in Hydra regeneration and on classical experiments performed by Hans Driesch, Thomas Hunt Morgan, Hildegard Stumpf, and others4. Wolpert postulated that pattern formation is initiated by the establishment of gradients of signaling molecules along the embryonic axes, known as morphogens5. Cells that are closer to the source of the morphogen receive a higher concentration of the molecule, while cells further away receive a lower concentration6. These gradients provide spatial context, allowing cells to adopt distinct molecular identities depending on their position1,3. This framework anticipated several landmark discoveries, including the identification of the first morphogen in Drosophila7 and the demonstration that extracellular signals can direct the activation of gene regulatory networks that determine regional cell identity8,9.

Much of the positional information required for the regionalization of the early embryo emerges during axial specification10. This refers to the establishment of molecular asymmetries along the anterior-posterior (AP), dorsal-ventral (DV), and medial-lateral (ML) axes in the developing embryo. In Drosophila melanogaster, gradients of maternally derived transcription factors establish the embryonic axes, resulting in a system of spatial coordinates across the syncytial embryo7,11,12. In vertebrates, which develop from a cellularized blastoderm, axial specification is driven primarily by intercellular signaling. WNTs, BMPs, FGFs, and retinoic acid have all been implicated in establishing embryonic asymmetries13–15. Ligands are produced by specific cells and generate concentration gradients via passive or active diffusion, and cells respond to these gradients by activating different sets of genes8. Varying concentrations of signals and patterning genes are thus thought to encode positional information2. However, recent studies have questioned whether morphogen gradients alone are sufficient to account for the precision and stability of early patterning, suggesting that additional regulatory mechanisms may be involved16.

A cellular process that may play a pivotal role in embryonic patterning is the regulation of the chromatin landscape. As early as 1969, long before the characterization of mechanisms of epigenetic regulation, Wolpert proposed that the establishment of positional information would require a coordinate system of spatial genomic responses in the embryo1. Yet, despite decades of research, how the epigenome is remodeled during early pattern formation remains poorly understood. A major obstacle has been the technical challenge of profiling chromatin organization in a spatially resolved manner. To address this, we leveraged the large size and accessibility of the avian embryo to perform spatially resolved genomic analyses during gastrulation. This approach enabled us to map chromatin accessibility and gene expression along the embryonic axes with spatial resolution. Our results demonstrate that axial specification involves the spatial patterning of the epigenome, and that chromatin accessibility changes in accordance with cellular position. We find that a subset of accessible cis-regulatory elements displays regionalized activity that can be used to infer where individual cells reside within the embryo and to reconstruct spatial patterns of gene expression. These findings suggest that early embryonic cells capture positional cues within their chromatin landscape through deployment of specialized cis-regulatory elements.

Results
Reconstruction of gene expression patterns in the early avian embryo
To map the transcriptional signatures that underlie axial specification, we conducted a genome-wide reconstruction of gene expression patterns in the early chick embryo. Like its human counterpart, the avian embryo develops as a flat, trilaminar disc, but its larger size makes it particularly well-suited for spatially resolved transcriptomic analysis. We performed Tomography-Seq (tomo-seq)17 on embryos at the end of gastrulation (Hamburger-Hamilton stage 5, HH5)18, which were cryosectioned along the AP and ML axes (Fig. 1a). Individual 100 um tissue sections were then subjected to 3’ RNA-seq (Fig. 1b). Tomo-seq allowed for genome-wide quantification of gene expression levels along embryonic axes, with normalized mRNA levels successfully recapitulating the axial domains of genes like OTX2, MEOX2 and TBXT (Fig. 1c). Integration of the AP and ML tomo-seq datasets allowed us to reconstruct gene expression patterns in the avian embryo (Fig. 1d) and to simulate in situ hybridization experiments in silico for hundreds of genes (Supplementary Fig. 1). Inferred gene expression patterns were consistent with those of known axial-restricted genes (Fig. 1c–e, Supplementary Fig. 2a), including OTX2, CHRD and WNT5A. Identification of genes with similar expression patterns also allowed us to delineate gene signatures for morphological landmarks like the anterior neural plate, primitive streak, and the notochord (Supplementary Data 1).

Next, we used the tomo-seq datasets to identify genes with spatially patterned expression along the AP axis. Transcriptomes of 100 μm tissue sections from two biological replicates showed strong concordance between individuals (Supplementary Fig. 2d, e). Hierarchical clustering of gene expression profiles identified distinct clusters of gene expression profiles along the AP axis (Supplementary Fig. 2b). We detected expression of approximately 13,500 genes (>50 counts across all slices), of which 1000 displayed a statistically significant association between gene expression and position along the AP axis. Of these, 706 genes were found to be highly enriched at the anterior region of the embryo (Supplementary Data 2), including previously characterized anterior factors such as OTX2, SIX3, and CYP26A (Fig. 1f, Supplementary Data 2)19–21. In contrast, only 294 genes were found to be posteriorly enriched. This cohort of genes contained HOX family members, and other posteriorizing factors like CDX2 and WNT5A (Fig. 1f, Supplementary Data 2)22,23. Gene ontology (GO) analysis of anterior-enriched gene sets showed an enrichment of terms associated with neurogenesis, head development, and sensory structure morphogenesis, while those associated with posterior-enriched genes included mesenchymal tissue formation, including heart and skeletal morphogenesis (Supplementary Fig. 2c). Averaging in silico expression patterns of all anterior genes showed that their expression was the strongest around the Hensen’s node, whereas posterior genes were enriched around the primitive streak (Fig. 1g). Together, these results allowed for high-throughput reconstruction of gene expression patterns along the AP and ML axes and led to the identification of AP gene signatures that are established during avian gastrulation.

Spatial patterning of the epigenome during axial specification
Next, we examined the organization of the epigenome along the AP axis of the gastrula. We divided HH5 avian embryos into seven transverse sections along the AP axis and subjected them to ATAC-seq analysis24 (Fig. 2a). Quality control analysis of our ATAC-seq datasets showed an expected nucleosomal distribution of ATAC-seq fragments, as well as enrichment at transcription start sites (Supplementary Fig. 3a). This tomoATAC-seq strategy led to the identification of a set of cis-regulatory elements (CREs) that display a graded pattern of chromatin accessibility along the AP axis (Fig. 2b). These genomic regions, which we termed regional cis-regulatory elements (rCREs), were particularly abundant in the loci of genes with regionalized expression patterns (Supplementary Fig. 3b). We identified a significant association between genes enriched for rCREs in their loci and those with anterior/posterior enriched expression patterns (P=3.635e-06, chi-square test). For instance, the OTX2 locus contained rCREs that were highly accessible at the anterior-most sections but displayed a gradually weaker ATAC signal in more posterior sections (Fig. 2c). The converse was observed in the locus of WNT5A, which contained many posterior rCREs (Fig. 2d). Quantification of normalized chromatin accessibility confirmed that, rather than recapitulating transcriptional domains, rCREs formed stereotypical gradients of accessibility along the AP axis (Fig. 2e, f, Supplementary Figs. 4, 5). Chromatin gradients were not uniform within a given locus, as rCREs were often interspersed by elements that were not patterned or that were patterned in the opposing direction along the AP axis (arrows, Fig. 2c, d). Hierarchical clustering of the ATAC-seq datasets revealed varying patterns of chromatin accessibility along the AP axis (Supplementary Fig. 3c). We next compiled anterior and posterior sets of rCREs to assemble AP epigenomic signatures. Differential accessibility analysis identified two sets of 10913 and 12532 elements that were significantly enriched at the anterior and posterior poles of the embryo, respectively, which we deemed anterior and posterior rCREs (Fig. 2b, Supplementary Fig. 3d; Supplementary Data 3). Gene onthology (GO) analysis of the genes associated with anterior or posterior rCREs revealed distinct functional enrichments corresponding to developmental programs operating in each region (Supplementary Fig. 3e). Anterior rCRE-associated genes were enriched for biological processes related to BMP signaling, axis specification, and sensory organ development. In contrast, posterior rCRE-associated genes were linked to processes such as cell motility, embryonic patterning, and tissue morphogenesis.

To identify the main regulators of AP identity in the embryo, we conducted transcription factor motif enrichment analysis with HOMER25 on anterior and posterior rCRE sets. Peaks active at the anterior end of the embryo were enriched for motifs for OTX2, GSC, SOX3, and FOXA2, while those active at the posterior end of the embryo contained HOX, CDX2, and TBXT motifs (Supplementary Fig. 3f). Analysis of transcription factor motif activity across the seven ATAC-seq sections via chromVAR26 confirmed high variability of anterior and posterior transcription factor binding motifs along the AP axis (Fig. 2g). Furthermore, we found a strong correlation between chromatin accessibility and deposition of H3K27AC at anterior and posterior rCREs, confirming that these genomic regions are not only open, but also display polarized chromatin activation in the embryo (Fig. 2b, Supplementary Fig. 3g). Similar to chromatin accessibility, H3K27AC deposition in rCREs was also graded along the AP axis (Fig. 2b, e, f, Supplementary Figs. 4, 5). Together, these results reveal that the chromatin landscape of the gastrula is spatially patterned across the avian gastrula. The two-dimensional expression of OTX2 is shown for reference (e). e In silico (left) reconstruction of the OTX2, CHRD, MEOX1, and TBXT 2D gene expression patterns, along with the corresponding in situ hybridization images (right). In silico projections recapitulate the gene expression profiles as seen via in situ hybridization. Scale bar: 400 μm. f Heatmap showing 706 anterior- and 294 posterior-enriched genes (p-value < 0.01, log fold change > 1/<−1) expressed in the avian gastrula, as determined via TradeSeq’s startVsEndTest function. g Spatial averages of anterior and posterior genes in the gastrula. Source data are provided as a Source Data file. AP anteroposterior, ML mediolateral.

Gradients of regional CREs are shaped by signaling pathways
To determine when rCRE gradients first appear during development, we analyzed chromatin patterning during blastula and gastrula stages (Fig. 3a). HH1 (specifically EGK XIII27), HH2, and HH3 embryos were divided into five sections along the AP axis and subjected to CUT&RUN for H3K27AC. When visualizing trends of chromatin activity at anterior and posterior rCREs, we observed little to no chromatin polarization at the early (EGK XIII) and late (HH2) blastula stages. Graded activation of rCREs became conspicuous at HH3, coinciding with the onset of gastrulation (Fig. 3a, Supplementary Fig. 6a). To verify how rCRE polarization correlates with transcription, we also performed RNA-seq along the AP axis in these early stages. The results showed that expression of the tomo-seq AP gene sets (Fig. 1) became regionalized following the establishment of rCRE gradients (Supplementary Fig. 6a). To confirm that embryos were positioned correctly prior to sectioning, we examined the NODAL gene locus. At all three stages assayed, NODAL transcripts were found to be enriched in the posterior-most sections (Supplementary Fig. 6b), consistent with its expression in Kohler’s sickle and the early primitive streak28,29. To further explore the relationship between early chromatin patterning and transcriptional activation, we examined individual loci of patterning genes that were transcriptionally inactive at the stages of rCRE activation. We observed that in such genes, including DMBX1 and HOXA11, rCRE gradients were established prior to regionalized transcription (Supplementary Fig. 7a, b). These observations indicate that chromatin accessibility gradients can form ahead of gene expression, suggesting that rCRE activation may help prime regulatory landscapes for subsequent transcriptional responses during axis formation.

We next examined the mechanisms that shape chromatin gradients during gastrulation. A well-known positional cue involved in the establishment of the anteroposterior axis is WNT signaling14. This highly conserved signaling system acts as a posteriorizing agent across metazoans30. We thus tested the hypothesis that rCREs are directly regulated by nuclear effectors of WNT signaling. Several lines of evidence from the analysis of our tomo-seq datasets point to a role for WNT signaling in the establishment of rCREs. First, a number of WNT ligands and WNT pathway-associated genes are restricted in expression along the AP axis (Supplementary Data 2). Moreover, motif enrichment analysis of posterior rCREs revealed robust enrichment of motifs for nuclear effectors of WNT signaling, including LEF1 and TCF-like factors (Fig. 2g, Supplementary Fig. 3f). To investigate the relationship between rCREs and WNT signaling, we performed CUT&RUN for LEF1/CTNNB1 in three sections along the AP axis of HH5 embryos (Fig. 3b). Analysis of our CUT&RUN datasets revealed 3729 peaks, representing genomic regions bound by LEF1/CTNNB1. Motif enrichment analysis of LEF1/CTNNB1 peaks showed TCF/LEF as the highest-ranked motifs, confirming our ability to isolate regions bound by these factors (Fig. 3c). Furthermore, consistent with the posteriorizing role of WNTs, we found that LEF1/CTNNB1 were predominantly associated with posterior rCREs (Fig. 3d), and that CUT&RUN signal displayed an anterior-to-posterior gradient (Fig. 3d, e). Analysis of earlier embryos (HH3) revealed that the majority of LEF1/CTNNB1 peaks are already present in the early gastrula (Supplementary Fig. 8a, b). Moreover, while we observed greater enrichment of LEF1/CTNNB1 at anterior rCREs at HH3 than HH5, the CUT&RUN signal was consistently stronger at the posterior end of the embryo at both stages (Fig. 3d, Supplementary Fig. 8c). Together, these results indicate that posterior rCREs are regulated by components of the WNT signaling pathway.

To test if WNT signaling is consequential for chromatin patterning, we next assessed the effect of modulation of WNT signaling on rCREs. We isolated cells of the anterior epiblast at stage HH2 (prior to rCRE patterning) and cultured the cells for 24 h under WNT gain-of-function (GOF) and loss-of-function (LOF) conditions (Fig. 3f). The media was supplemented with either a WNT pathway agonist, CHIR9902131, a WNT antagonist, XAV93932, or DMSO as a control, and cells were subjected to ATAC-seq following treatment. Analysis of anterior- and posterior-enriched rCREs revealed a significant change in accessibility upon WNT modulation. We observed a decrease in the accessibility of anterior rCREs upon treatment with CHIR99021, and a concomitant increase in accessibility following XAV939 treatment (Fig. 3g, Supplementary Fig. 8d). The opposite trend was observed for posterior rCREs, which displayed an increase in accessibility following CHIR treatment and a decrease in accessibility upon treatment with XAV939. These results demonstrate that posteriorizing signals enhance chromatin accessibility at posterior rCREs while suppressing accessibility at anterior sites.

We next examined how peaks bound by LEF1/CTNN1B responded to manipulations in WNT signaling. We projected LEF1/CTNN1B-bound peaks onto our WNT modulation epigenomic dataset (Fig. 3h, i). As expected, genomic regions associated with LEF1/CTNN1B binding were positively regulated by WNT signaling (i.e., their accessibility is promoted by CHIR99021 and repressed by XAV939). The integration of these datasets allowed us to compile a list of genomic targets of canonical WNT signaling that respond to the manipulation of the pathway (Supplementary Data 4). Finally, we examined the spatial patterns of accessibility in WNT-responsive genomic regions to examine if they display features of rCREs. Projection of the accessibility of WNT-responsive peaks (Fig. 3j, k) in our tomoATAC-seq datasets showed that they display AP polarization and graded accessibility along the AP axis (Fig. 3l). As expected, peaks activated by WNTs formed posterior to anterior gradients, while peaks inhibited by the signaling system are anteriorly polarized (Fig. 3l). These results show that signaling systems modulate AP gradients of chromatin accessibility, linking signaling pathway activity to the spatial organization of the epigenome.

Chromatin gradients enable inference of cellular position along the embryonic axis
The finding that rCREs are established during gastrulation led us to postulate that these elements capture positional cues in embryonic cells. To test this, we first examined if rCRE signatures (Supplementary Data 3) can be used to infer the position of embryonic cells along the AP axis. We thus performed single cell multiome analysis33 of the avian gastrula using the 10X genomics platform to measure both chromatin accessibility and gene expression in individual cells (Fig. 4a). This analysis allowed us to estimate cellular AP position based on rCRE accessibility, and to subsequently test the accuracy of this prediction based on the gene expression signatures from the tomo-seq analysis (Fig. 1). Upon filtering for high quality cells (n=7750) and performing dimensionality reduction, we identified 16 distinct clusters within the multiome dataset (Fig. 4b). Clusters were distributed similarly within the UMAP projection regardless of which dataset was used for dimensionality reduction (RNA-seq, ATAC-seq, or both) (Supplementary Fig. 9).

Cellular identities were assigned to individual clusters using a combination of cluster-specific marker gene expression from scRNA-seq, along with motif enrichment analysis from scATAC-seq (Supplementary Fig. 10A–C). Cluster annotation showed separation of cell identity primarily based on their germ layer of origin. Of the 16 identified clusters, six clusters were mesodermal, four clusters were ectodermal, and one cluster was endodermal (Fig. 4b). Extraembryonic tissues, including blood islands, extraembryonic endoderm, and the area opaca were also identified in the dataset. After cluster assignment, we attempted to infer the position of individual cells using our anterior and posterior epigenomic signatures. First, we examined whether anterior and posterior rCREs were associated with regionalized gene expression. We employed peak-gene co-accessibility analysis34 in the multiome dataset to assign rCREs to their cognate target genes (Fig. 4c) and visualized the expression profiles of these genes along the AP axis. As expected, genes found to be co-accessible with AP rCREs displayed spatial enrichment of anterior or posterior tomo-seq gene signatures (Fig. 4d). Consistent with these findings, GO enrichment analysis of genes linked to anterior and posterior rCREs revealed region-specific developmental signatures (Supplementary Fig. 10D), such as sensory organ and brain development in anterior regions and somitogenesis and mesoderm formation in posterior regions.

Next, we attempted to order cells along the AP axis based on epigenomic information. We utilized a quadratic programming approach35,36 to calculate the fractional identity of each cell based on its similarity to the anterior vs. posterior rCRE signatures (Fig. 4e). This strategy allowed us to assign the approximate position of each cell based on rCRE score (Fig. 4f). To test the accuracy of rCRE-based cell ordering, we projected transcriptomic information from the multiome onto the ordered cells. We first assessed the expression of the anterior- and posterior-enriched gene sets identified in the tomo-seq analysis (Fig. 1f, Supplementary Data 2). We found that these gene sets were enriched in the anterior and posterior positioned cells, validating our cell ordering analysis (Fig. 4g). Reconstruction of expression patterns of individual AP genes, including OTX2, WNT5A, and TCF7, were consistent with their expression domains (Fig. 4h). This was also true for genes expressed in intermediate AP positions within the embryo, such as CHRD and SHH, indicating that rCRE gradients contain positional information along the entire embryonic axis. Cell ordering was also able to reconstruct the early HOX code in the gastrula (Fig. 4h). Expression patterns were consistent with chromVAR analysis, which identified enrichment of motifs for regulators at their estimated expression domains (Fig. 4i). Together, these results demonstrate that rCREs can be used to infer the axial position of embryonic cells in silico and support the conclusion that spatial information can be detected in the chromatin landscape.

Regional CREs provide a spatial register of cells within germ layers
Our multiome analysis demonstrated that the avian gastrula is composed of 16 different cell types (Fig. 4b). We next tested whether we could use the combination of rCRE signatures and quadratic programming to define how these cell types are organized along the AP axis. We first examined the cluster identity of cells enriched for anterior vs. posterior rCREs. We found that the anterior cells are largely ectodermal and endodermal, while the posterior cells consist mostly of mesoderm (Fig. 5a). Indeed, projection of the anterior rCREs on the single-cell multiome UMAP shows robust signal in the ectodermal and endodermal clusters, whereas the posterior signal is strongest in the mesodermal cluster (Fig. 5b). To examine AP composition of the embryo in more detail, we generated 50-cell bins ordered from the anterior to the posterior poles of the embryo. We then assigned the identity of the cells in each bin to determine the regional abundance of distinct cell types (Fig. 5c). This confirmed that a larger percentage of ectodermal and endodermal cells were positioned in the anterior half of the embryo, while mesodermal clusters were posteriorly enriched.

To verify the accuracy of this inference, we performed immunostaining in sagittal sections of the avian gastrula with antibodies for markers for the three germ layers (Fig. 5d)37–39. Consistent with the cellular composition analysis, endodermal (SOX17+) and ectodermal (SOX2+) progenitors were enriched in the anterior region of the embryo, whereas the posterior gastrula is predominantly mesodermal (TBXT+). Closer examination revealed that the anterior accumulation of SOX17+ endoderm reflects the presence of the prechordal plate, where the endoderm forms a multi-layered structure. In contrast, the posterior region lacks definitive endodermal cells at this stage, with the underlying layer being largely composed of hypoblast cells (Fig. 5e). Furthermore, we also observed that cells assigned to distinct embryonic structures, including the primitive streak and the notochord, occupied the expected position along the AP axis (Fig. 5c).

The asymmetric distribution of cells along the AP axis led us to test if rCREs could be used to infer the spatial distribution of cells within individual germ layers. We subdivided our multiome dataset into ectodermal, mesodermal, and endodermal clusters (2481, 2936, and 440 cells, respectively) (Supplementary Fig. 11a–c). Projection of the anterior- and posterior-enriched rCREs onto germ layer-specific subclusters revealed distinct subpopulations of cells enriched for axial-restricted peaks in each of the three germ layers (Supplementary Fig. 11a–c). Thus, although unsupervised clustering largely groups cells based on cellular origin, a closer look at intra-cluster organization suggests additional separation based on positional identity. To test this, we examined the axial identity of cells from neural ectoderm, which undergoes early regionalization along the AP axis during the formation of the central nervous system40. We employed rCRE-based cell ordering in the neural ectoderm subcluster (Fig. 5f, 708 cells), and assessed the expression of the forebrain, midbrain, and hindbrain markers (OTX2, GBX2, and HOXA2, respectively) (Fig. 5f)41,42. Ordered cells displayed the expected expression profiles of these marker genes along the AP axis (Fig. 5g).

Next, to determine the extent to which axial-restricted genes are shared across cell types, we conducted trajectory-based differential expression analysis in the ectodermal and mesodermal subclusters (2481 cells and 2936 cells, respectively). This analysis revealed a set of 64 genes significantly associated with AP-restricted expression in both germ layers (Supplementary Fig. 11d, e). Together, these results support the existence of a regulatory code associated with AP patterning that is independent of cellular identity.

Conservation of rCRE signatures across vertebrates
Spatial epigenomic analysis of the avian gastrula allowed for the identification of anterior and posterior rCREs, which provide a molecular signature of cellular position in the embryo. To determine if this mechanism is present on other vertebrate embryos, we analyzed multiomic data from the E7.5 mouse embryos43. Liftover44 of anterior and posterior rCREs identified 1560 anterior and 1724 posterior rCREs conserved in mouse (Fig. 5h). We then used quadratic programming to assign the relative position of individual cells from the E7.5 mouse dataset along the AP axis (Fig. 5i). Upon ordering cells based on their chromatin accessibility at conserved rCREs, we assessed the expression patterns of anterior- and posterior-restricted genes as 200-cell bins along the AP axis. This analysis revealed spatial restriction of anterior genes OTX2 and POU3F1, as well as posterior genes TBXT and WNT3 in the developing mouse gastrula (Fig. 5j)45. Therefore, even a small subset (~15%) of rCREs can be used to estimate anterior/posterior gene expression patterns across species. These findings suggest that regional chromatin organization associated with positional cues is conserved among amniotes and may represent a general strategy for spatial patterning during early development.

Methods
Tissue collection
Fertilized Leghorn White chicken eggs were obtained from the Department of Animal Science, University of Connecticut. Eggs were incubated at 37 °C and 50% humidity until embryos reached the desired developmental stage. Embryos were collected in Ringer’s solution, staged according to the Eyal-Giladi and Kochav27 or Hamburger and Hamilton18 staging systems, and prepared for genomic protocols as described below. EGK XIII embryos were collected before primitive streak formation and staged by assessing epiblast expansion and Koller’s sickle position. Embryos were examined and dissected under a ZEISS Stemi 508 stereo microscope. The anterior-posterior and lateral axes were measured, and embryos of similar size were grouped for processing.

Tomo-seq sample and library preparation.
HH5 embryos were collected, staged, and cleaned in RNase-free solution. Similar-sized embryos were then transferred to freezing molds, positioned depending on the axis to be sectioned, and embedded in frozen section medium (Thermo Scientific™ Richard-Allan Scientific Neg 50™, #6502). Before freezing, each embryo had its Hensen’s node, anterior-posterior and lateral ends marked with dye to serve as landmarks for the sectioning process. In an RNase-free environment, frozen samples were trimmed to avoid extra-embryonic tissues and frozen medium excess, and cryosectioned into 100 μm sections. Each anterior-posterior or medial-lateral section was then collected for total RNA extraction with the RNAqueous-Micro Total RNA Isolation Kit (Thermo Fisher Scientific, #AM1931) according to the kit’s protocol and processed for 3’ RNA-seq library preparation. Libraries were prepared using the QuantSeq 3’ mRNA-seq Library Prep Kit (FWD) for Illumina (Lexogen, #K01596), following the manufacturer’s recommendations for low RNA input samples. Final libraries were quantified, pooled, and sequenced with single-end 75bp reads on an Illumina NextSeq500 instrument.

Tomo-seq data analysis
Trimming of adapter sequences and poly A tails was performed using BBDuk (version 38.93) using the following parameters: ref=polyA.fa.gz, truseq_rna.fa.gz k=13 ktrim=r useshortkmers=t mink=5 qtrim=r trimq=10 minlength=20. Next, STAR (version 2.7.3a)57 was used to align trimmed reads to the galGal6 chicken reference genome using the following arguments: --outFilterType BySJout --outFilterMultimapNmax 20 --alignSJoverhangMin 8 --alignSJDBoverhangMin 1 --outFilterMismatchNmax 999 --outFilterMismatchNoverLmax 0.1 --alignIntronMin 20 --alignIntronMax 1000000 --alignMatesGapMax 1000000 --outSAMattributes NH HI NM MD --outSAMtype BAM SortedByCoordinate. Bam alignment files were indexed using samtools (version 1.10)58 and reads were counted with htseq-count (version 0.13.5)59 using a GTF file containing chicken gene annotations and using the following parameters: -m intersection-nonempty -s yes -f bam -r pos. AP section 14 from replicate 1 was omitted from downstream analyses because it received very few counts and was deemed a failed library. Differential gene expression analysis was performed on two embryo replicates using TradeSeq (version 1.8)60. Anterior- and posterior-enriched genes were identified using TradeSeq’s startVsEndTest function, and genes with a p-value < 0.01 and a log fold change > 1 or < −1 were considered significantly enriched. For plotting of AP and 2D gene expression profiles, raw counts were normalized to the total number of spike-in reads, and the moving average of the normalized counts was plotted along each axis. Gene ontology analysis was performed using PANTHER’s61,62 Fisher’s Exact test with Bonferroni correction.

In situ hybridization
For in situ hybridization, embryos were fixed in phosphate buffered saline (PBS) containing 4% paraformaldehyde (PFA) for 2 h at room temperature or overnight at 4 °C. Following fixation, embryos were dissected, washed with PBS containing 0.1% Tween (PBST), dehydrated, and stored in methanol at −20 °C. Whole-mount in situ hybridization was performed as previously described63. Embryos were imaged using a ZEISS Stemi 508 stereo microscope.

HCR RNA-FISH
Hybridization chain reaction RNA fluorescent in situ hybridization (Molecular Instruments HCR v3.0) was performed in HH5 (19-22h) whole mount chick embryos64. Briefly, embryos were collected and fixed in 4% paraformaldehyde for 2 h at room temperature. Embryos were dissected to remove the vitelline membrane, washed with PBST on ice and then stored in methanol at −20 °C until use. Next, embryos were rehydrated with a series of graded Methanol/PBST washes on ice. For the detection step, embryos were pre-hybridized with probe hybridization buffer for 30 min at 37 °C, then incubated with 4 pmol of each probe set overnight. The excess of probe was removed by washing the embryos with probe wash buffer at 37 °C. The embryos were then pre-amplified with amplification buffer and incubated with hairpin solution overnight in the dark at room temperature. Excess hairpins were removed by washing the embryos with 5X SSCT at room temperature. Embryos were incubated with DAPI for 1 h followed by sample mounting for microscopy. Embryos were imaged at the Nikon Imaging Center at Harvard Medical School using a Nikon AXR point scanning confocal. The target mRNA probe sets, amplifiers and buffers were purchased from Molecular Instruments (https://www.molecularinstruments.com/).

Immunohistochemistry
For whole-mount immunohistochemistry, embryos were collected and fixed in 4% PFA-PB for 20 min at room temperature. Post fixation, embryos were dissected from the filter paper and washed in TBS containing 0.1% Triton and 1% DMSO (TBTD). Embryos were blocked at room temperature for 2 h in TBTD supplemented with 10% donkey serum and incubated in anti-SOX2 (Abcam, #ab97959) and anti-SOX17 (R&D, #AF1924) primary antibodies diluted in blocking solution, overnight at 4 °C. Following the primary antibody incubation, embryos were washed, blocked for 30 min at room temperature, and stained with appropriate secondary antibodies and anti-Brachyury Alexa Fluor 488 Conjugated (R&D, #IC2085G) for 2 h at room temperature. Secondary antibodies used included donkey anti-rabbit/goat IgG conjugated with Alexa Fluor 594/647 (Molecular Probes). Following the secondary antibody step, the embryos were washed, stained with DAPI and post-fixed with 4% PFA for 1 h, prior to imaging. Whole-mount images were collected using an upright ZEISS Axio Imager fluorescent microscope.

TomoATAC-seq
HH5 chicken embryos were collected in Ringer’s solution and stage- and size-matched embryos were dissected out of any extraembryonic tissue. Then, using Hensen’s node as a landmark, seven sections were obtained along the anteroposterior axis, three anterior to the node and four posterior to it. For EGK XIII, HH2, and HH3 embryos, staging was based on established morphological criteria, and embryos were similarly matched for size and developmental progression. These earlier embryos were also cryosectioned into five evenly spaced sections along the embryonic axis. Five position-matching sections were pooled and dissociated in Accumax (Innovative Cell Technologies, #AM193). Cell suspensions were washed and processed using the OMNI-ATAC-Seq protocol65. Briefly, dissociated cells were washed in resuspension buffer containing 10 mM Tris-HCl pH 7.4, 10 mM NaCl, and 3 mM MgCl2. After washing, cells were resuspended in lysis buffer composed of the resuspension buffer plus 0.1% NP-40 (Sigma, I8896), 0.1% Tween-20 (Sigma, P1379), and 0.01% Digitonin (Millipore, 300410). The samples were incubated on ice for 3 min, followed by washing with resuspension buffer containing only 0.1% Tween-20. Nuclei were then resuspended in 50 μL of transposition mixture containing 10 mM Tris-HCl, pH 7.4, 5 mM MgCl2, 10% Dimethyl Formamide (Sigma, D4551), 1X PBS, 0.1% Tween-20, 0.01% Digitonin, and 5 μL Tn5 transposase (Illumina, 20034210). Transposition was performed for 1 h at 37 °C in a thermal mixer. The pre-amplification DNA was purified using the Qiagen Minelute Reaction Cleanup Kit (Qiagen, 28204). Library amplification PCR was performed using the Q5 High-Fidelity 2X Master Mix (NEB, M0492S). Library size selection and cleanup was carried out using Ampure XP beads (Beckman Coulter, A63881). Final libraries were quantified, pooled, and sequenced with paired-end 37bp reads on an Illumina NextSeq500 instrument.

TomoATAC-seq data analysis
Trimming of Nextera adapter sequences was performed using cutadapt (version 3.4)66 using the following parameters: -a CTGTCTCTTATACACATCT -A AGATGTGTATAAGAGACAG --minimum-length=25. Trimmed reads were aligned to the galGal6 chicken reference genome using bowtie2 (version 2.4.4)67 with the following arguments: --local --very-sensitive-local --no-unal --no-mixed --no-discordant -X 2000. Bam alignment files were indexed using samtools (version 1.10), followed by marking of duplicate reads using Picard’s (version 2.26.2) MarkDuplicates function. Duplicate reads were then removed using samtools view with the following arguments: -F 1804 -f 2. Bigwigs were generated from Bam files using DeepTools’ bamCoverage function (version 3.5.1)68 with the following arguments: --binSize 5 --normalizeUsing RPKM --extendReads. All downstream analyses were conducted using peaks called using cluster-specific peakcalling on the HH5 multiome dataset (316,412 peaks, see 10X multiome data analysis). By calling peaks on individual cell clusters, this increased the sensitivity of peak calling and ensured rare, cell-type-specific regulatory elements were not omitted from the analysis. TomoATAC-seq reads in these 250,772 peaks were counted using DiffBind (version 3.4.11)69. Peaks with anterior- and posterior-enriched accessibility were identified using TradeSeq’s (version 1.16.0) startVsEndTest function and peaks with a p-value < 0.01 and a log fold change > 1 or < −1 were considered significantly enriched. For plotting of accessibility profiles, raw counts were library size normalized and the moving average of the normalized counts was plotted along the AP axis. Gene enrichment analysis was performed by assigning each peak to the nearest annotated transcription start site using GenomicRange’s70 (version 1.54.1) nearest() function. Next, a binomial test was conducted using the total number of cis-regulatory elements assigned to a gene as the number of trials and the proportion of all rCREs over the total number of cis-regulatory elements as the probability of success. Genes returning a p-value of </= 0.05 were considered significantly enriched for rCREs. Motif enrichment analysis was performed by calculating motif deviations using ChromVar (version 1.16)26, as well as using HOMER’s25 findMotifsGenome.pl function. Functional enrichment analysis of anterior and posterior rCREs was conducted using GREAT71 (version 2.4.0).

Anterior epiblast ATAC-seq
As described previously, chick embryos were obtained from fertile eggs, collected, and staged in Ringer’s solution. HH2 embryos were dissected off any extra-embryonic tissue and had the anterior hypoblast removed. Then, the anterior, naive, epiblasts from three embryos were microdissected, pooled, and transferred to 96-well plates previously coated with fibronectin and containing 100 uL of 10% FBS DMEM culture media containing the WNT agonist CHIR99021 (Tocris, #4423), or antagonist XAV939 (Sigma, #X3004-5MG) at the final concentrations of 10 uM (the same volume of vehicle, DMSO, was used for control explants). The explants were incubated for 24 h at 37 °C in a CO2 incubator. After this period, explants were washed in PBS and collected with Accumax. Cell suspensions were washed and processed using the OMNI-ATAC-Seq protocol65 as described above.

Epiblast ATAC-seq data analysis
Initial processing of ATAC-seq libraries, including trimming, alignment, removal of duplicate reads, peak calling, and generation of bigwigs, was performed in the same manner as the tomoATAC-seq experiments. After initial processing, ATAC-seq reads were counted at 250,772 peaks obtained from the analysis of the 10X multiome data DiffBind (version 3.4.11). Pairwise differential accessibility analysis between the DMSO vs CHIR99021 and DMSO vs XAV939 conditions was performed via DESeq272 (version 1.42.1) using three biological replicates. Wnt-responsive (increased) peaks were defined as those enriched in the CHIR99021 condition (log2 fold change > 0.5) and depleted in the XAV939 condition (log2 fold change < −0.5). Wnt-responsive (decreased) peaks were defined as those enriched in the XAV939 condition (log2 fold change > 0.5) and depleted in the CHIR99021 condition (log2 fold change < −0.5).

CUT&RUN
As described above for the TomoATAC-seq experiments, the embryos processed for CUT&RUN were first collected and stage- and size-matched. For each experiment, position-matched sections from 3 embryos were pooled together. H3K27AC CUT&RUN experiments were performed with EGK XIII, HH2, HH3, and HH5 embryos sectioned into 5, 5, 6 and 7 sections, respectively. Embryonic structures like Hensen’s node and primitive streak were used as landmarks. The sectioning was performed in a way to preserve the amount of tissue per fragment. For the CTNNB1/LEF1 CUT&RUN, embryos from three biological replicates were dissected into 3 anterior-posterior sections. The rostral segment was cut right after the regressing Hensen’s node and the remaining embryo was cut in half, dividing the middle and caudal segments. CUT&RUN experiments were conducted as previously described73. First, the tissue was dissociated into a single cell suspension using Accumax at room temperature for 20 min under mild agitation. Next, cells were bound to BioMag Plus Concavalin A magnetic beads (Bangs Laboratories, BP531) and incubated with rabbit anti-H3K27AC (AbCam, #ab177178), or rabbit anti-CTNNB1 (AbCam, #ab32572) plus anti-LEF1 (AbCam, #ab137872) antibodies overnight at 4 °C, rocking. All antibodies were used at a 1:50 concentration. After washing away the unbound antibody, protein A-MNase was added to a final concentration of 700 ng/mL and incubated on the tube rotator for 1 h at 4 °C. Cells were chilled down to 0 °C, and CaCl2 was added to a final concentration of 2 mM to activate the MNase enzyme. MNase digestion was conducted for 45 min and terminated by adding the 2X STOP buffer containing heterologous Saccharomyces cerevisiae spike-in DNA at a concentration of 2 pg/mL. The protein-DNA complexes were released by centrifugation and digested with proteinase K for 10 min at 70 °C. DNA fragments were extracted using phenol-chloroform and ethanol precipitation. CUT&RUN libraries were prepared using the NEBNext Ultra II DNA Library Prep Kit (New England Biolabs, #E7645) following the manufacturer's protocol. Fragment analysis was performed with Agilent 4200 TapeStation to perform quality control for the libraries. Equimolar concentrations of the libraries were pooled and sequenced with paired-end 37bp reads on an Illumina NextSeq500 instrument.

CUT&RUN data analysis
Trimming of TruSeq adapter sequences was performed using cutadapt (version 3.4) using the following parameters: -a AGATCGGAAGAGCACGTCTGAACTCCAGTCA -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT --minimum-length=25. Trimmed reads were aligned to the galGal6 chicken reference genome using bowtie2 (version 2.4.4) with the following arguments: --local --very-sensitive-local --no-unal --no-mixed --no-discordant -I 10 -X 1000. Bam alignment files were indexed using samtools (version 1.10), followed by marking of duplicate reads using Picard’s (version 2.26.2) MarkDuplicates function. Duplicate reads were then removed using samtools view with the following arguments: -F 1804 -f 2. Peaks were called with macs2 (version 2.2.7.1) using the parameters: -f BAMPE -g 1e9 -q 0.05 --call-summits. Bigwigs were generated from Bam files using DeepTools’ bamCoverage function (version 3.5.1) with the following arguments: --binSize 5 --normalizeUsing RPKM --extendReads. Reads in peaks were counted using DiffBind (version 3.4.11). For plotting of chromatin activity profiles, raw counts were library size normalized and the moving average of the normalized counts were plotted along the AP axis. Motif enrichment analysis and identification of the top enriched transcription factor binding motifs were performed using HOMER’s findMotifsGenome.pl function. For LEF1/CTNNB1 CUT&RUN analysis, consensus peaksets from macs2-called peaks were obtained using DiffBind’s dba.count() function with minOverlap=3. Average bigwigs were generated from three replicates using DeepTools’ (version 3.5.1) bigwigAverage function. Plotting of LEF1/CTNNB1 binding scores at anterior-/posterior-enriched peaks was accomplished using DeepTools’ computeMatrix and plotHeatmap functions, using averaged bigwigs as input. LEF1/CTNNB1 peaks were mapped onto ATAC-seq peaks using the GenomicRanges findOverlaps() function.

Tissue isolation and 10X multiome library preparation
The nuclei of three whole chick embryos at stage HH5 were isolated following the Chromium Next GEM Single Cell Multiome ATAC + Gene Expression (GEX) protocol. Briefly, the embryos were dissected and dissociated using Accumax for 15-20 min at room temperature to generate a single cell suspension. Next, the cell suspension was filtered using a Flowmi cell strainer (Millipore-Sigma), porosity 40uM to remove cell debris and clumps. To isolate the nuclei, cells were centrifuged at 500 x g for 5 min at 4 °C and then lysed using a 0.1X lysis buffer for 5 min on ice. The isolated nuclei were washed and resuspended in Diluted nuclei buffer. Finally, the nuclei concentration was determined using an automated cell counter and 15,000 nuclei were used to target a 10,000 nuclei recovery. The ATAC and gene expression library preparation and sequencing were performed at Cornell University BRC Genomics Core Facility according to the manufacturer’s instructions (10X genomics). Final libraries were sequenced on an Illumina NextSeq500 instrument.

10X multiome data analysis
Initial processing and alignment of paired scRNA-seq and scATAC-seq datasets from the 10X multiome kit was performed using the cellranger-arc pipeline (version 2.0.0)74. Upon generation of the UMI count matrix via cellranger-arc’s count function, downstream processing was performed using Seurat (version 4.3.0)75 and Signac (version 1.8.0)34. In addition to cellranger’s initial cell quality filtering, cells were additionally filtered in Seurat using the following criteria: nCount_ATAC < 8e4 & nCount_ATAC > 4e3 & nCount_RNA < 5000 & nCount_RNA > 300 & percent.mt < 10. In total, 7750 cells passed cell quality filtering. RNA- and ATAC-seq datasets underwent independent preprocessing and dimensionality reduction using the standard Seurat and Signac pipelines, respectively. Upon initial dimensionality reduction and clustering of the ATAC-seq dataset, peaks were re-called on individual clusters using MACS376 resulting in a final set of 316,412 peaks. The two datasets were then integrated by calculating a weighted nearest neighbor (WNN) graph, which was used for the final UMAP visualization and clustering. Clusters were empirically assigned using marker genes identified through gene expression and motif enrichment across UMAP clusters. Ordering of cells along the AP axis was accomplished using a quadratic programming approach developed by Yu and colleagues35. Extraembryonic clusters were excluded from the cell ordering analysis, as they were not predicted to be restricted along the AP axis. Cells were ranked according to their accessibility at anterior-/posterior-enriched peaks derived from TradeSeq analysis of the bulk tomoATAC-seq datasets. Cell ranks were assigned using custom R scripts from Yu and colleagues35. Smoothed gene expression profiles along the AP axis were computed using TradeSeq’s predictSmooth() function. Peak-gene linkage analysis was performed by computing the correlation between gene expression and accessibility at nearby peaks using Signac’s LinkPeaks function. Only high-confidence peak-gene links (p-value < 1 × 10-6) were used for downstream analysis. Analysis of individual germ layers was performed by subsetting the Seurat object to include only the clusters of interest, followed by repeating dimensionality reduction and WNN integration on each germ layer-specific data subset. AP differentially expressed genes in individual germ layers were identified using TradeSeq’s associationTest() function, and genes with an FDR < 0.05 were considered significantly enriched.

Mouse multiome data analysis
For the analysis of conserved rCREs in the mouse gastrula, 10913 anterior and 12532 posterior rCREs were lifted over to the mouse mm10 genome using UCSC genome browser’s LiftOver44 tool. This resulted in 1560 and 1724 conserved anterior and posterior rCREs, respectively. E7.5 mouse multiome data was obtained from Argelaguet et al.43 (GSE121708 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE121708]) and conserved rCREs were mapped onto ATAC-seq peaks in the multiome dataset using GenomicRanges findOverlaps() function. After omitting extraembryonic clusters, quadratic programming was used to order cells along the AP axis. Smoothed gene expression profiles along the AP axis were computed using TradeSeq’s predictSmooth() function.

Statistics & Reproducibility
Statistical analyses were performed using R (version 4.3.0). Statistical details for individual experiments have been provided in the figure legends and in the methods details. No statistical methods were used to predetermine sample size. In general, three biological replicates were used for genome-wide sequencing experiments.

Ethics statement
This work was conducted with ethical regulations and guidelines approved by Boston Children’s Hospital (IBC-RN00001327-1).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Azambuja and Rothstein et al present an interesting resource mapping the spatial epigenome of the gastrulating avian embryo. The methodologies seem sound and the data could be of great value to the community.\nMajor concern:\nMy main concern is whether the epigenetic information along AP axis accurately reflects the 'real' AP axis or rather the varied proportions of different cell types and their epigenome along the AP axis. This relates to the fact that tomoATAC-seq was performed only at seven transverse sections along the AP axis to establish positional CREs. The single cell 10X multiomic analysis was then compares to these pCREs to infer position. However, the tomoATAC-seq is a combination of all cells in the section. An alternative strategy would instead be to use the TOMO-seq reference which contains both AP and ML information to position the single cells epigenome.\nMy concern is raised from the observation that the SOX17 positive endoderm does not seem to be positioned correctly as they are exclusively placed in the anterior peak-set (Figure 5A). This is not the expected position and indeed it does not match the SOX17 IF in Figure 5D which clearly positions endoderm along the full AP axis.\nAn alternative but less elegant strategy to reduce the influence caused by different cell types is to check whether the pCREs display regional expression within the same cell types. This means that within each cell type, can a regional gradient be found and are these gradients the same in several cell types?\nMinor points:\nLine 145: \" Furthermore, we found a strong correlation between chromatin accessibility and deposition of H3K27AC at anterior and posterior pCREs, confirming that these genomic regions were not only open, but also displayed polarized chromatin activation in the embryo (Fig. S4, Fig. S5).\"\nNot clear to me that it is appropriate to draw this conclusion solely based on results from OTX2 and WNT5A? A systematic examination of the correlation between the gradient of chromatin accessibility and the deposition of H3K27AC needs to be conducted.\nLine165 \"Analysis of loci individual genes support that chromatin patterning precedes gene expression, as formation of pCRE gradients can already be detected in the loci of genes of DMBX1 and HOX11A at HH3, before transcription begins (Fig S7A-B). These results indicate pCRE patterning is not a consequence of regionalized gene expression\"\nAgain, it requires a statistical analysis to determine the relevance of gene expression, how many genes like DMBX1 and HOX11A, irrelevant to pCRE (or proceed to RNA expression), and how may genes such as OTX2 and WNT5B and their association with pCRE. And the reverse relationships.\nLine230:\"We employed peak-gene co-accessibility analysis in the multiome dataset to assign pCREs to their cognate target genes (Fig. 4C) and visualized the expression profiles of these genes along the AP axis. \"\nhow many total and pCRE peaks were recalled in scATAC-seq compared with the ones identified in TOMO-ATAC-seq?\nLine 600,\n\"Extraembryonic clusters were excluded from the cell ordering analysis, as they were not predicted to be restricted along the AP axis\"\nBased on the deconvolution strategy, it would be important check the proportion of extramembryonic tissues in TOMO-seq section. If it does have uneven distribution in different AP. How author exclude this difference from the analysis?", "reviewer_2": "This manuscript aims to generate a spatial map of gene expression (using tomo-seq) and epigenetic features such as chromosome accessibility (ATAC-seq) and bioinformatically predicted TF binding sites in chick embryos at primitive streak and early neurulation stages, to try to uncover spatial patterns of gene regulation in the embryo. The paper then explores some effects of exposing cells to activation of Wnt, as a putative posteriorizing factor. The data are important and interesting and should be a very useful resource for the community and to prompt further studies, coupled to embryological knowledge and future cell biology studies.\nThe title promised a very interesting and provocative story and given the lab's long track record of beautiful papers, I expected that it would contain a major discovery about gene regulatory patterns in space and time that would presage the acquisition of cell positional identity at a later moment in developmental time. Unfortunately, the paper does not deliver this. The main problem is some misinterpretation of the topology of the embryo at these stages in relation to time (ie. the space-time relationships of cells), which seriously affects the interpretation of the key findings and generating incorrect conclusions.\nMost of the studies have been done at HH stage 5, when the primitive streak has reached its full length and the head process (cranial part of the notochord) has emerged, with the prechordal mesendoderm at its tip (by the way this is not the \"gastrula\" stage but early neurula, because the future neural plate is already thickening at this stage and will start to elevate soon afterwards - but this is not so important). What is more important is that if an embryo at this stage is sectioned transversely (at right angles to the primitive streak and head process), these coordinates do not entirely correspond to the \"anterior-posterior\" axis in terms of cell identities or fates. From the tip of the primitive streak forwards, where the neural plate starts to form, yes indeed, this is the future head neuroectoderm and there the positions of cells do roughly correspond to their future rostrocaudal identities (forebrain to mid-hindbrain). But from the node down, the positions of cells at this stage do not correspond to the future \"anterior-posterior\" axis but much more closely to a temporal axis. It is true that the future tail will form at the most extreme posterior level of this, but the cells occupying this position in the epiblast as well as mesoderm at this stage are mainly lateral (mesoderm) and extraembryonic rather than posterior. Gradually as the streak shortens, the a/p axis of both the epiblast and the underlying mesoderm become more coincident, again the region in front of the node. In fact, some cells just emerging from the middle of the primitive streak at HH5 will end up in quite anterior positions in the mesoderm, as this tissue moves both away from the midline and anteriorly concomitantly to the narrowing and extension of the whole embryo. For these reasons, the sequential set of sections analyzed only correspond to the future a-p/cranial-caudal axis (in terms of cell identity) in the most anterior parts. Thus, the finding that in the most anterior sections there are many genes with differential expression corresponding to their final positions, whereas in the posterior sections there are very few, is not surprising - it is exactly what would be expected. It does not, as the text states, reveal an \"unexpected regulatory complexity in the anterior gastrula\". Genes like Cdx4 and Wnt5A \"move\" back with the elongating axis, marking a moving state rather than a cell identity - as cells move out of this domain they switch off these markers and acquire their positional characteristics. In other words, the a-p axis of the embryo at this stage, at least from the node downwards, is more like a temporal axis than a spatial one. So a series of cartesian coordinates does not accurately represent the embryonic organization or the future identity of the cells located in different positions, and different sections are not equivalent to each other.\nSimilar problems apply to the earlier embryos studied. At what the authors call \"HH stage 1\" [which actually no longer exists - this stage has been replaced by a series of 14 stages described by Eyal-Giladi and Kochav in 1976 so HH1 is no longer used/useful], cells are undergoing major reorganization within the epiblast and the most lateral cells in the disc will come to occupy the most \"posterior\" positions at the time of primitive streak formation. These are the well described \"polonaise\" movements. So again, a series of cartesian coordinates of embryos at these stages does not reflect future positional identities.\nThe \"a-p patterning\" literature, mainly in amphibians, contains evidence for at least 3 separate signals involved in \"posteriorization\": Wnt, FGF and retinoic acid. They probably are all involved in this, but in different regions of the axis and possibly in different germ layers. (in fact, such \"posteriorization\" is in principle incompatible with the authors' chosen alternative model, Mangold's idea of separate organizers for different parts of the axis). The paper looks at aspects of Wnt and its relation to their proposed pCREs (positional regulatory elements). Interesting and useful data using CUT&RUN with the Wnt target LEF/CTNNB1 are provided, in 3 regions of the embryo. Wnt signals are enriched in the most posterior but again this is likely to be related to the comparatively younger state of the epiblast and mesoderm in these regions relative to the maturing early neural plate at the other end and the authors' experiments do not allow these to be disentangled. The remaining experiments are done on dissociated cells in vitro (from the anterior part of the embryo at HH2, corresponding to the early gastrulation stage when the primitive streak starts to form), treated with Wnt modulators. There are differences between gain and loss of function, not surprisingly, but again the biological meaning of these changes cannot be resolved in terms of temporal/maturation versus axial identity. Moreover the cells appear to include all germ layers which complicates things further. The choice of tissue as the source of these cells is rather odd - this is not the region of the embryo that is likely to be patterned (posteriorized) by Wnt, certainly not the mesoderm, since at HH2 the anterior half of the embryo has very few if any mesoderm cells. One possibility is that the treatments can generate some mesoderm from epiblast. The hypoblast is also likely to be present in these dissociated samples. So very difficult to interpret unambiguously.\nFinally, the authors use single cell multi-omics using 10X genomics for RNAseq and ATAC. Another useful dataset, but also again a problem with interpretation. It is not surprising that there are clusters of like-cells, some of them probably identifiable from their signatures, but these results fit with many different interpretations. The authors claim that they support the proposal in the title.\nI think the authors have been led astray by these results. Interestingly, a similar problem affects another relatively recent paper (Metzis et al., cited in this paper) partly for the same reasons (\"anterior\" in the paper is mainly neural plate, \"posterior\" is mainly a mixture of things, and what the paper calls \"neural induction\" is in fact the start of the later events of neuronal differentiation, characterized by the expression of neuronal, rather than neural plate, markers). I bring this up in this review because it is a good example of how the title of a paper can be taken for granted and propagate through the literature.\nI do feel that this study does contain some very interesting data, and that the authors have indeed opened up an interesting question that is certainly worth pursuing. But I think that in order to do this appropriately, experiments need to be designed that take into account the temporal and spatial organization of the embryo itself as a guide, rather than imprint upon it human-generated systems of orthogonal coordinates that do not correspond to the underlying biological organization. I encourage the authors to re-think this story. The results are worth publishing in some form, but the interpretation needs to be seriously revised. The question of whether there is a general pre-pattern of chromatin organization that prefigures axial patterning, as already demonstrated 40 years ago for the Hox clusters, is indeed a good one and the authors are encouraged to explore this further in future studies.", "reviewer_3": "To discover whether and how positional information is recorded in the epigenome during embryo pattern formation, the authors took advantage of the large size of avian embryos, generated abundant data from spatial transcriptomic and epigenomic techniques, and performed comprehensive analyses. In the avian gastrula, the authors observed that some positional CREs (pCREs) formed gradients of chromatin accessibility (assessed by ATAC-seq) along the anterior-posterior axis. By unbiased clustering, they identified two sets of pCREs: anterior and posterior pCREs, on which the chromatin accessibility and the chromatin activity (assessed by H3K27ac CUT&RUN) are significantly higher at the anterior and posterior poles of the gastrula, respectively. They claimed that polarization of chromatin activity on these two pCRE sets precedes the spatial patterning of gene expression, and the pCRE gradients are shaped by signaling pathways, such as the WNT pathway, at the onset of gastrulation. Interestingly, the authors found that anterior and posterior pCREs can be used to predict the location of cells following the anterior-posterior axis not only in the entire embryo but also within individual germ layers. Overall, this manuscript is well prepared and nicely illustrated. However, there are several points the authors need to address to enhance the organization of the manuscript and better substantiate their conclusions.\n1. line 143-148: Data indicated here (Fig. S4, Fig. S5, Fig. 2G) are merely generated from eight peaks in the OTX2 locus and eight peaks in the WNT5A locus. The authors need to include data from all anterior and posterior pCREs, such as Fig. S6. Moreover, can the authors explain what the criteria were for selecting these 16 peaks?\n2. line 148-150: Can the authors show evidence of this restriction? If not, this conclusion should be removed.\n3. line 160-162: Without more data or statistical analyses, it’s difficult to tell if AP gene sets get regionalized following/after the establishment of pCRE gradients in Fig. S6. Please elaborate.\n4. line 165-166: This conclusion is overstated if only data on individual gene loci are shown.\n5. line 177: Why were HH5 embryos used for LEF1/CTNNB1 CUT&RUN? To test whether WNT signaling factors regulate pCRE gradients, CUT&RUN needs to be performed no later than when pCRE is spatially activated, which is the HH3 stage according to line 158. Otherwise, it’s difficult to distinguish between direct effects and indirect effects.\n6. line 183-184: This conclusion is overstated considering my previous comment.\n7. line 199-200: Statistical analysis is needed to support this claim.\n8. Fig. 2E: Please highlight anterior and posterior pCRE clusters and move other clusters to supplemental figures. Moreover, the authors are encouraged to include additional analyses to strengthen the biological significance of anterior and posterior pCRE clusters, such as the expression positioning of pCRE-associated genes (current Fig. 4C-D), the relationship between pCRE-associated genes and the 706/294 anteriorly/posteriorly enriched genes, as well as GREAT GO analysis of pCREs.\n9. Fig. 3J,K: Please separately indicate WNT-responsive and -inhibitory peaks/profiles.\n10. Fig. 4C,D: These two panels are suggested to move to Fig. 2 or Fig. 2 related supplemental figures.\n11. Fig. 5D: Mislabeled markers.\n12. Fig. 5F,G: These two panels do not belong to this figure.\n13. Fig. S6: Mislabeled y-axis in the right part."}, "author_response": "We thank the reviewer for underscoring the value of our manuscript for the community. We appreciate the comments regarding the potential confounding influence of varied cell-type proportions within each tomoATAC-seq section. The reviewer accurately notes that there is a higher number of cells assigned to the endodermal cluster in the anterior regions. To test whether this was due to an artifact in the pCRE- based positioning - or if it reflects the genuine spatial distribution of these cells - we took a closer look at the distribution of endodermal cells in sagittal sections of the gastrula- stage embryo.\nThis analysis, based on immunostaining of HH5 embryos using markers of the three germ layers (SOX17 for endoderm, SOX2 for ectoderm, and TBXT for mesoderm), is now presented in Figure 5D-E. It shows that the anterior enrichment of endodermal cells is due to the presence of the prechordal plate, an anteriorly located structure where the endoderm forms a thicker, multi-layered tissue compared to the single-layer endoderm observed in more posterior sections. The prechordal plate thus contributes to the high number of endodermal cells to the anterior of the embryo. In addition, we found that the posterior-most regions of the embryo at this stage (HH5) do not yet contain endodermal cells, as the endoderm has not completed its posterior migration, further contributing to this anterior bias.\nThus, the distribution patterns we observed in the embryo are consistent with the distribution of cells positioned based on the pCRE signature. We hope this added context and clarification address the reviewer’s concerns regarding AP axis positioning and the distribution of SOX17-positive cells.\nAn alternative but less elegant strategy to reduce the influence caused by different cell types is to check whether the pCREs display regional expression within the same cell types. This means that within each cell type, can a regional gradient be found and are these gradients the same in several cell types?\nWe thank the reviewer for this suggestion. To address this, we elaborated our analysis to check whether positional cis-regulatory elements (pCREs) exhibit region- specific expression within a single cell type, specifically focusing on neural ectoderm cells. To test for regional gradients within this cell type, we conducted a pCRE-based cell ordering in the neural ectoderm subcluster. We then assessed the expression of key CNS axial markers. This analysis revealed that ordered cells displayed the expected expression profiles for these markers along the AP axis, confirming that pCREs reliably encode positional information within a single tissue or germ layer (Figure 5G). Our results thus support the conclusion that pCREs maintain cellular positional identity within specific cell types.\nFurthermore, we visualized the distribution of pCREs in germ layer-specific subclusters of our multiome data. Projection of the anterior- and posterior-enriched pCREs onto germ layer-specific subclusters revealed distinct subpopulations of cells enriched for axial-restricted peaks in each of the three germ layers (Figure S11A-C). This observation indicates that spatial chromatin patterning can be observed consistently across different germ layers, suggesting that positional information is encoded in the epigenome in a manner that is largely independent of cell identity.\nMoreover, we conducted an analysis in which we performed trajectory-based differential expression analysis via Tradeseq in both the ectodermal and mesodermal subclusters of our multiome dataset (2481 cells and 2936 cells, respectively). This analysis revealed a set of 64 genes significantly associated with anterior-posterior position (fdr < 0.05) in both germ layers (Figure S11D-E).\nTogether, these findings indicate that pCREs encode region-specific positional information within cell types and that similar AP information is present across the germ layers.\nMinor points:\nLine 145: \" Furthermore, we found a strong correlation between chromatin accessibility and deposition of H3K27AC at anterior and posterior pCREs, confirming that these genomic regions were not only open, but also displayed polarized chromatin activation in the embryo (Fig. S4, Fig. S5).\"\nNot clear to me that it is appropriate to draw this conclusion solely based on results from OTX2 and WNT5A? A systematic examination of the correlation between the gradient of chromatin accessibility and the deposition of H3K27AC needs to be conducted.\nWe thank the reviewer for this suggestion. To address this, we have included plots showing chromatin accessibility and H3K27ac occupancy across the AP axis for all anterior and posterior pCREs (Figure 2B). Both data modalities exhibit clear anterior/posterior enrichment patterns. In addition, we performed a statistical analysis examining the correlation between chromatin accessibility and H3K27ac deposition across all peaks in our dataset. This analysis confirmed a strong global correlation between chromatin accessibility and H3K27ac enrichment in both anterior and posterior pCREs (Figure S3G, p-value<2x10-16).\nLine165 \"Analysis of loci individual genes support that chromatin patterning precedes gene expression, as formation of pCRE gradients can already be detected in the loci of genes of DMBX1 and HOX11A at HH3, before transcription begins (Fig S7A-B). These results indicate pCRE patterning is not a consequence of regionalized gene expression\"\nAgain, it requires a statistical analysis to determine the relevance of gene expression, how many genes like DMBX1 and HOX11A, irrelevant to pCRE (or proceed to RNA expression), and how may genes such as OTX2 and WNT5B and their association with pCRE. And the reverse relationships.\nWe chose to highlight DMBX1 and HOXA11 in Figure S7 because these are among the very few AP-restricted genes that are not yet transcriptionally active at stage HH3, a developmental timepoint where we already detect robust graded activation of pCREs. The reviewer is correct in suggesting that a systematic analysis correlating chromatin patterning with the onset of gene expression would strengthen the conclusion. However, most AP-restricted genes in our dataset—such as OTX2 and WNT5A—are already expressed at HH3, limiting our ability to assess, on a genome- wide scale, whether chromatin activation precedes transcriptional activation for this specific class of genes. Since only a small subset of genes displays delayed activation relative to chromatin remodeling, we lack the statistical power to perform a broader analysis of this phenomenon.\nIn light of the reviewer’s comment, we have revised the text to soften our interpretation and avoid overgeneralizing this observation. We now present these cases as specific examples supporting the notion that pCRE activation can, in certain cases, precede gene expression, while acknowledging that this is not uniformly applicable across all loci.\nLine230:\"We employed peak-gene co-accessibility analysis in the multiome dataset to assign pCREs to their cognate target genes (Fig. 4C) and visualized the expression profiles of these genes along the AP axis. \"\nhow many total and pCRE peaks were recalled in scATAC-seq compared with the ones identified in TOMO-ATAC-seq?\nFor analysis of all bulk and single-cell ATAC-seq datasets, a master set of 316,412 peaks were used. This set was originally obtained from calling peaks on individual cell clusters of the scATAC-seq portion of the HH5 multiome dataset and concatenated to generate a comprehensive set of all possible open chromatin regions in the HH5 chicken embryo. By calling peaks on individual cell clusters, this increased the sensitivity of peak calling and ensured rare, cell type-specific regulatory elements were not omitted from the analysis. This ‘master’ peak set was subsequently used for all bulk ATAC-seq datasets, including the tomo-ATAC-seq in Figure 2, as well as the epiblast ATAC-seq in Figure 3. To ensure that pCREs were, in fact, identified as accessible in our bulk tomo-ATAC-seq datasets, we performed de novo peak calling on these datasets using MACS3. Of the 10913 and 12532 anterior and posterior pCREs originally identified in our initial analysis, 8525 and 10271 were identified using de novo peak calling of the bulk tomo-ATAC-seq datasets. Thus, the majority of pCREs were independently identified as accessible in both the bulk and scATAC-seq datasets. We have now clearly stated our peak calling strategy and rationale in the methods, as well as the body of the manuscript.\nLine 600,\n\"Extraembryonic clusters were excluded from the cell ordering analysis, as they were not predicted to be restricted along the AP axis\"\nBased on the deconvolution strategy, it would be important check the proportion of extramembryonic tissues in TOMO-seq section. If it does have uneven distribution in different AP. How author exclude this difference from the analysis?\nBased on our multiome data, we estimate that approximately 13.8% of the cells are extraembryonic. Importantly, we observed that these extraembryonic cells exhibit low levels of pCRE accessibility compared to embryonic cells (as shown on Figure 5B), indicating that they do not substantially contribute to or interfere with the detection and quantification of AP-axis-specific pCREs.\nWe thank the reviewer for their thoughtful and constructive comments. We would like to clarify that the primary aim of our study is not to define the final cell identities or static positional outcomes of embryonic cells, but rather to uncover how the genome is pre-patterned at early stages of axial specification. Our findings reveal that chromatin accessibility and activation are spatially organized across the embryo, indicating that the chromatin landscape encodes positional cues well before cells acquire their definitive fates.\nWe propose that this early spatial patterning of chromatin represents a regulatory mechanism that may prime cells for subsequent identity acquisition along the anterior-posterior axis. Our approach focuses on this foundational level of genomic organization, and we agree that this does not directly equate to final cell positions or identities. In light of the reviewer’s feedback, we have clarified this distinction throughout the manuscript, as it underscores the core contribution of our work: the identification of an early spatial chromatin code that likely informs subsequent patterning events.\nWe thank the reviewer for this helpful clarification regarding staging terminology. We have updated the manuscript to describe HH stage 5 as the “late gastrula” stage, acknowledging that this timepoint marks the transition toward early neurulation, as the neural plate begins to thicken.\nWe thank the reviewer for this thoughtful and detailed discussion of cell positioning dynamics during early development. We fully agree that, due to extensive tissue remodeling and cell movements during gastrulation, the transverse sections used in our tomoATAC-seq analysis do not strictly map onto the final, mature anterior- posterior (AP) axis - particularly in more posterior regions where the relationship between initial position and final cell fate is less direct.\nHowever, we would like to clarify that our study does not aim to assign definitive rostrocaudal fates to cells at this stage, nor to equate our sections directly with the final AP axis. Although we refer to anterior and posterior coordinates at HH5, we recognize that these positions are transient in the rapidly changing gastrula/early neurula. Thus, our use of the term ‘positional identity’ here refers to the chromatin state and the cellular position at HH5, rather than a definitive or final fate. Our analysis of the three developmental stages inform upon how epigenetic landscapes are initially configured, even though we agree that ongoing regulatory processes will continue to reshape the embryo’s spatial coordinates.\nThe central contribution of our work is the identification of an early chromatin- based spatial code, which precedes and likely contributes to the later establishment of gene expression domains along the axis. We now emphasize in the revised manuscript that this chromatin organization should be interpreted as an early regulatory blueprint for positional information, rather than a direct representation of terminal cell fates or final AP patterning. We appreciate the reviewer’s insights, which have prompted us to better contextualize our findings within the dynamic spatial-temporal framework of early embryonic development.\nWe have clarified in the manuscript that we used EGK XIII stage embryos for our early genomic analyses. We fully recognize the reviewer’s point regarding the dynamic tissue movements and rearrangements, such as polonaise movements and primitive streak elongation, which decouple early spatial coordinates from future positional identities, especially in the posterior embryo.\nTo address this, we have expanded the discussion to more explicitly state that our analysis pertains to the spatial organization of chromatin accessibility as it relates to cell position at HH5, rather than to the ultimate fate or final positioning of cells along the anterior-posterior axis. Our data indicate that pCRE gradients at HH5 reflect the positional information of cells at this specific stage, prior to the extensive reorganization that will define the mature body plan.\nAdditionally, our findings show that pCRE patterning is shaped by signaling pathways such as WNT, supporting a model where chromatin accessibility integrates both positional cues and dynamic signaling inputs as cells migrate and reposition within the embryo. While we agree that axial elongation introduces further complexity, we do not interpret our data as reflecting the definitive future axial identity of cells, but rather as capturing an early, stage-specific chromatin landscape that encodes positional information at HH5.\nWe appreciate the reviewer’s insights, which have helped us better articulate this distinction and to more carefully frame our conclusions in the revised manuscript.\nWe thank the reviewer for their thoughtful comments and the opportunity to clarify our experimental approach. In these experiments, we specifically isolated epiblast cells from the anterior region of the embryo at HH2, with the hypoblast carefully removed to ensure that it was not present in our dissociated samples. We selected epiblast cells at this stage because they do not yet display chromatin patterning, making them ideal for examining the initial effects of Wnt signaling on chromatin accessibility.\nImportantly, we did not observe activation of mesodermal or other germ layer- specific genes following Wnt treatment. Rather, we see expression of posterior genes, such as components of the HOX code. Our results therefore demonstrate that posteriorization by Wnt signaling occurs alongside an increase in accessibility at posterior pCREs, rather than to infer specific germ layer differentiation. We appreciate the reviewer’s insights, which have helped us clarify the interpretation and significance of these findings.\nOur analysis demonstrates that gradients of chromatin accessibility can accurately predict the anterior-posterior (AP) position of individual cells within the embryo. This predictive power supports our main proposal - that chromatin accessibility patterns establish a spatial framework reflective of positional information along the AP axis.\nIn response to the reviewer’s comments, we have revised the discussion section to more thoroughly address additional factors influencing our observations, including cell movements and the temporal progression of cell identities during early embryogenesis. While we recognize the complexity these dynamics introduce, we believe that the discovery of a spatial code embedded in chromatin accessibility along the embryonic axis constitutes a significant finding.\nWe appreciate the reviewer’s valuable perspective on the need to more closely align our interpretations with the embryo's intrinsic spatial and temporal organization. Accordingly, we have better contextualized our conclusions regarding chromatin pre- patterning and clarified the scope and limitations of our study. Nonetheless, we believe that our findings offer new insights into how chromatin organization may interface with embryonic patterning mechanisms.\nWe share the reviewer’s view that this opens a compelling line of investigation, and we are actively considering follow-up experiments that integrate additional spatial and temporal datasets to further refine our understanding of this relationship. We are grateful for the reviewer’s encouragement to continue exploring this important area.\nWe thank the reviewer for their positive assessment of our manuscript and for highlighting the key findings and strengths of our study. We appreciate the reviewer’s comments regarding areas for improvement and clarity, and we have carefully addressed each of these points to enhance the organization of the manuscript and to further substantiate our conclusions. We believe these revisions strengthen our study and improve its presentation, and we are grateful for the reviewer’s constructive feedback.\n1. line 143-148: Data indicated here (Fig. S4, Fig. S5, Fig. 2G) are merely generated from eight peaks in the OTX2 locus and eight peaks in the WNT5A locus. The authors need to include data from all anterior and posterior pCREs, such as Fig. S6. Moreover, can the authors explain what the criteria were for selecting these 16 peaks?\nWe thank the reviewer for this suggestion. In response, we have now included profiles of all anterior and posterior pCREs to provide a more comprehensive view of chromatin accessibility across all anterior/posterior pCREs (see Figure 2B). We now explain in the manuscript that pCREs associated with OTX2 and WNT5A shown on Figure 2 were selected based on their proximity to these genes' transcriptional start sites. These loci were chosen to exemplify how specific gene regions exhibit spatial chromatin patterning during gastrulation. We hope this additional data helps to further clarify the spatial organization of chromatin accessibility across the anterior-posterior axis.\n2. line 148-150: Can the authors show evidence of this restriction? If not, this conclusion should be removed.\nTo address this, we have conducted a gene enrichment analysis of pCREs (Figure S3B). Through this analysis, we identified a significant association between genes enriched for pCREs in their loci and those with anterior/posterior-enriched expression patterns (P=3.635e-06, chi-square test). Notably, while there is a significant association between pCRE patterning and AP restricted expression, this is not always the case. There are many genes with AP restricted expression patterns that are not enriched for pCREs. In addition, pCREs can also be enriched in the loci of genes that do not display clear axial patterning. To better reflect this observation, we have revised the manuscript to clarify that pCREs are abundant in the loci regionally expressed genes, rather than being restricted to such loci. We believe this revision more accurately captures the complexity of our data and fully addresses the reviewer’s concern.\n3. line 160-162: Without more data or statistical analyses, it’s difficult to tell if AP gene sets get regionalized following/after the establishment of pCRE gradients in Fig. S6. Please elaborate.\n4. line 165-166: This conclusion is overstated if only data on individual gene loci are shown.\nAs mentioned above (see rebuttal to Rev. 1), DMBX1 and HOXA11 represent a small subset of AP-restricted genes that are transcriptionally inactive at the stage (HH3) when pCRE gradients first emerge. The majority of AP-patterned genes in our dataset, such as OTX2 and WNT5B, are already transcriptionally active at this stage, limiting the pool of genes available to systematically assess whether chromatin patterning broadly precedes transcription.\nDue to the scarcity of genes that are both AP-restricted and not expressed before pCRE patterning, a genome-wide statistical analysis to test whether pCRE gradients generally predate gene activation lacks sufficient power. As a result, we focused on DMBX1 and HOXA11 as representative cases where early chromatin accessibility clearly prefigures transcriptional onset. In response to the reviewer’s comment, we have edited the text to clarify this limitation and to explain why only a small number of loci are suited to this specific analysis. We have also softened our conclusions, now presenting these findings as case studies that suggest early chromatin organization can, in some instances, precede and potentially prime gene expression, while avoiding an overgeneralization across the genome.\n5. line 177: Why were HH5 embryos used for LEF1/CTNNB1 CUT&RUN? To test whether WNT signaling factors regulate pCRE gradients, CUT&RUN needs to be performed no later than when pCRE is spatially activated, which is the HH3 stage according to line 158. Otherwise, it’s difficult to distinguish between direct effects and indirect effects.\nWe thank the reviewer for this valuable suggestion. In response, we have now performed LEF1 and CTNNB1 CUT&RUN at stage HH3, when pCRE gradients first emerge, to better assess potential direct regulation of these elements by WNT signaling factors. The results are consistent with our model, showing that LEF1/CTNNB1 binding is spatially enriched at early-activated pCREs. This new data is now included in Figure S8 and supports the notion that WNT signaling contributes directly to the early establishment of pCRE gradients.\n6. line 183-184: This conclusion is overstated considering my previous comment.\nThe results from the new experiment requested by the reviewer (Figure S8) support our conclusion that WNT signaling factors directly contribute to the regulation of pCRE gradients at this early stage.\n7. line 199-200: Statistical analysis is needed to support this claim.\nWe have now included additional statistical analysis demonstrating that LEF1/CTNN1B-bound peaks are significantly associated with increased accessibility upon WNT GOF, as well as decreased accessibility upon WNT LOF (Figure 3H). Statistical analysis of WNT-responsive peaks (peaks increasing accessibility upon WNT GOF and decreasing accessibility upon WNT LOF) also showed a significant association with LEF1/CTNN1B binding (Figure 3I).\n8. Fig. 2E: Please highlight anterior and posterior pCRE clusters and move other clusters to supplemental figures. Moreover, the authors are encouraged to include additional analyses to strengthen the biological significance of anterior and posterior pCRE clusters, such as the expression positioning of pCRE-associated genes (current Fig. 4C-D), the relationship between pCRE-associated genes and the 706/294 anteriorly/posteriorly enriched genes, as well as GREAT GO analysis of pCREs.\nWe have now replaced the former Figure 2E with a focused view of the anterior and posterior pCREs (now Figure 2B) and moved the other clusters to the supplemental figures for clarity. Additionally, we have performed a GREAT GO analysis of pCREs to reveal functional enrichments of genes associated with these elements (Figure S3E). These added analyses provide greater context for the biological significance of the anterior and posterior pCRE clusters. We note that the analysis presented in Figure 4C- D remain in place, as they are dependent upon the multiome analysis presented in Figure 4.\n9. Fig. 3J,K: Please separately indicate WNT-responsive and -inhibitory peaks/profiles.\nWe have modified Figures 3J and 3K to separately indicate WNT-responsive and WNT-inhibitory peaks/profiles. This adjustment provides a clearer visualization of the differential effects of WNT signaling on pCRE accessibility.\n10. Fig. 4C,D: These two panels are suggested to move to Fig. 2 or Fig. 2 related supplemental figures.\nWe thank the reviewer for this suggestion. However, we believe that these panels are best kept in their current location, as they directly depend on the multiome analysis introduced in Figure 4.\n11. Fig. 5D: Mislabeled markers.\nWe thank the reviewer for catching this mistake. It has been fixed.\n12. Fig. 5F,G: These two panels do not belong to this figure.\nWe thank the reviewer for this suggestion. We believe the comparative analysis showing that pCREs can also predict cell position in the mouse embryo is a key component of the study. This result demonstrates that spatial chromatin patterning is conserved in vertebrate embryos and extends the relevance of our findings beyond the chick model. Given the importance of highlighting this evolutionary conservation, we think it is appropriate to retain this analysis in a main figure instead moving it to supplemental data.\n13. Fig. S6: Mislabeled y-axis in the right part.\nWe have fixed this mistake."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have sufficiently addressed my concerns.", "reviewer_3": "The majority of my concerns have been adequately addressed. I would like to congratulate the authors for a job well done!"}, "author_response": "We are pleased that the revised manuscript addresses the concerns raised and we appreciate the reviewer’s support for the current version of the study.\n[TEXT REDACTED]\n[TEXT REDACTED]\nWe thank the reviewer for their kind words and constructive feedback throughout the review process. We are glad that the revised manuscript addresses the majority of the concerns raised and appreciate the encouraging assessment."}, {"round": 3, "reviewer_comments": {"reviewer_4": "The authors have addressed most of the concerns/issues raised by Reviewer 2 in the rebuttal and with the revision of the manuscript.\nThere are issues on “the significance of these putative coordinates” and “toned down interpretation” that may need further clarification.\nSignificance of the putative coordinates (or “positional information is encoded within the chromatin landscape” – Abstract).\nAs previously pointed out by Reviewer 1, the information (regionalised gene expression pattern and active CREs in the accessible chromatin region) collated from the TomoRNA-seq and TomoATAC-Seq datasets reflects the composition of cell types of a multitude of cell populations localized in different regions of the HH5 avian embryo. The CRE profile does not by default carry any positional information of the descendants of the cell population along the A-P axis in the final body plan. The spatial pattern of the expression of the TomoRNA-seq gene set that is characteristic of the cell types in the regionalized cell populations of interest and together with the sets of CREs that display peak-gene co-accessibility (e.g. Fig 4C, D) provided the spatial register (positional address) of the different cell types in the whole embryo or the germ layer (e.g., Fig. 5). Based on this register, single cells from the embryos at the same developmental stage can be mapped to the inferred positions in the embryo (re Fig. 4- the legend may be revised as “region-related CREs enable inferring the position of cells in the embryo”). The primary attribute (“significance”) of the spatial register, or measurable positional values information (line 67) is for inferring (not “can predict” – Line 67) where in the germ layers individual cells reside within the embryo (e.g., Fig 4B, 5C). This approach meets the requirement of identifying the spatial position of cell clusters in the embryo, without resorting to high-resolution spatial multi-omics analysis. Mapping of single cells to the specific germ layer population in the embryo facilitates the clustering of the single cells by markers genes (e.g. Fig S10), and the follow-on elucidation of the enrichment of CREs in these populations (e.g. Fig. S11).\nThe terminology anterior pCREs and posterior pCREs potentially lead to misconstruing the measurable positional values reflects the final positional address of the specific cell population in the body axes (as articulated in the first two paragraphs of Introduction). To avoid such misconstruing, it may be appropriate to revise the Introduction by removing the first two paragraphs and go straightly with the third paragraph. Since the “anterior: a” and “posterior: p” descriptors are based on the anatomical/topographical position of the cell population with refence to a morphological landmark (Hensen’s Node of HH5 embryo), it may be appropriate to name these CREs as regional CRE (rCRE) with designation of aN-rCREs and pN-rCREs (aN and pN stand for anterior to node: AP1, AP2, and posterior to node: AP3-7 Fig 2A). This would help distinguishing these regional CREs from those that are truly positional CREs (pCRE) of axial identity genes (e.g., pCREs of members of HOX gene clusters displaying different expression domain in the A-P axis, Fig. 1F, 2G, 3E), in the discussion about the axial information recorded in the chromatin landscape / epigenome (line 176-180, line 269-274, Fig 4F, G - HH10 embryo)\nToned down the interpretation of the CRE functionality.\nThe most impactful findings of this study are the characterization of the epigenome landscape/architecture of the germ layers of HH5 embryo at spatial resolution, and the identification of the CREs of different functionality (e.g., Supp Fig S3). In this context, while arguably that functional analysis of HH5 CREs is outside the scope of this study (apart from the demonstration of changes in CRE activity of HH2 epiblast cells in response to modulated WNT signalling activity), further functional annotation of the CRE sets would enhance the merit of this study on CRE functions. In a similar study of mouse gastrula embryo, a range of CRE functionality has been identified, including driving cell lineage development, anterior-posterior patterning, response to long-range signalling activity and short-range signalling activity and cell-cell interaction, midline mesendoderm formation, mesoderm patterning, cell cycle control and others. Of interest is the scope of functionality of the germ layer-specific CREs, as well as specifically those enriched CRs that are related to genes that display overlapping spatially-restricted expression across different germ layers (e.g., Fig S11D for ectoderm-mesoderm, yet missing the data for ectoderm-endoderm and endoderm-mesoderm). Together, knowledge gleaned from this study would deeper insights into the molecular architecture and functional attributes of the gene regulatory networks of lineage specification and differentiation and tissue patterning in the developing embryo.\nThis additional information may help to justify the interpretation of CREs of developmental genes that displayed spatially restricted expression domain actively record the positional information shared by distinct cell types and across populations in different germ layers (line 345-354) and that CREs regulated by nuclear effectors of signalling pathways provide the readout of the mechanistic guidance for tissue patterning in the A-P axis (line 359-374).\nRegarding the assertion that results of this study inform how cells are allocated to positional (anterior or posterior) fate during development, further clarification may be required to identify which results are showing the spatiotemporal changes in epigenomic states for both anterior and posterior axis identity (or cell fates) are established early in development (line 375-382, line 390-392), and that the epigenome stores the positional information, as records of morphogen-derived signals into chromatin accessibility landscape or positional references for downstream GRN activity (for early specification of correct axial cell types, line 396-399), during early development. Of note: the developmental series of embryos analysed for this study was incomplete (Fig. S6), which was missing HH4 when further activation of CREs (both those that are concomitant with expression of target genes and preceding gene expression, e.g. Fig. S7) may happen.\nIt is of interest to demonstrate conservation of CRE signatures, but the study may go further to demonstrate conservation of functionality to account for species-specific difference in functionality of certain conserved CREs. As the comparative study covered only two species (both are amniote), it would not be appropriate to extend the findings to conservation across all vertebrates. Reviewer 3 previously queried the merit/relevance of the conservation study. Given that the study was limited to two species and it is unclear if the tissue constitution of the germ layers and the spatial pattern of gene expression are comparable/conserved between HH5 avian embryo and E7.5 mouse embryo, it may be appropriate to remove these results (and discussion, line339-403) and include them for a more in-depth study across a range of developmental stages of these two species as well as including other vertebrate model organisms.\nSpecific points:\nLine 231: Suggest revising to Chromatin gradients infer axial position ….\nLine 276: Suggest revising to Regional CREs record spatial register of cells in the germ layers.\nLine 294-296: The absence of “endoderm” in the posterior region may not be due to incomplete migration. It is more likely that the definitive endoderm has yet to be specified and the “endoderm” in the posterior region is the descendant of the hypoblast.\nFig. 1E, Fig S2A: Could the concordance of in silico gene expression pattern and in situ Hybridization pattern be evaluated for statistical significance to determine how well the in-silico data recapitulate the in-situ data?\nFig 2B: Define the colour code (of solid lines) for data of chromatin accessibility and histone enrichment and use the same colour code for data presented in other panels and Fig. 2A.\nFig 5: Revise the legend to not using the terms pCREs and positional information. 5C: Is this schematic only depicting cell types in the midline? Which structures are prechordal plate and head process? Should the “notochord” be the head process?\nFig S3C: Only the data of six clusters were shown but not seven clusters as indicated in the legend."}, "author_response": "We appreciate the reviewer’s assessment that the main issues raised in previous rounds have been resolved. In this revision, we have further edited the manuscript to clarify the significance of the inferred spatial coordinates and to temper the interpretation of CRE functionality. We have also revised the abstract to soften the interpretive language regarding positional information in the chromatin landscape.\nWe have revised the legend of Figure 4 to adopt the suggested phrasing, clarifying that “Regional CREs enable inferring the position of cells in the embryo.” We also ensured that the corresponding text consistently reflects that these elements are used to infer, rather than predict, spatial position based on chromatin accessibility patterns.\nWe have implemented the requested revisions, replacing “can predict” with “inferring” throughout the manuscript.\nWe appreciate the reviewer’s suggestion. We tested the removal of the first two paragraphs but found that they provide important conceptual background for the study. We have therefore retained them and revised the language to emphasize early embryonic patterning and the generation of spatial cues - avoiding any reference to the final positional fate of cells.\nWe thank the reviewer for this thoughtful suggestion. We have adopted the term regional CREs (rCREs) as recommended, to more accurately describe these elements. Since our analysis did not use the Hensen’s node as a direct anatomical landmark, we have designated these regions as anterior and posterior rCREs. The terminology has been updated accordingly in the figures, legends, and main text throughout the manuscript.\nWe thank the reviewer for these insightful comments. We fully agree that a comprehensive functional analysis of the identified CREs would provide valuable insights into their developmental roles. However, such an analysis is beyond the scope of the present study, which focuses on defining the spatial architecture of the epigenome at gastrulation. Our primary objective was to establish how chromatin accessibility is regionally organized across the germ layers rather than to functionally characterize individual CREs.\nTo address the reviewer’s suggestion and provide additional context for CRE functionality, we performed a gene ontology (GO) enrichment analysis of genes associated with anterior and posterior regulatory regions (included in Figure S10D). This analysis revealed biologically relevand trends, with anterior rCREs enriched for terms related to sensory organ development, telencephalon morphogenesis, and cell adhesion, processes consistent with anterior ectodermal and neuroepithelial patterning. In contrast, posterior rCREs showed enrichment for mesodermal developmental programs, including somitogenesis and segmentation. These results provide additional functional context for the spatial organization of the chromatin landscape during gastrulation\nImportantly. we have also toned down the language throughout the manuscript to ensure that our interpretation remains strictly supported by the data. We now emphasize that the identified rCREs reflect spatial organization of the chromatin landscape during gastrulation rather than directly mediate tissue patterning or positional information.\nWe thank the reviewer for pointing out the need for these clarifications. The spatiotemporal emergence of anterior and posterior chromatin states is shown in Fig. S6, which illustrates the progressive establishment of graded H3K27ac patterns at anterior and posterior rCREs during gastrulation. These data demonstrate that chromatin polarization is already evident at HH3, indicating that both anterior and posterior epigenomic states are established early in development. We have revised the relevant text to clarify that by stating the “epigenome stores positional information,” we refer to the effects of morphogen signaling being recorded as localized changes in chromatin accessibility, rather than long-term epigenetic memory. Finally, we acknowledge that our developmental series does not include HH4. However, rCRE polarization is already detectable at HH3, suggesting that the key transitions in chromatin accessibility precede HH4. Thus, the inclusion of HH4 would likely capture further refinement - but not the initial establishment of these epigenomic states.\nWe thank the reviewer for this comment. The aim of the comparative analysis was not to broadly examine species-specific differences, but rather to demonstrate that the identified rCREs are not unique to the chick embryo and display a notable degree of conservation within amniotes. This supports the interpretation that aspects of regional chromatin organization are evolutionarily preserved. We have revised the text to clarify this point and have toned down the language in the Results and Discussion to avoid overextending the conclusions beyond the scope of the comparison.\nThis has been revised according to the reviewer’s suggestion.\nWe have revised according to the reviewer’s suggestion.\nWe have revised the text to indicate that the absence of endoderm in the posterior region likely reflects that definitive endoderm has not yet been specified at this stage and that the posterior layer corresponds to hypoblast descendants, rather than incomplete migration.\nWe thank the reviewer for this suggestion. While we agree that assessing concordance between in silico and in situ datasets is important, a statistical comparison is not feasible in this context because in situ hybridization data is not quantitative. Hence the two data types differ in modality, precluding a direct one-to-one quantitative comparison. We are, however, developing approaches to improve the statistical assessment of spatial concordance in future analyses by integrating quantitative image analysis with spatial transcriptomic profiles.\nHere, we opted to preserve the original color code so that the ATAC-seq and CUT&RUN datasets could be plotted together in the same charts, allowing for direct visual comparison between chromatin accessibility and histone enrichment profiles.\nWe thank the reviewer for catching this oversight. Figure 5C depicts a sagittal section through the embryonic midline, and we have clarified this point in the figure legend. We also corrected the labeling to indicate the head process rather than the notochord.\nWe thank the reviewer for catching this error. The figure correctly displays six clusters, and we have updated the legend to reflect this."}]

Decision letter