Divergent combinations of enhancers encode spatial gene expression

Abstract

Spatial transcriptomics and epigenomics have enabled mapping gene regulation in the tissue context. However, it remains poorly understood how spatial gene expression patterns are orchestrated by enhancers. Here we build eSpatial, a computational framework that deciphers spatially resolved enhancer regulation of gene expression by integrating spatial profiles of gene expression and chromatin accessibility. Applying eSpatial to diverse spatial datasets, including mouse embryo and brain, as well as human melanoma and breast cancer, we reveal a "spatial enhancer code", in which divergent combinations of enhancers regulate the same gene in spatially segregated domains. We validate the spatial enhancer code using public spatial datasets such as VISTA, Allen in situ hybridization (ISH), and H3K27ac MERFISH. Moreover, we conduct transgenic reporter assays and in vivo CRISPR/Cas9-mediated perturbation experiments to confirm the Atoh1 spatial enhancer code in determining Atoh1 spatial expression in mouse embryonic spinal cord and brain. Our study establishes the spatial enhancer code concept, revealing how combinations of enhancers dynamically shape gene expression across diverse biological contexts, providing insights into tissue-specific regulatory mechanisms and tumor heterogeneity.

Fig. 1. Spatial domains and spatial-specific cis-regulatory elements in the P22 mouse brain.

a Schematic representation of the workflow in this study, including data collection, spatial enhancer code, and in vivo validation. b Spatial distribution of major cell types in spatial-ATAC-RNA-seq of the P22 mouse brain by integration of RNA pixels and scRNA-seq data. Oligo, oligodendrocyte; OPC, Oligodendrocyte precursor cell. Pixel size, 20 µm; scale bars, 1 mm. c Spatial distribution of RNA (c) and ATAC (d) clusters in spatial-ATAC-RNA-seq of the P22 mouse brain. Pixel size, 20 µm; scale bars, 1 mm. d Spatial distribution of identified spatial domains based on combined spatial profiles of chromatin accessibility and gene expression. Anatomical annotation defined from Nissl staining was obtained from ref. . Pixel size, 20 µm; scale bars, 1 mm. e Spatial mapping of denoised gene expression for selected marker genes across different spatial domains. f Heatmap showing gene expression (left) and enhancers chromatin accessibility (right) for 59,125 gene-enhancer pairs, clustered into 9 modules (K1–K9) based on gene expression. The relevant markers and representative enhancers were highlighted. The major cell types across spatial domains are shown in Supplementary Fig. 2b. g Spatial mapping of denoised chromatin accessibility for representative enhancers across different spatial domains. Source data are provided as a Source Data file.

Fig. 2. Spatial enhancer code in the P22 mouse brain.

a Bar plot illustrating the number of spatial domains in which genes are expressed. Pie chart showing the proportion of genes expressed across multiple spatial domains (≥2) versus those specific to a single domain. b Heatmap displaying the spatial enhancer units regulating K1 genes in Fig. 1f. These enhancers were clustered into 17 patterns (M1–M17, columns) based on their spatial specificity across K1 genes expressed domains (D0, D1, D2, D3, D13, and D4). c Schematic representation of divergent regulation and coherent regulation by enhancer clusters. d Bar plot showing the number of enhancer clusters (y-axis) containing enhancers from various numbers of enhancer units (x-axis). Pie chart showing the proportion of divergent regulation (≥2 enhancer units within a cluster) versus coherent regulation (only 1 enhancer pattern). e The combination of enhancers encodes Nrn1 spatial expression. Genome tracks presenting chromatin accessibility, peak site (marked by open regions, enhancer clusters, and enhancer units), gene tracks (left), and gene expression (right) for the Nrn1 locus in cortex regions (D0, D1, D2, and D3). Spatial mapping of denoised chromatin accessibility of representative enhancers from three distinct enhancer units (M12, M3, and M1) and gene expression (top). f Spatial distribution of cell types in spatial-ATAC-RNA-seq of the P22 mouse brain by integration of ATAC pixels and scATAC-seq data. Arrows indicate the dorsal-ventral (D-V) and deep-superficial (D-S) axes in the coronal section. Pixel size, 20 µm; scale bars, 1 mm. g Directionality analysis of chromatin accessibility along D-V or D-S axes for the top 1500 cell-type-specific accessible peaks (ranked by log₂[fold change]) across individual cell types. Each point represents a putative enhancer, positioned by Spearman’s correlation coefficient (R) between peak accessibility and spatial coordinates (n = cell-type-specific counts, coordinates normalized to 0–1 scale), with color indicating statistical significance (−log10 P-value) of correlation. Positive or negative R values denote ventral/superficial or dorsal/deep directional biases, respectively. The top 5 enhancers with the highest absolute R values were explicitly labeled. h Spatial mapping of denoised Hlf gene expression and denoised chromatin accessibility of representative directional enhancers (from g) associated with Hlf showing ITL23GL-specific chromatin accessibility and directionality preference. Schematic representation of the location of Hlf gene and Hlf putative enhancers. i Spatial mapping of denoised Bean1 gene expression and denoised chromatin accessibility of representative directional enhancers (from g) associated with Bean1 showing CTGL-specific chromatin accessibility and directionality preference. Schematic representation of the location of Bean1 gene and Bean1 putative enhancers. Source data are provided as a Source Data file.

Fig. 3. Spatial enhancer code in human tumors.

a Schematic representation (top) and spatial distribution (bottom) of two distinct spatial compartments defined by Slide-tags multiome which was adapted from the original literature. b Spatial distribution of tumor cells and T cells in the human melanoma sample. c Heatmap showing gene expression (left) and enhancer chromatin accessibility (right) for 15,641 gene-enhancer pairs, clustered into 6 modules (K1–K6) based on gene expression in tumor cells and T cells from two compartments. d Bar plot illustrating the number of compartments in which genes are expressed. Pie chart showing the proportion of genes expressed in two compartments versus those specific to one compartment. e Bar plot showing the number of enhancer clusters (y-axis) composed of enhancers displaying various numbers of enhancer units (x-axis). Pie chart depicting the proportion of divergent regulation (≥2 enhancer units within a cluster) versus coherent regulation (only 1 enhancer pattern). f Heatmap displaying the spatial enhancer units regulating module K1 genes in (c). These enhancers were clustered into 3 enhancer units (M1–M3, rows) based on their binarized spatial specificity in two compartments. g Motif enrichment rank of enhancers with spatial patterns (M2 and M3) from (f). Top 10 motifs are highlighted. P values calculated by a one-sided hypergeometric test. h Spatial mapping of denoised DNMT3A gene expression and denoised chromatin accessibility of three DNMT3A enhancers from different enhancer units (M1, M2, and M3). Schematic representation of the location of DNMT3A gene and DNMT3A enhancers. i Spatial distribution of tumor, immune-rich, myeloid, and normal tissue regions in a human breast cancer sample. Schematic representation of spatial-ATAC was from a published study. j A zoom-in window of spatial distribution of tumor, immune-rich, myeloid, and normal tissue regions in the human breast cancer sample from (i). k Left: Heatmap showing gene scores clustered into 9 modules (K1–K9). Right: Heatmap displaying the spatial enhancer units regulating K1 genes in (k). These enhancers were clustered into 3 enhancer units (M1–M3, rows) based on their binarized spatial specificity. l Spatial mapping of denoised NOTCH1 gene expression and denoised chromatin accessibility of three NOTCH1 enhancers from different enhancer units (M1, M2, and M3) in the tumor and myeloid regions. Schematic representation of the location of NOTCH1 gene and NOTCH1 enhancers. Source data are provided as a Source Data file.

Fig. 4. Spatial enhancer code in e11 mouse embryos.

a Spatial distribution of spatial domains in e11 mouse embryos. Pixel size, 20 µm; spatial-ATAC-seq data was from ref. . b Left: Allen RNA ISH images in e11.5 mouse embryos. Right: spatial mapping of gene scores in e11 mouse embryos for selected marker genes across different spatial domains. c Heatmap showing gene expression (left) and chromatin accessibility (right) for 29,485 gene-enhancer pairs. Gene expression was clustered using k-means clustering (k  =  7). The relevant markers and representative enhancers were highlighted. d Spatial mapping of denoised chromatin accessibility (left) and the corresponding VISTA enhancer reporter activity for e11.5 mouse embryos (right) for representative enhancers. e Enrichment of VISTA-validated enhancers among enhancers in (c) and all CREs. *p < 0.05, **p < 0.01, ***p < 0.001; n.s. not significant (binomial test; exact p-value < 2.2e-16). f The percentage of VISTA-validated enhancers in brain, spinal cord, and heart regions overlapped with sCREs from spatial-ATAC data. g Pie chart displaying the proportion of gene expression patterns (top: multiple domains and single domain) and enhancer cluster regulation patterns (bottom: divergent regulation and coherent regulation). h Heatmap showing enhancer units for enhancers regulating genes in K1 as depicted in (c). i The combination of enhancer units encodes Cmtm8 spatial expression. Genome tracks showing the chromatin accessibility, peak site (marked by open region, enhancer cluster, and enhancer unit), gene tracks (left), and gene expression (right) for Cmtm8 in D0 and D1. Spatial mapping of denoised chromatin accessibility of three representative enhancers from different enhancer units (M2, M4, and M6) and gene expression (top). Source data are provided as a Source Data file.

Fig. 5. Spatial enhancer code in e13 mouse brain.

a Spatial distribution of spatial domains in e13 mouse embryos. Pixel size, 50 µm; scale bars, 1 mm. Spatial data was from ref. . b Schematic highlighting different brain regions (forebrain, midbrain, and hindbrain in solid color shades, while the cortex, diencephalon, and prosomere in dotted color lines) of an imaged sagittal slice of an e13.5 mouse brain from the H3K27ac MERFISH data. The background shows the DAPI signal. Scale bars: 1 mm. The picture was adapted from ref. . c Spatial mapping of denoised gene expression (left) and Allen RNA ISH images (right) for Neurog2, Tbr1, and Foxg1. d Spatial mapping of denoised chromatin accessibility (left) and epigenomic MERFISH images of H3K27ac signals (right) for the enhancers controlling Neurog2, Tbr1, and Foxg1. e Heatmap showing gene expression (left) and chromatin accessibility (right) for 14,504 gene-enhancer pairs. Gene expression clustered via k-means clustering (k  =  7). The relevant markers and representative enhancers are highlighted. f Pie chart demonstrating the proportion of gene expression patterns (left: multiple domains and single domain) and enhancer cluster regulation patterns (right: divergent regulation and coherent regulation). g Heatmap showing the enhancer units for enhancers regulating K1 genes as depicted in (e). h The combination of enhancer units encodes Bcl11a spatial expression. Genome tracks showing the chromatin accessibility, peak site (labeled by open region, enhancer cluster, and enhancer unit), gene tracks (left), and gene expression (right) for Bcl11a in D0 and D2. Spatial mapping of denoised chromatin accessibility of three representative enhancers in distinct enhancer units (M2, M4, and M6) and gene expression (top). Source data are provided as a Source Data file.

Fig. 6. Experimental validation of Atoh1 spatial enhancer code in e11 mouse embryos.

Spatial mapping of Atoh1 gene activity scores (a) in the spinal cord and hindbrain in e11 mouse embryos (b), where spatial-ATAC-seq data were from a published study. c Fluorescence in situ hybridization (FISH) images depicting Atoh1 expression at e11 mouse embryos. The experiment was repeated independently 3 times with similar results. Scale bars, 500 µm. d Heatmap showing the accessibility levels of Atoh1 enhancers. The spatial specificity scores of Atoh1 enhancers in the e11 spinal cord and hindbrain were labeled on the heatmap. e Genome tracks showing the chromatin accessibility, peak coordinates, gene activity sores (left), and gene tracks (right) for the Atoh1 locus in the spinal cord and hindbrain of e11 mouse. Spatial mapping of denoised chromatin accessibility of four Atoh1 enhancers (E0, E1, E2, and E3) in the e11 mouse spinal cord and hindbrain. f Whole-mount X-gal staining of E0 (top), E2 (middle), and E3 (bottom) reporter embryos at e11. The blue signal displays the β-gal activity of the LacZ reporter driven by the indicated enhancers, representing the activities of the Atoh1 enhancers. The experiment was repeated independently 3 times with similar results. Scale bars, 500 µm. g, h Schematic diagram represented sampled embryonic tissue (left), imaging perspectives (middle), and eSpatial prediction (right). The blue dot line indicates the section direction. The pink box denotes the RNAscope imaged region in the spinal cord (g). The pale green box denotes the RNAscope imaged region of the hindbrain (h). i–l RNAscope detection of Atoh1 mRNA in E0/E2/E3/E0 + 3 KO or WT mice in e11 spinal cord. Left: Representative images show staining for Atoh1 RNAscope probes (red) and DAPI (blue). Right: Quantification of Atoh1 RNAscope probe signal. Values shown are the mean ± standard error of the mean; n = 5 independent experiments. Statistical significance was determined by two-sided Student’s t-test: E0-WT vs. E0-KO (p < 0.0001), E2-WT vs. E2-KO (p = 0.8792), E3-WT vs. E3-KO (p = 0.0296), and E0 + 3-WT vs. E0 + 3-DKO (p < 0.0001). *p < 0.05, **p < 0.01, ***p < 0.001; n.s., not significant. Scale bar, 50 μm. WT wild type, KO knockout. m–p RNAscope detection of Atoh1 mRNA in E0/E2/E3/E0 + 3 KO or WT mice in e11 hindbrain. Left: Representative images show staining for Atoh1 RNAscope probes (red) and DAPI (blue). Right: Quantification of Atoh1 RNAscope probe signal. Values shown are the mean ± standard error of the mean; n = 5 independent experiments. Statistical significance was determined by two-sided Student’s t-test: E0-WT vs. E0-KO (p < 0.0001), E2-WT vs. E2-KO (p = 0.9910), E3-WT vs. E3-KO (p = 0.8914), and E0 + 3-WT vs. E0 + 3-DKO (p < 0.0001). *p < 0.05, **p < 0.01, ***p < 0.001; n.s., not significant. Scale bar, 50 μm. WT wild type, KO knockout. Source data are provided as a Source Data file.