Spatial epigenome–transcriptome co-profiling of mammalian tissues

Emerging spatial technologies, including spatial transcriptomics and spatial epigenomics, are becoming powerful tools for profiling of cellular states in the tissue context1–5. However, current methods capture only one layer of omics information at a time, precluding the possibility of examining the mechanistic relationship across the central dogma of molecular biology. Here, we present two technologies for spatially resolved, genome-wide, joint profiling of the epigenome and transcriptome by cosequencing chromatin accessibility and gene expression, or histone modifications (H3K27me3, H3K27ac or H3K4me3) and gene expression on the same tissue section at near-single-cell resolution. These were applied to embryonic and juvenile mouse brain, as well as adult human brain, to map how epigenetic mechanisms control transcriptional phenotype and cell dynamics in tissue. Although highly concordant tissue features were identified by either spatial epigenome or spatial transcriptome we also observed distinct patterns, suggesting their differential roles in defining cell states. Linking epigenome to transcriptome pixel by pixel allows the uncovering of new insights in spatial epigenetic priming, differentiation and gene regulation within the tissue architecture. These technologies are of great interest in life science and biomedical research.

Emerging spatial technologies, including spatial transcriptomics and spatial epigenomics, are becoming powerful tools for profiling of cellular states in the tissue context [1][2][3][4][5] . However, current methods capture only one layer of omics information at a time, precluding the possibility of examining the mechanistic relationship across the central dogma of molecular biology. Here, we present two technologies for spatially resolved, genome-wide, joint profiling of the epigenome and transcriptome by cosequencing chromatin accessibility and gene expression, or histone modifications (H3K27me3, H3K27ac or H3K4me3) and gene expression on the same tissue section at near-single-cell resolution. These were applied to embryonic and juvenile mouse brain, as well as adult human brain, to map how epigenetic mechanisms control transcriptional phenotype and cell dynamics in tissue. Although highly concordant tissue features were identified by either spatial epigenome or spatial transcriptome we also observed distinct patterns, suggesting their differential roles in defining cell states. Linking epigenome to transcriptome pixel by pixel allows the uncovering of new insights in spatial epigenetic priming, differentiation and gene regulation within the tissue architecture. These technologies are of great interest in life science and biomedical research.
Single-cell multiomics allows uncovering of the mechanisms of gene regulation across different omics layers [6][7][8][9] but lacks spatial information, which is crucial to understanding cellular function in tissue. Spatial epigenomics, transcriptomics and proteomics emerged recently [1][2][3][4][5] but most of these can capture only one layer of the omics information. Although computational methods can integrate data from multiple omics 10 , they cannot readily uncover the mechanistic links between different omics layers. Previously we developed deterministic barcoding in tissue for spatial omics sequencing for spatially resolved comeasurement of the transcriptome and a panel of proteins 1 , and recently it was also realized in the 10X Visium platform 11 . Herein, to further investigate the epigenetic mechanisms underlying gene expression regulation, we developed spatially resolved, genome-wide comapping of the epigenome and transcriptome by simultaneous profiling of chromatin accessibility and messenger RNA expression (spatial assay for transposase-accessible chromatin and RNA using sequencing (spatial ATAC-RNA-seq)), or histone modifications and mRNA expression (spatial assay of cleavage under targets and tagmentation and RNA using sequencing (spatial CUT&Tag-RNA-seq); applied to H3K27me3, H3K27ac or H3K4me3 histone modifications) on the same tissue section at the cellular level via deterministic cobarcoding, to integrate the chemistry of spatial-ATAC-seq 3 or spatial-CUT&Tag 2 with that for spatial transcriptomics. We applied these techniques to embryonic and juvenile mouse brain, as well as adult human brain hippocampus, to dissect epigenetic and transcriptional states and their role in the regulation of cell types dynamics in tissue. This work opens new frontiers in spatial omics and may provide unprecedented opportunities to biological and biomedical research.

Article
adaptor containing a universal ligation linker that can be inserted into transposase-accessible genomic DNA loci. The same tissue section was then incubated with a biotinylated DNA adaptor containing a universal ligation linker and a poly-T sequence that binds to the poly-A trail of mRNAs to initiate reverse transcription (RT) in tissue. A microfluidic channel array chip was then placed on the tissue section to introduce spatial barcodes Ai (i = 1−50 or 100) that were covalently conjugated to the universal ligation linker via templated ligation. Next, another microchip with microchannels perpendicular to the first flow direction was used to introduce spatial barcodes Bj ( j = 1−50 or 100) that were then ligated to barcodes Ai (i = 1−50 or 100), resulting in a two-dimensional grid of spatially barcoded tissue pixels, each defined by the unique combination of barcodes Ai and Bj (i = 1−50 or 100, j = 1−50 or 100; barcoded pixels, n = 2,500 or 10,000). Finally, barcoded complementary DNA and genomic DNA fragments were released after reverse crosslinking. cDNAs were enriched with streptavidin beads and gDNA

Fig. 1 | Design and evaluation of spatial epigenome-transcriptome
cosequencing with E13 mouse embryo. a, Schematic workflow. b, Comparison of number of unique fragments and fraction of reads in peaks (FRiP) in spatial ATAC-RNA-seq and spatial CUT&Tag-RNA-seq. c, Gene and UMI count distribution in spatial ATAC-RNA-seq and spatial CUT&Tag-RNA-seq. Number of pixels in E13, 2,187; in human brain, 2,500; in mouse brain (ATAC), 9,215; in mouse brain (H3K27me3), 9,752; in mouse brain (H3K27ac), 9,370; in mouse brain (H3K4me3), 9,548. Box plots show the median (centre line), the first and third quartiles (box limits) and 1.5× interquartile range (whiskers). d, Spatial distribution and UMAP of all clusters for ATAC, RNA and joint clustering of ATAC and RNA data. Overlay of clusters with the tissue image shows that spatial clusters precisely match anatomic regions. Pixel size, 50 µm; scale bars, 1 mm. e, Spatial mapping of GAS and gene expression for selected marker genes in different clusters for ATAC and RNA in spatial ATAC-RNA-seq. f, Pseudotime analysis from radial glia to postmitotic premature neurons visualized at the spatial level. g, Heatmaps delineating gene expression and GAS for marker genes. h, Dynamic changes in GAS and gene expression across pseudotime. fragments were retained in the supernatant. The libraries of gDNA and cDNA were constructed separately for next-generation sequencing (NGS). Spatial CUT&Tag-RNA-seq was performed by applying an antibody against a specific histone modification (H3K27me3, H3K27ac or H3K4me3) to the tissue section and then using a protein A-tethered Tn5-DNA complex to perform CUT&Tag. The remaining steps were similar to those for spatial ATAC-RNA-seq ( Fig. 1a and Extended Data Fig. 1a,d,e), but resulted in spatial co-profiling of genome-wide histone modification occupancy and transcriptome.

Spatial comapping of mouse embryo
Spatial ATAC-RNA-seq on E13 mouse embryo identified eight major ATAC clusters and 14 RNA clusters, suggesting that, at this stage of development, chromatin accessibility may not allow discrimination of all cell types defined by transcriptional profiles. Spatial distribution of these clusters agrees with tissue histology (see haematoxylin and eosin staining of an adjacent tissue section in Fig. 1d and Extended Data Fig. 2a). In the ATAC data, cluster A3 represents the embryonic eye field with open chromatin accessibility at the loci of genes including Six6 (Fig. 1d (left) and Fig. 1e). Clusters A4-A5 are associated with several developing internal organs. Clusters A6-A7 cover the central nervous system (CNS). To benchmark the ATAC result, we aggregated the chromatin accessibility profiles in pixels within specific organs and compared these with organ-specific ENCODE E13.5 ATAC-seq reference data ( Supplementary Fig. 3b). Additionally, the peaks obtained from our spatial ATAC data are also consistent with the ENCODE reference (Supplementary Fig. 3c). We integrated spatial ATAC data with single-cell RNA-seq 12 data for cell type assignment in each cluster (Extended Data Fig. 2c,d). As expected, radial glia (neural stem/progenitor cells) were observed predominantly in the ventricular layer, and differentiated cell types such as postmitotic premature neurons and inhibitory interneurons were enriched in those regions distant from the ventricular layer (Extended Data Fig. 2d).
Cell type-specific marker genes were identified for individual clusters, and their expression was inferred from chromatin accessibility ( Fig. 1e and Extended Data Fig. 2b,e) and predicted by gene activity score (GAS 13 ; Methods). Sox2, which is involved in the development of nervous tissue and optic nerve formation 14,15 , showed high chromatin accessibility in the embryonic eye field and in the ventricular layer containing neural stem/progenitor cells. Pax6 exhibited a similar spatial pattern of chromatin accessibility. Myt1l, which encodes myelin transcription factor 1-like protein, presented a higher ATAC signal in the embryonic brain and neural tube 16 . Six6, a key gene involved in eye development 17 , showed highest GAS in the eye region. Nrxn2, which encodes Neurexin 2, a key gene in the vertebrate nervous system 18 , was extensively accessible in most neural cell regions. Accessibility at Rbfox3, which encodes RNA binding protein fox-1 homolog 3 (a splicing factor known as NeuN) 19 , was observed in neurons 19 whereas Sox1 and Sox2 presented enriched chromatin accessibility in the ventricular layer ( Fig. 1e and Extended Data Fig. 2b). Cell type-specific enrichment of transcription factor (TF) regulators was also examined using Chrom-VAR 20 analysis of deviation in TF motifs (Sox2 and Nfix), and identified the positive TF regulators (Extended Data Fig. 2g,h). We observed that the Sox2 motif was enriched in cluster A7, consistent with its function in embryonic brain development 21 . GREAT analysis 22 further verified the strong concordance between gene regulatory pathways and anatomical annotation (Extended Data Fig. 2i,j).
For the RNA spatial data, 14 distinct clusters were identified and characterized by specific marker genes ( Fig. 1d (middle) and Extended Data Fig. 2f). For example, cluster R10 (Six6) was correlated with the embryonic eye, clusters R2, R6 and R8 were related to the CNS and cluster R7 was associated with the formation of cartilage. To evaluate data quality we performed correlation analysis with organ-specific ENCODE E13.5 RNA-seq reference data, which showed high concordance ( Supplementary Fig. 3a). Cell type-specific marker genes were also identified for individual RNA clusters-for example, Mapt (R2) may function in establishing and maintaining neuronal polarity 23 . Epha5 (R6) is involved in axon guidance, whereas Slc1a3 (R8) is involved in regulation of excitatory neurotransmission in the CNS 18 . Myh3 (R3) is associated with muscle contraction 24 . Col2a1 (R7), which shows specific expression in cartilaginous tissues, has an essential role in normal embryonic skeleton development 18 . Pathway analysis 25 results agree with anatomical annotation (Supplementary Fig. 4a).
Integration of spatial RNA data with scRNA-seq mouse organogenesis data 12 was performed to determine cell identities in each pixel (Extended Data Fig. 2c,d). We observed that radial glia, postmitotic premature neurons and inhibitory interneurons (clusters R8, R6 and R2) were present in the same major clusters as shown in the ATAC analysis (clusters A7 and A6), which verified the use of multiple omics information for more robust cell type identification. We attempted joint clustering of spatial ATAC and RNA data to refine the spatial patterns. A new neuronal cluster ( J10) was identified in the joint clustering analysis, which was not readily resolved by single modalities alone (Fig. 1d, right). This result highlights the value of using joint multiomics profiles to improve spatial cell type mapping 26 .
Co-profiling of the epigenome and transcriptome facilitates investigation of the correlation between accessible peaks and expressed genes pixel by pixel in the tissue context. We observed distinct signals at predicted enhancers for some genes (Sox2, Pax6 and Sox1) (Extended Article Data Fig. 2e and Supplementary Fig. 4b). For example, enhancers for Sox2 and Pax6 had higher chromatin accessibility in clusters A7 and A3, respectively, suggesting their roles in these tissue regions in the regulation of Sox2 and Pax6 expression. Although spatial RNA distribution of these genes corresponded with their chromatin accessibility, the expression level may differ significantly from the degree of accessibility ( Fig. 1e and Extended Data Fig. 2b). Some marker genes (Pax6, Sox2 and Myt1l) were highly accessible in some regions of the embryonic brain but showed modest or low RNA expression ( Fig. 1e and Extended Data Fig. 2b), which may indicate lineage priming 10,27 of these genes in embryonic brain development. Despite a strong correlation between replicates ( Supplementary Fig. 2d), this observation could be due in part to the sequencing depth and RNA detection sensitivity. Nevertheless, these results still highlight the potential to link genome-wide epigenetic regulation to transcription in the spatial tissue context.
To investigate the spatiotemporal relationship between chromatin accessibility and gene expression in embryonic development, we analysed the differentiation trajectory from radial glia to various types of neurons such as postmitotic premature neurons 28 . Pseudotime analysis 29 was conducted under the ATAC pseudotime coordinate system, and developmental trajectories were directly visualized in the spatial tissue map (Fig. 1f). Chromatin accessibility GAS and gene expression along this trajectory showed dynamic changes in selected marker genes (Fig. 1g,h). As expected, the expression levels of Sox2, Pax6 and other genes involved in progenitor maintenance and proliferation (Gene Ontology (GO) Fabp7 to Pax6 in Extended Data Fig. 3a) were downregulated during the transition to postmitotic neurons. The loss of chromatin accessibility at the Pax6 and radial glia marker Fabp7 loci preceded downregulation of the corresponding RNAs (Fig. 1h). In turn, genes involved in neuronal identity, axonogenesis and synapse organization (GO Myt1l to Dnm3 in Extended Data Fig. 3a), such as Dcx and Tubb3, showed increased expression in spatial pseudotime but chromatin accessibility at their loci was already elevated at earlier stages, suggesting lineage priming of these genes in expression 10,27 . We also found a cohort of genes whose expression rapidly declined during spatial pseudotime but whose chromatin accessibility was maintained throughout pseudotime or declined only at very late stages. Many of these genes, such as Ptprz1, Bcan and Luzp2, are characteristic of oligodendrocyte precursor cells, and the corresponding GO biological processes are negative regulation of myelination and regulation of gliogenesis (GO Cp to Sparcl1 in Extended Data Fig. 3a), suggesting that the neuronal lineage might retain the potential to acquire an oligodendrocyte identity even when it has already migrated away from the ventricular zone in the embryonic brain. Monocle2 pseudotime analysis showed a bifurcation in chromatin accessibility but not in RNA expression; one path led to regions close to the ventricular zone (green pixels, Extended Data Fig. 3d,e) and the other terminated in regions distal to the ventricular zone (blue pixels). In contrast to the green path, the blue path presented an increase in chromatin accessibility for genes involved in axonogenesis and dendrite formation (Extended Data Fig. 3e (red box), 3f), suggesting that the chromatin state of neural cells distal to the ventricle corresponds to a more differentiated neuronal state. Thus, our spatial ATAC-RNA-seq can be used to decipher the gene regulation mechanism and spatiotemporal dynamics during tissue development.

Spatial ATAC-RNA-seq of mouse brain
We developed a microchip with 100 serpentine microchannels to barcode 100 × 100 or, in total, 10,000 pixels per tissue section and up to five samples per chip (Fig. 2a, 20-µm pixel size). This was applied to P22 mouse brain coronal sections (at bregma 1) for joint profiling of chromatin accessibility with transcriptome. In contrast to E13 mouse embryonic brain, we found a higher number of ATAC clusters (14) compared with RNA clusters (11), indicating terminal differentiation of most major cell types in the juvenile brain relative to that in embryo, which consists of many undifferentiated or multipotent cell states. Spatial distribution of clusters agreed with the anatomical annotation defined by Nissl staining, reflecting arealization of the juvenile brain (Fig. 2b). Moreover, spatial clusters between ATAC and RNA showed concordance in cluster assignment ( Supplementary Fig. 5a). Chromatin accessibility of marker genes identified major regions such as striatum (Pde10a and Adcy5, markers of medium spiny neurons, cluster A1), corpus callosum (Sox10, Mbp and Tspan2, markers of oligodendroglia, cluster A3), cortex (Mef2c, Neurod6 and Nrn1, clusters A0 and A4, markers of excitatory neurons; Cux2, cluster A4) and lateral ventricle (Dlx1, Pax6, Notch1 and Sox2, markers of ependymal/neural progenitor cells, cluster A11) ( Fig. 2e and Extended Data Fig. 4a,b). GREAT analysis confirmed that the major pathways correlated with tissue functions in different anatomical regions ( Supplementary Fig. 5b,c). Spatial RNA clusters also showed region-specific gene expression such as striatum (Pde10a, cluster R2, medium spiny neurons), corpus callosum (Mbp and Tspan2, cluster R5, oligodendroglia) and cortex (Mef2c, clusters R0 and R1, excitatory neurons) ( Fig. 2e and Extended Data Fig. 4a).
Integration of single-cell ATAC-seq mouse brain atlas data 30 with spatial ATAC-seq data facilitated identification of all major cell types, and label transfer was then used to assign cell types to spatial locations where the epigenetic state may control specific cell type formation ( Fig. 2c and Extended Data Fig. 4c). For example, we observed immature oligodendrocytes and oligodendrocytes in cluster A3 corresponding to corpus callosum. A thin layer of radial glia-like cells was found in the lateral ventricle, medium spiny neurons (D1MSN and D2MSN) in striatum and inhibitory neurons in the lateral septal nucleus 30 (Extended Data Fig. 4c) based on chromatin accessibility. Furthermore, the subclasses of excitatory neurons in different layers of cortex (ITL23GL, Cortex L2/3; ITL4GL, Cortex L4; ITL5GL, Cortex L5; PTGL, Cortex PT; NPGL, Cortex NP; ITL6GL, Cortex L6; CTGL, Cortex CT; and L6bGL, Cortex L6b) were clearly identified 31 (Extended Data Fig. 4c). Integration of spatial RNA data and the juvenile CNS scRNA-seq atlas 32 allowed for visualization of dominant transcriptional cell types in tissue ( Fig. 2d and Extended Data Fig. 4d). We observed a high degree of concordance in spatial cell type distribution as determined by ATAC versus RNA, suggesting a general congruence between chromatin accessibility and transcriptome in defining cell identities in tissue. For instance, neuronal intermediate progenitor cells (SZNBL, included in radial glial-like cells from scATAC-seq 30 ), mature oligodendrocytes (MOL, included in oligodendrocytes from scATAC-seq), newly formed oligodendrocytes (included in immature oligodendrocytes from scATAC-seq) and medium spiny neurons (MSN) shown by the RNA data were enriched in the same regions as identified by the ATAC data (Extended Data Fig. 4c,d). Furthermore, the detection of a thin layer of SZNBL in the lateral ventricle and vascular leptomeningeal cells (VLMC) in the regions with preserved meninges demonstrated a high spatial resolution for our technology in identifying low-abundance cells, even at near-single-cell resolution (for VLMC) (Extended Data Fig. 4c,d).
Joint profiling of ATAC and RNA facilitates inferring the gene regulatory landscape by searching correlated peak accessibility and gene expression. We detected 21,417 significant peak-to-gene linkages between regulatory elements and target genes (Extended Data Fig. 5c,d). Some potential enhancers with dynamically regulated promoter interactions were found to be enriched in specific clusters such as Dlx1 and Sox2 (cluster A11), Tspan2 (cluster A3) and Adcy5 (cluster A1) (Extended Data Fig. 4b), indicating the ability of spatial ATAC-RNA-seq to identify the key regulatory regions for target genes. Spatial patterns of TF motif enrichment (Dnajc21, Pax6, Sox2, Sox4 and Mef2c) provided further insights into the TF regulators 33 visualized in tissue (Extended Data Fig. 5a,b). For example, Sox2 was examined simultaneously for chromatin accessibility, gene expression, putative enhancers and TF motif enrichment, enabling a more comprehensive understanding of gene regulation dynamics. As observed in the embryonic mouse brain, some marker genes with open chromatin accessibility (Sox10, Sox2, Neurod6, Pax6 and Notch1) were either not or lowly expressed ( Fig. 2e and Extended Data Fig. 4a). A biological replicate experiment on P21 mouse brain was also performed (Extended Data Fig. 6; 20 µm pixel size, 50 × 50 barcodes). Many of these genes encode for transcription factors involved in neural development, suggesting the possibility of epigenetic-but not transcriptional-memory of these genes in brain development.

Spatial CUT&Tag-RNA-seq of mouse brain
In addition to chromatin accessibility, histone modifications are also a key aspect of epigenetic regulation. Spatial CUT&Tag-RNA-seq was performed to co-profile transcriptome and H3K27me3 (repressing loci), H3K27ac (activating promoters and/or enhancers) or H3K4me3 (active promoters) histone modifications, respectively, in P22 mouse brain (20-µm pixel size, 100 × 100 barcodes). We identified 13 and 15 specific clusters for H3K27me3 and RNA, 12 and 13 clusters for H3K27ac and RNA and 11 and 12 clusters for H3K4me3 and RNA, respectively (Fig. 3a-c). These clusters agreed with the anatomical annotation by Nissl staining and showed good concordance between CUT&Tag and RNA in spatial patterns (Fig. 3a-c), which was further confirmed by Belayer 34 analysis (Supplementary Fig. 7c).
Integration and co-embedding of spatial CUT&Tag and scCUT&Tag data 36,37 (Fig. 3d,e and Extended Data Fig. 7a) showed that the epigenetic states in our data for H3K27me3, H3K27ac and H3K4me3 agreed with the corresponding projection in scCUT&Tag. Integration with the corresponding juvenile CNS scRNA-seq atlas 32 allowed for label transfer to assign transcriptional cell types to the spatial location of epigenetic identities/states (Extended Data Fig. 7g,h, bottom). For example, we observed an enrichment of MOL within the corpus callosum, a thin layer of ependymal cells in the lateral ventricle, excitatory neurons in cerebral cortex and MSN in the striatum. We also integrated paired RNA data with the scRNA-seq atlas 32 to identify dominant transcriptional cell types via label transfer (Fig. 3f and Extended Data Fig. 7b-h). MOL, a thin layer of ependymal cells, MSN and excitatory neurons were enriched in the same spatial regions identified by CUT&Tag (H3K27ac and H3K4me3) (Extended Data Fig. 7g,h, top). In particular for spatial H3K27me3, whereas integration with scCUT&Tag could not clearly show the identity of several clusters in the epigenomic modalities (clusters C0, C1 and C3), label transfer of scRNA-seq data with paired RNA data from the same tissue section clearly showed cell identities in these clusters (Fig. 3a,d and Extended Data Fig. 7f), highlighting the power of combining CUT&Tag and RNA-seq (spatial CUT&Tag-RNA-seq) in the same tissue section to more accurately identify cell types or states. To directly infer interactions between genome-wide gene expression and the corresponding enhancers across all clusters, we identified a total of 19,468 significant peak-to-gene linkages from spatial CUT&Tag (H3K27ac)-RNA-seq ( Supplementary Fig. 7a,b). A biological replicate was also performed on P21 mouse brain (Extended Data Fig. 8; 20-µm pixel size, 50 × 50 barcodes).

Region-specific gene expression regulation
To further understand the spatial epigenetic regulation of gene expression at the genome scale, we compared the GAS or CSS obtained from ATAC and CUT&Tag (H3K27me3, H3K27ac and H3K4me3) with the corresponding RNA expression in all major tissue regions of P22 mouse brain (Fig. 4a-c and Supplementary Figs. 9a, 10a and 11a). For example, in oligodendrocyte-abundant corpus callosum we observed a robust anticorrelation between H3K27me3 and RNA (Fig. 4a)  Pixel size, 20 µm; scale bars, 1 mm. d, Integration of H3K27me3 data with scCUT&Tag (H3K27me3) data 37 from mouse brain. e, Integration of H3K27ac data with scCUT&Tag (H3K27ac) data 37 from mouse brain. f, Integration of RNA data in spatial CUT&Tag (H3K27me3)-RNA-seq, spatial CUT&Tag (H3K27ac)-RNA-seq and spatial CUT&Tag (H3K4me3)-RNA-seq with scRNA-seq data 32 from mouse brain.
to GO biological processes such as myelination and regulation of oligodendrocyte differentiation (Extended Data Fig. 9a-c and Supplementary Fig. 8; GO for quadrant IV). By contrast, genes such as Grin2b and Syt1 with high CSS and low RNA expression are related to neuronal processes such as synaptic transmission and neurotransmitter release (Extended Data Fig. 9d,e and Supplementary Fig. 8; GO for quadrant II). These results are in accordance with ATAC/H3K27me3/ H3K27ac Nano-CUT&Tag analysis 37 of oligodendroglial differentiation in the juvenile brain, which showed two waves of H3K27me3 repression during this process 37 . We also found a small subset of genes ( Fig. 4a and Supplementary Fig. 8; GO for quadrant I) with higher levels of both H3K27me3 and RNA in the corpus callosum, such as Ptprz1, a marker of oligodendrocyte precursor cells 38,39 (Extended Data Fig. 9g), which might indicate transcriptional poising of some of these genes. This could also be due in part to the presence in proximity of both oligodendrocyte precursor cells (expressing Ptprz1) and mature oligodendrocytes (where Ptprz1 is repressed) in the corpus callosum. We performed a similar correlation analysis for RNA and H3K27me3 in other regions of the juvenile brain, including the striatum and superficial and deeper cortical layers (Supplementary Figs. 9a, 10a  and 11a). We also observed strong anticorrelation with a cohort of genes involved in neuronal processes being activated or repressed in a region-specific manner. For instance, in the superficial cortical layer, genes for GABAergic regulation of synaptic transmission were enriched in H3K27me3 whereas glutamatergic synapse transmission genes had high expression and low H3K27me3 (Supplementary Figs. 9b, 10b and 11b; GOs for striatum, superficial and deeper cortical layers). In contrast to corpus callosum, only a limited number of genes positive for both H3K27me3 and RNA were found in these regions. Despite a general anticorrelation between RNA expression and H3K27me3 deposition, we also observed regional differences (blue-shaded columns in Fig. 4d); Nav3 and Sncb, with low levels of H3K27me3 in both the cortex and striatum, were expressed in the former but not the latter (Extended Data Fig. 9h). The opposite pattern was observed for genes such as Ablim2 and Gng7 (Extended Data Fig. 9h). Car2, a marker gene expressed in oligodendrocytes and abundant in the corpus callosum, presented H3K27me3 occupancy in the cortical layers, in contrast to the corpus callosum and, surprisingly, the striatum (Extended Data Fig. 9c). Thus, mechanisms other than Polycomb-mediated H3K27me3 deposition might be involved in the transcriptional repression of these genes in different areas of the CNS.
In contrast to H3K27me3, we observed a robust correlation between RNA and activating marks H3K4me3 or H3K27ac in the corpus callosum, with genes highly expressed and with high deposition in these two modalities regulating processes in oligodendroglia (Fig. 4b,c and Supplementary Fig. 12; GO for H3K27ac quadrant I and H3K4me3 quadrant I). For instance, Mal showed the highest activity or expression in the corpus callosum for ATAC, H3K27ac, H3K4me3 and paired RNA (Extended Data Fig. 9a). Gpr88, a marker gene of MSN, was enriched in striatum for ATAC, H3K27ac and H3K4me3 but with low CSS in striatum for H3K27me3 (Extended Data Fig. 9f). Furthermore, we examined collective regulation among different histone modifications in the Article corpus callosum ( Fig. 4e and Supplementary Figs. 13 and 14). In general H3K4me3 and H3K27ac showed high correlations, with their signals anticorrelated with H3K27me3 (Fig. 4e). We also observed H3K4me3 or H3K27ac and H3K27me3 co-occupancy in several gene loci in this tissue region ( Fig. 4e and Supplementary Figs. 13 and 14). These genes are involved in oligodendrocyte differentiation and myelination and other processes such as protein catabolism and localization (Supplementary Figs. 13 and 14). Some of these genes, such as Ptprz1 and Fnbp1, also presented higher levels of H3K27me3 and RNA (Fig. 4a). As mentioned above, this potential H3K4me3/H3K27me3 bivalency may reflect transcriptional poising.

Spatial comapping of human brain
We performed spatial ATAC-RNA-seq on adult (31-year-old male) human brain hippocampal formation, which is a complex brain region involved in cognitive functions and diseases such as major depression disease 40 and Alzheimer's disease 41 . We identified seven and eight major clusters for ATAC and RNA, respectively, and the spatial pattern agreed with major anatomical landmarks (Fig. 5a) 3 , including specific marker genes (Fig. 5d). The ATAC cluster A4 represents the granule cell layer (GCL) (THY1, BCL11B) and cluster A6 corresponds to the choroid plexus (Fig. 5a,d). TF motifs (NEUROD1 and SNAI1) and their spatial patterns were visualized in different tissue regions, and positive TF regulators were identified (Extended Data Fig. 10b,c). For the RNA data, we also detected distinct clusters with unique marker genes (Fig. 5a,d) such as PROX1 and BCL11B enriched in cluster R4 (GCL). We also conducted integration of our ATAC data and scATAC-seq data from human brain samples 42 (Fig. 5b), and our RNA data with adult human brain snRNA-seq data 43 , to show the dominant cell identities and states. Cell types identified by snRNA-seq 43 were assigned to each cluster via label transfer ( Fig. 5c and Extended Data Fig. 10a). Granule cells (GC) were detected in the GCL, cornu ammonis (CA) neurons were enriched in CA3-4 regions and VLMC were strongly distinguished in choroid plexus as opposed to other regions.
In general, both ATAC and RNA can readily resolve all major tissue features in this region, but spatial cosequencing of epigenome and transcriptome may provide new insights into the dynamic gene regulation mechanism that cannot be realized by single modalities. For example, PROX1, a signature gene defining granule neuron identity during pyramidal neuron fate selection 44 , is indeed highly expressed in GCL but showed modest chromatin accessibility (Fig. 5d). This might be attributed to a minimal demand for synthesis of new PROX1 transcripts in postmitotic mature granule cells, and is thus not required to maintain an active open chromatin state.

Discussion
Spatial omics technologies (spatial epigenomics, transcriptomics and proteomics), based on either NGS [1][2][3][4][45][46][47] or imaging 5,48 , offer an unprecedented opportunity to generate new and rich insights into gene regulation in the spatial tissue context. However, to comprehensively understand the mechanism of gene regulation, different layers of omics information need to be profiled simultaneously. This was first demonstrated in dissociated single cells [6][7][8][9] but is yet to be realized directly in tissue. Imaging-based DNA sequential fluorescence in situ hybridization combined with RNA sequential fluorescence in situ hybridization detected spatial chromatin and gene expression for ATAC pixels in spatial ATAC-RNA-seq (this work) Single-nucleus RNA-seq RNA pixels in spatial ATAC-RNA-seq (this work) Oligo.
DG GC CA Astro.
Spatial ATAC-RNA-seq (RNA) c, Integration of our RNA data with snRNA-seq data from human brain 43 . d, Spatial mapping of GAS and gene expression for selected marker genes in different clusters for ATAC and RNA. Oligo, oligodendrocytes; astro, astrocytes; DG, dentate gyrus. target genes and selected genomic loci 49 . As of today, current technologies have not been able to perform unbiased genome-wide comapping of epigenome and transcriptome on the same tissue section at the cellular level. We developed spatial ATAC-RNA-seq and spatial CUT&Tag-RNA-seq (applied to H3K27me3, H3K27ac and H3K4me3) for co-profiling of genome-wide chromatin accessibility or histone modifications in conjunction with whole transcriptome on the same tissue section at near-single-cell resolution (available for spatial browsing, see 'Data availability'). These techniques were applied to comapping of embryonic and juvenile mouse brain, as well as adult human brain. Spatial epigenome-transcriptome cosequencing yielded two layers of spatial omics information directly in tissue, with data quality similar to that previously obtained by single modalities. These technologies allow us to examine spatiotemporal dynamics and genome-wide gene regulation mechanisms in the tissue context. We have demonstrated either ATAC or CUT&Tag in conjunction with RNA for spatial multiomics mapping. It might be possible to combine all three-chromatin accessibility, histone modifications and transcriptome-to delineate a more comprehensive landscape of gene regulation network in tissue. It might be also possible to simultaneously measure multiple histone modifications to assess the multivalency effect on spatial gene expression regulation. Previously we demonstrated spatial profiling of transcriptome and a large panel of proteins 1 . We envision that it is possible and highly valuable to further combine epigenome, transcriptome and proteome 50 to dissect spatial patterns of cell type, state and gene regulation. Finally, spatial multiomics as reported herein may find widespread applications beyond neuroscience and developmental biology. For example, multiple modalities for spatial mapping of human disease tissues not only cross-validate one another but also better elucidate the mechanisms driving abnormal cell states that could not be readily discerned using single-modality methods. In addition, spatial information may further show how the local tissue environment affects cell state, dynamics and function across all layers of the central dogma. In this work, we developed a device with a larger mapping area and higher throughput that allows the mapping of a fourfold larger tissue area, using 100 × 100 barcodes (total 10,000 pixels, 20-µm pixel size) to cover almost the entire hemisphere of a juvenile P22 mouse brain coronal section. The serpentine channel design allows simultaneous processing of five tissue samples, each mapped for 10,000 pixels. It is feasible to further increase the mapping area by increasing the number of barcodes and further enhance throughput using automated liquid handling.
In summary, spatially resolved, genome-wide cosequencing of epigenome and transcriptome at the cellular level represents one of the most informative tools in spatial biology and can be applied to a wide range of biological and biomedical research.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-023-05795-1. Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Juvenile mouse brain tissue (P21−P22) was obtained from the Sox10:Cre-RCE:LoxP enhanced green fluorescent protein (eGFP) line on a C57BL/6xCD1 mixed genetic background maintained at Karolinska Institutet. The line was first generated by crossing Sox10:Cre animals (The Jackson Laboratory, no. 025807) with RCE:loxP (eGFP) animals (The Jackson Laboratory, no. 032037-JAX). This was established and maintained by breeding males lacking the Cre allele with females carrying a hemizygous Cre allele; the reporter allele was maintained in homozygosity or hemizygosity in both males and females. This resulted in specific labelling of oligodendrocyte lineage with eGFP. All animals were free from mouse bacterial and viral pathogens, ectoparasites and endoparasites. The following light/dark cycle was maintained for the mice: dawn, 06:00-07:00; daylight, 07:00-18:00; dusk, 18:00-19:00; night, 19:00-06:00. Mice were housed in individually ventilated cages at a maximum number of five per cage (IVC sealsafe GM500, tecniplast). General housing parameters including temperature, ventilation and relative humidity followed the European Convention for the Protection of Vertebrate Animals used for experimental and other scientific purposes. Air quality was controlled using stand-alone air-handling units equipped with a high-efficiency particulate air filter. Relative air humidity was consistently 55 ± 10%, at a temperature of 22 °C. Husbandry parameters were monitored with ScanClime (Scanbur) units. The cages contained a card box shelter, gnawing sticks and nesting material (Scanbur), placed on hardwood bedding (TAPVEI). Mice were provided with a regular chow diet, and water was supplied by water bottle and changed weekly. Cages were changed every 2 weeks in a laminar air-flow cabinet.
All experimental procedures were conducted following European directive no. 2010/63/EU, local Swedish directive no. L150/SJVFS/2019:9, Saknr no. L150 and Karolinska Institutet complementary guidelines for procurement and use of laboratory animals, no. Dnr. 1937/03-640. All procedures described were approved by the local committee for ethical experiments on laboratory animals in Sweden (Stockholms Norra Djurförsöksetiska nämnd, nos. 1995/2019 and 7029/2020). Human brain tissue was obtained from the brain collection of the New York State Psychiatric Institute at Columbia University, which includes brain samples from the Republic of Macedonia. Brain tissue collection was conducted with New York State Psychiatric Institute Institutional Review Board approval, and informed consent obtained from next of kin who agreed to donate the brains and participate in psychological autopsy interviews.
We analysed brain hippocampus tissue from a 31-year-old Caucasian male individual, with no psychiatric or neurological diagnosis, who had died of a traumatic accident and had a high level of global functioning before death as measured by global assessment scale 51 score, which was 90 (scoring 1-100, with 100 the highest function), and with toxicology negative for psychotropic medications and drugs. Postmortem interval (time from demise to brain collection) was 6.5 h.
The anterior hippocampal region was dissected from a fresh-frozen coronal section (20-mm thickness) of the right brain hemisphere. The dentate gyrus region (around 10 × 10 mm 2 ) of the anterior hippocampal region was selected. Cryosections of 10 µm were collected on poly-l-lysine-coated glass slides (no. 63478-AS, Electron Microscopy Sciences). Samples were stored at −80 °C until further use.

Fabrication of PDMS microfluidic device
Chrome photomasks were purchased from Front Range Photomasks, with a channel width of either 20 or 50 µm. The moulds for polydimethylsiloxane (PDMS) microfluidic devices were fabricated using standard photolithography. The manufacturer's guidelines were followed to spin-coat SU-8-negative photoresist (nos. SU-2025 and SU-2010, Microchem) onto a silicon wafer (no. C04004, WaferPro). The heights of the features were about 20 and 50 µm for 20-and 50-µm-wide devices, respectively. PDMS microfluidic devices were fabricated using the SU-8 moulds. We mixed the curing and base agents in a 1:10 ratio and poured the mixture into the moulds. After degassing for 30 min the mixture was cured at 70 °C for 2 h. Solidified PDMS was extracted for further use. We have published a detailed protocol for the fabrication and preparation of the PDMS device 52 . For RT, the following mixture was used: 12.5 µl of 5× RT buffer, 4.5 µl of RNase-free water, 0.4 µl of RNase inhibitor, 0.8 µl of Superase In RNase inhibitor, 3.1 µl of 10 mM deoxynucleotide triphosphate each, 6.2 µl of Maxima H Minus Reverse Transcriptase, 25 µl of 0.5× PBS-RI and 10 µl of RT primer. Tissues were incubated for 30 min at room temperature, then at 42°C for 90 min in a wet box. After the RT reaction, tissues were washed with 1× NEBuffer 3.1 and 1% RNase inhibitor for 5 min.

Spatial ATAC-RNA-seq
For first barcode (barcode A) in situ ligation, the first PDMS chip was covered to the tissue region of interest. For alignment purposes, a 10× objective (Thermo Fisher Scientific, EVOS FL Auto microscope no. AMF7000, EVOS FL Auto 2 Software Revision 2.0.2094.0) was used to take a brightfield image. The PDMS device and tissue slide were clamped tightly with a home-made acrylic clamp. Barcode A was first annealed with ligation linker 1, 10 µl of 100 µM ligation linker, 10 µl of 100 µM each barcode A and 20 µl of 2× annealing buffer (20 mM Tris pH 7.5-8.0, 100 mM NaCl, 2 mM EDTA) and mixed well. For each channel, 5 µl of ligation master mixture was prepared with 2 µl of ligation mix (27 µl of T4 DNA ligase buffer, 72.4 µl of RNase-free water, 5.4 µl of 5% Triton X-100, 11 µl of T4 DNA ligase), 2 µl of 1× NEBuffer 3.1 and 1 µl of each annealed DNA barcode A (A1−A50/A100, 25 µM). Vacuum was used to load the ligation master mixture into 50 channels of the device, followed by incubation at 37 °C for 30 min in a wet box. The PDMS chip and clamp were removed after washing with 1× NEBuffer 3.1 for 5 min. The slide was washed with water and dried in air.
For second barcode (barcode B) in situ ligation, the second PDMS chip was covered to the slide and a further brightfield image taken with the 10× objective. An acrylic clamp was applied to clamp the PDMS and tissue slide together. Annealing of barcodes B (B1−B50/ B100, 25 µM) and preparation of the ligation master mix were carried out as for barcodes A. The device was incubated at 37 °C for 30 min in a wet box. The PDMS chip and clamp were removed after washing with 1× DPBS with SUPERase In RNase inhibitor for 5 min. The slide was washed with water and dried in air. A brightfield image was then taken for further alignment.
For tissue lysis, the region of interest was digested with 100 µl of reverse crosslinking mixture (0.4 mg ml -1 proteinase K, 1 mM EDTA, 50 mM Tris-HCl pH 8.0, 200 mM NaCl, 1% SDS) at 58 °C for 2 h in a wet box. The lysate was collected in a 1.5 m tube and incubated at 65 °C overnight.
For DNA and cDNA separation, the lysate was purified with Zymo DNA Clean & Concentrator-5 and eluted to 100 µl of RNase-free water. The 1× B&W buffer with 0.05% Tween-20 was used to wash 40 µl of Dynabeads MyOne Streptavidin C1 beads three times. Then, 100 µl of 2× B&W buffer with 2.5 µl of SUPERase In RNase inhibitor was used to resuspend the beads, which were then mixed with the lysate and allowed to bind at room temperature for 1 h with agitation. A magnet was used to separate beads and supernatant in the lysate.
The supernatant was removed for ATAC library construction, then purified with Zymo DNA Clean & Concentrator-5 and eluted to 20 µl of RNase-free water. PCR solution (25 µl of 2× NEBNext Master Mix, 2.5 µl of 25 µM indexed i7 primer, 2.5 µl of 25 µM P5 PCR primer) was added and mixed well. PCR was first performed with the following programme: 72 °C for 5 min, 98 °C for 30 s and cycling at 98 °C for 10 s, 63 °C for 30 s and 72 °C for 1 min, five times. To determine additional cycles, the pre-amplified mixture (5 µl) was mixed with quantitative PCR (qPCR) solution (5 µl of 2× NEBNext Master Mix, 0.24 µl of 25× SYBR Green, 0.5 µl of 25 µM new P5 PCR primer, 3.76 µl of nuclease-free H 2 O, 0.5 µl of 25 µM indexed i7 primer). The qPCR reaction was then carried out with the following programme: 98 °C for 30 s with cycling at 98 °C for 10 s, 63 °C for 30 s and 72 °C for 1 min, 20 times. The remaining pre-amplified DNA (45 µl) was amplified by running additional cycles as determined by qPCR (to reach one-third of saturated signal). The final PCR product was purified by 1× Ampure XP beads (45 µl) and eluted in 20 µl of nuclease-free H 2 O.
The beads were used for cDNA library construction. They were first washed twice with 400 µl of 1× B&W buffer with 0.05% Tween-20 and once with 10 mM Tris pH 8.0 containing 0.1% Tween-20. Streptavidin beads with bound cDNA molecules were resuspended in a TSO solution (22 µl of 10 mM deoxynucleotide triphosphate each, 44 µl of 5× Maxima RT buffer, 44 µl of 20% Ficoll PM-400 solution, 88 µl of RNase-free water, 5.5 µl of 100 uM template switch primer (AAGCAGTGGTATCAACGCA GAGTGAATrGrG+G), 11 µl of Maxima H Minus Reverse Transcriptase, 5.5 µl of RNase Inhibitor). The beads were incubated at room temperature for 30 min and then at 42 °C for 90 min, with gentle shaking. After washing beads once with 400 µl of 10 mM Tris and 0.1% Tween-20 and once with water, they were resuspended in a PCR solution (110 µl of 2× Kapa HiFi HotStart Master Mix, 8.8 µl of 10 µM primers 1 and 2, 92.4 µl of RNase-free water). PCR thermocycling was carried out using the following programme: 95 °C for 3 min and cycling at 98 °C for 20 s, 65 °C for 45 s and 72 °C for 3 min, five times. After five cycles, beads were removed from the PCR solution and 25× SYBR Green was added at 1× concentration. Samples were again placed in a qPCR machine with the following thermocycling conditions: 95 °C for 3 min, cycling at 98 °C for 20 s, 65 °C for 20 s and 72 °C for 3 min, 15 times, followed by 5 min at 72 °C. The reaction was removed once the qPCR signal began to plateau. The PCR product was purified with 0.8× Ampure XP beads and eluted in 20 µl of nuclease-free H 2 O.
A Nextera XT Library Prep Kit was used for library preparation. Purified cDNA (1 ng) was diluted in RNase-free water to a total volume of 5 µl, then 10 µl of Tagment DNA buffer and 5 µl of Amplicon Tagment mix were added with incubation at 55 °C for 5 min; 5 µl of NT buffer was then added, with incubation at room temperature for 5 min. PCR master solution (15 µl of PCR master mix, 1 µl of 10 µM P5 primer, 1 µl of 10 µM indexed P7 primer, 8 µl of RNase-free water) was added and the PCR reaction performed with the following programme: 95 °C for 30 s, cycling at 95 °C for 10 s, 55 °C for 30 s, 72 °C for 30 s and 72 °C for 5 min, for 12 cycles. The PCR product was purified with 0.7× Ampure XP beads to obtain the library.
An Agilent Bioanalyzer High Sensitivity Chip was used to determine size distribution and concentration of the library before sequencing. NGS was conducted on an Illumina NovaSeq 6000 sequencer (paired-end, 150-base-pair mode).

Spatial CUT&Tag-RNA-seq
The frozen tissue slide was thawed for 10 min at room temperature. Tissue was fixed with formaldehyde (0.2%, with 0.05 U µl -1 RNase inhibitor) for 5 min and quenched with 1.25 M glycine for a further 5 min. After fixation, the tissue was washed twice with 1 ml of wash buffer (150 mM NaCl, 20 mM HEPES pH 7.5, one tablet of protease inhibitor cocktail, 0.5 mM Spermidine) and dipped in DI water. The tissue section was permeabilized with NP40-digitonin wash buffer (0.01% digitonin, 0.01% NP40 in wash buffer) for 5 min. The primary antibody (1:50 dilution with antibody buffer (0.001% BSA, 2 mM EDTA in NP40-digitonin wash buffer) was added with incubation at 4 °C overnight. The secondary antibody (guinea pig anti-rabbit IgG, 1:50 dilution with NP40-digitonin wash buffer) was added with incubation for 30 min at room temperature. The tissue was then washed with wash buffer for 5 min. A 1:100 dilution of pA-Tn5 adaptor complex in 300-wash buffer (one tablet of Protease inhibitor cocktail, 300 mM NaCl, 0.5 mM Spermidine, 20 mM HEPES pH 7.5) was added with incubation at room temperature for 1 h, followed by a 5 min wash with 300-wash buffer. Tagmentation buffer (10 mM MgCl 2 in 300-wash buffer) was added with incuation at 37 °C for 1 h. Next, 40 mM EDTA with 0.05 U µl -1 RNase inhibitor was added with incubation at room temperature for 5 min to stop tagmentation. Tissue was washed twice with 0.5× DPBS-RI for 5 min for further use.
For RT, two ligations and bead separation the protocols were as for spatial ATAC-RNA-seq. For construction of the CUT&Tag library, the supernatant was purified with Zymo DNA Clean & Concentrator-5 and eluted to 20 µl of RNase-free water. PCR solution (2 µl each of 10 µM P5 PCR primer and indexed i7 primer, 25 µl of NEBNext Master Mix) was added and mixed well. PCR was performed with the following programme: 58 °C for 5 min, incubation at 72 °C for 5 min and 98 °C for 30 s then cycling at 98 °C for 10 s, with incubation at 60 °C for 10 s 12 times and a final incubation at 72 °C for 1 min. The PCR product was purified by 1.3× Ampure XP beads using the standard protocol and eluted in 20 µl of nuclease-free H 2 O. cDNA library construction followed the spatial ATAC-RNA-seq protocol.
An Agilent Bioanalyzer High Sensitivity Chip was used to determine size distribution and concentration of the library before sequencing. NGS was conducted on an Illumina NovaSeq 6000 sequencer (paired-end, 150-base-pair mode).

Data preprocessing
For ATAC and CUT&Tag data, linkers 1 and 2 were used to filter read 2 and sequences were converted to Cell Ranger ATAC v.1.2 format (10X Genomics). Genome sequences were in the newly formed read 1, and barcodes A and B were included in newly formed read 2. Human reference (GRCh38) or mouse reference (GRCm38) was used to align fastq files. The BED-like fragments thus obtained were used to conduct downstream analysis. The fragments file includes fragments of information on spatial locations (barcode A × barcode B) and the genome.
For RNA data, read 2 was refined to extract barcode A, barcode B and UMI. ST pipeline v.1.7.2 (ref. 53 ) was used to map the processed read 1 against the mouse genome (GRCm38) or human genome (GRCh38), which created the gene matrix for downstream analysis that contains information on genes and spatial locations (barcode A × barcode B).

Data clustering and visualization
We first identified the location of pixels on tissue from the brightfield image using MATLAB 2020b (https://github.com/edicliuyang/Hip-lex_proteome).
Signac v.1.8 (ref. 54 ) was loaded in R v.4.1. ATAC, CUT&Tag and RNA matrices were read into Signac v.1.8 (ref. 54 ). The 'DefaultAssay' function was used for the RNA assay. For RNA data visualization, the feature was set to 3,000 with the 'FindVariableFeatures' function, then data were normalized using the 'SCTransform' function. Normalized RNA data were clustered and RNA UMAP was built. The DefaultAssay function was applied to the ATAC/CUT&Tag assay. For ATAC/CUT&Tag data visualization, minimum cutoff was set with the 'FindTopFeatures' function. Data were normalized and dimensionally reduced using latent semantic indexing, then ATAC/CUT&Tag data were clustered and ATAC/ CUT&Tag UMAP was built.
The DefaultAssay function was used for the joint ATAC/CUT&Tag and RNA assay. For visualization of joint ATAC/CUT&Tag and RNA data 55 , the 'FindMultiModalNeighbors' function was used. The reduction list was set to ('pca', 'lsi'), the dimensions list was set to that for RNA and ATAC/CUT&Tag, the modality weight.name was set to RNA weight and the joint UMAP was built.
To plot the above-generated UMAPs together, DefaultAssay was set to RNA and the UMAPs for ATAC/CUT&Tag, RNA or joint ATAC/CUT&Tag and RNA were visualized separately using 'DimPlot'.
In regard to RNA spatial data visualization, the gene matrix obtained from RNA was loaded into Seurat v.4.1 (ref. 56 ) as a Seurat object, and RNA metadata obtained from Signac were read into the Seurat object. All spatial maps were then plotted with the 'SpatialPlot' function.
In regard to ATAC/CUT&Tag spatial data visualization, the fragment file obtained from ATAC/CUT&Tag was read into ArchR v.1.0.1 (ref. 13 ) as an ArchRProject and the ATAC/CUT&Tag metadata obtained from Signac were read into the ArchRProject. The data from ArchRProject were normalized and dimensionally reduced using iterative latent semantic indexing. For GAS and CSS calculation we used the Gene Score model in ArchR. A gene score matrix was obtained for downstream analysis. The 'getMarkerFeatures' and 'getMarkers' functions in ArchR (testMethod = "Wilcoxon", cutOff = "FDR <= 0.05", groupBy = "seurat_ cluster") were used to find the marker genes/regions for each cluster. To visualize spatial data, results obtained from ArchR were input to Seurat v.4.1 to map the data back to the tissue. Pixel size was scaled using the 'pt.size.factor' parameter in the Seurat package for better visualization.
For peak-to-gene links we input RNA Seurat object using the 'addGe-neIntegrationMatrix' function in ArchR, then peak-to-gene links were drawn with the 'addPeak2GeneLinks' function. Co-accessibility of peaks was calculated using the 'addCoAccessibility' function in ArchR.

Integrative data analysis and cell type identification
Seurat v.4.1 (ref. 56 ) was used for RNA data integration and cell type identification, and the 'SCTransform' function to normalize our spatial RNA and scRNA-seq data. The 'SelectIntegrationFeatures' function was used to obtain features common to the two datasets. The 'Find-IntegrationAnchors' function was applied to find anchors, and the 'IntegrateData' function to create an integrated dataset through the identified anchors. The obtained integrated dataset was clustered, showing a good match between our spatial RNA and scRNA-seq data. The 'FindTransferAnchors' function was used to find transfer anchors, which were then used to conduct label transfer with the 'TransferData' function (if more than one cell type was presented in one pixel, the major cell type was assigned).
Signac v.1.8 and Seurat v.4.1 were used for integration of our ATAC/CUT&Tag and scATAC-seq/scCUT&Tag data. The scATAC-seq/ scCUT&Tag data were quantified according to our ATAC/CUT&Tag data to ensure that there were features common across both datasets. The FindIntegrationAnchors function (reduction = "rlsi") was used to identify anchors between the two datasets. The 'IntegrateEmbeddings' function was used to obtain an integrated dataset through the identified anchors. The obtained integrated dataset was clustered, showing a good match between our spatial ATAC/CUT&Tag and scATAC-seq/ scCUT&Tag data. For ATAC data, the FindTransferAnchors function was used to find transfer anchors, which were then used to map scATAC-seq to our spatial ATAC data with the 'MapQuery' function.
ArchR v.1.0.1 was used for cell type identification for our ATAC/ CUT&Tag data from scRNA-seq data. The gene score matrix of our ATAC/CUT&Tag was compared with the gene expression matrix from scRNA-seq, and aligned pixels from our ATAC/CUT&Tag data with cells from scRNA-seq. The function 'GeneIntegrationMatrix' was used to add pseudo-scRNA-seq profiles and cell identities.

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability
Raw and processed data reported in this paper are deposited in the Gene Expression Omnibus with accession code GSE205055. These datasets are available as web resources and can be browsed within the tissue spatial coordinates in the UCSC Cell and Genome Browser (https:// brain-spatial-omics.cells.ucsc.edu), and in our own data portal generated with AtlasXplore (https://web.atlasxomics.com/visualization/