Spatial joint profiling of DNA methylome and transcriptome in tissues

Tissue slide preparation

Mouse C57 embryo sagittal frozen sections (MF-104-11-C57 and MF-104-13-C57) were purchased from Zyagen. Freshly collected E11 or E13 mouse embryos were snap frozen in optimal cutting temperature (O.C.T.) compounds and sectioned at 7–10 μm thickness. Tissue sections were collected on poly-l-lysine-coated glass slides (Electron Microscopy Sciences, 63478-AS).

Juvenile mouse brain tissue (P21) was obtained from the C57BL/6 mice housed in the University of Pennsylvania Animal Care Facilities under pathogen-free conditions. All procedures used were approved by the Institutional Animal Care and Use Committee.

Mice were euthanized at P21 using CO₂ inhalation, followed by transcranial perfusion with cold Dulbecco’s PBS (DPBS). After isolation, brains were embedded in Tissue-Tek O.C.T. compound and snap frozen on dry ice and a 2-methylbutane bath. Coronal cryosections of 8–10 μm were mounted on the back of Superfrost Plus microscope slides (Fisher Scientific, 12-550-15).

Preparation of transposome

Unloaded Tn5 transposome (C01070010) was purchased from Diagenode and the transposome was assembled following the manufacturer’s guidelines. The oligonucleotides used for transposome assembly were: Tn5ME-B, 5′-/5Phos/CATCGGCGTACGACTAGATGTGTATAAGAGACAG-3′; Tn5MErev, 5′-/5Phos/CTGTCTCTTATACACATCT-3′.

DNA barcode sequences, DNA oligonucleotides and other key reagents

DNA oligonucleotides used for PCR and library construction are shown in Supplementary Table 1. All DNA barcode sequences are provided in Supplementary Tables 2 (barcode A) and 3 (barcode B) and all other chemicals and reagents are listed in Supplementary Table 4.

Fabrication of the polydimethylsiloxane 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 (Microchem, SU-2025 and SU-2010) onto a silicon wafer (WaferPro, C04004). 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 onto the moulds. After degassing for 30 min, the mixture was cured at 66–70 °C for 2–16 h. Solidified PDMS was extracted from the moulds for further use. The detailed protocol for the fabrication and preparation of the PDMS device can be found in our previous research²⁴.

Spatial joint profiling of DNA methylation and RNA transcription

Frozen tissue slides were quickly thawed for 1 min in a 37 °C incubator. The tissue was fixed with 1% formaldehyde in PBS containing 0.05 U ml⁻¹ RNase inhibitor (Enzymatics) for 10 min and quenched with 1.25 M glycine for another 5 min at room temperature. After fixation, tissue was washed twice with 1 ml of DPBS–RNase inhibitor and cleaned with deionized H₂O.

The tissue was subsequently permeabilized with 100 μl of 0.5% Triton X-100 plus 0.05 U ml⁻¹ RNase inhibitor for 30 min at room temperature, then washed twice with 200 μl DPBS–RNase inhibitor for 5 min each. After permeabilization, the tissue was treated with 100 μl of 0.1 N HCl for 5 min at room temperature to disrupt histones from the chromatin, then washed twice with 200 μl of wash buffer (10 mM Tris-HCl pH 7.4, 10 mM NaCl, 3 mM MgCl₂, 1% BSA and 0.1% Tween 20) plus 0.05 U ml⁻¹ RNase inhibitor for 5 min at room temperature. Next, 50 μl of transposition mixture (5 μl of assembled transposome, 16.5 μl of 1× DPBS, 25 μl of 2× Tagmentation buffer, 0.5 μl of 1% digitonin, 0.5 μl of 10% Tween 20, 0.05 U ml⁻¹ RNase inhibitor (Enzymatics) and 1.87 μl nuclease-free water) was added and incubated at 37 °C for 60 min. After 60 min incubation, the first round of transposition mixture was removed and a second round of 50 μl of fresh transposition mixture was added and incubated for another 60 min at 37 °C. To stop the transposition, 200 μl of 40 mM EDTA with 0.05 U ml⁻¹ RNase inhibitor was added with incubation for 5 min at room temperature. After that, 200 μl 1× NEB3.1 buffer plus 1% RNase inhibitor was used to wash the tissue for 5 min at room temperature. The tissue was then washed again with 200 μl of DPBS–RNase inhibitor for 5 min at room temperature before proceeding with the in situ reverse transcription reaction.

In situ reverse transcription

For the in situ reverse transcription, the following mixture was added: 12.5 μl 5× reverse transcription buffer, 4.05 μl RNase-free water, 0.4 μl RNase inhibitor (Enzymatics), 1.25 μl 50% PEG-8000, 3.1 μl 10 mM dNTPs, 6.2 μl 200 U μl⁻¹ Maxima H Minus Reverse Transcriptase, 25 μl 0.5× DPBS–RNase inhibitor and 10 μl 100 μM reverse transcription primer (biotinylated-dT oligo). The tissue was incubated for 30 min at room temperature, then at 45 °C for 90 min in a humidified container. After the reverse transcription reaction, tissue was washed with 1× NEB3.1 buffer plus 1% RNase inhibitor for 5 min at room temperature.

In situ barcoding

For in situ ligation with the first barcode (barcode A), the first PDMS chip was covered at the tissue ROI. For alignment purposes, a 10× objective (KEYENCE BZ-X800 fluorescence microscope, BZ-X800 Viewer Software) was used to take the bright-field image. The PDMS device and tissue slide were clamped tightly with a custom acrylic clamp. Barcode A was first annealed with ligation linker 1 by mixing 10 μl of 100 μM ligation linker, 10 μl of 100 μM individual barcode A and 20 μl of 2× annealing buffer (20 mM Tris-HCl pH 7.5–8.0, 100 mM NaCl₂ and 2 mM EDTA). For each channel, 5 μl of ligation master mixture was prepared with 4 μl of ligation mixture (27 μl T4 DNA ligase buffer, 0.9 μl RNase inhibitor (Enzymatics), 5.4 μl 5% Triton X-100, 11 μl T4 DNA ligase and 71.43 μl RNase-free water) and 1 μl of each annealed DNA barcode A (A1–A50, 25 μM). Vacuum was applied to flow the ligation master mixture into the 50 channels of the device and cover the ROI of the tissue, followed by incubation at 37 °C for 30 min in a humidified container. Then the PDMS chip and clamp were removed after washing the tissue with 1× NEB 3.1 buffer for 5 min. The slide was then washed with deionized water and dried using compressed air.

For in situ ligation with the second barcode (barcode B), the second PDMS chip was covered at the ROI and a bright-field image was taken with the 10× objective. An acrylic clamp was applied to clamp the PDMS and tissue slide together. Annealing of barcode B (B1–B50, 25 μM) and preparation of the ligation mixture are the same as barcode A. The whole device was incubated at 37 °C for 30 min in a humidified container. The PDMS chip and clamp were then removed, and the slide was washed with deionized water and dried using compressed air. A bright-field image was then taken for further alignment.

Reverse crosslinking

For tissue lysis, the ROI was digested with 100 μl of the reverse crosslinking mixture (0.4 mg ml⁻¹ proteinase K, 1 mM EDTA, 50 mM Tris-HCl pH 8.0, 200 mM NaCl and 1% SDS) at 58–60 °C for 2 h in a humidified container. The lysate was then collected in a 0.2-ml PCR tube and incubated on a 60 °C shaker overnight.

gDNA and cDNA separation

For gDNA and cDNA separation, the lysate was purified with Zymo DNA Clean and Concentrator kit and eluted with 100 μl nuclease-free water. Next, 40 μl of Dynabeads MyOne Streptavidin C1 beads were used and washed three times with 1× B&W buffer containing 0.05% Tween 20 (50 μl 1 M Tris-HCl pH 8.0, 2,000 μl 5 M NaCl, 10 μl 0.5 M EDTA, 50 μl 10% Tween 20 and 7,890 μl nuclease-free water). After washing, beads were resuspended in 100 μl of 2× B&W buffer (50 μl 1 M Tris-HCl pH 8.0, 2,000 μl 5 M NaCl, 10 μl 0.5 M EDTA and 2,940 μl nuclease-free water) containing 2 μl of SUPERase In RNase inhibitor, then mixed with the gDNA–cDNA lysate and allowed to bind for 1 h with agitation at room temperature. A magnet was then used to separate the beads, which bind to the cDNA that contains dT, from the supernatant that contains the gDNA.

gDNA library generation

Supernatant (200 μl) was collected from the above separation process for further methylated gDNA detection and library construction. Next, 1 ml of DNA binding buffer was added to the 200 μl supernatant and purified with the Zymo DNA Clean and Concentrator kit again, then eluted in 84 μl (3 × 28 μl) nuclease-free water. The NEBNext enzymatic methyl-seq conversion module (EM-seq) was then used to detect methylated DNA in the sample by converting unmethylated cytosines to uracil; the manufacturer’s guidelines were followed. Then, 28 μl of DNA sample was aliquoted into each PCR tube, TET2 reaction mixture (10 μl TET2 reaction buffer containing reconstituted TET2 reaction buffer supplement, 1 μl oxidation supplement, 1 μl DTT, 1 μl oxidation enhancer and 4 μl TET2) was added to the DNA sample on ice. In brief, 5 μl of diluted 1:1,300 of 500 mM Fe (II) solution was added to the mixture and incubated for 1 h at 37 °C in a thermocycler. After the reaction, the sample was transferred to ice and 1 μl of stop reagent from the kit was added. The sample was then incubated for another 30 min at 37 °C. TET2 converted DNA was then purified with 90 μl of solid-phase reversible immobilization (SPRI) beads and eluted with 16 μl nuclease-free water. The thermocycler was preheated to 85 °C, 4 μl formamide was added to the converted DNA and incubated for 10 min at 85 °C in the preheated thermocycler. After the reaction, the heated sample was immediately placed on ice to maintain the open chromatin structure, then reagents from the kit were added (68 μl nuclease-free water, 10 μl APOBEC reaction buffer, 1 μl BSA and 1 μl APOBEC) to deaminate unmethylated cytosines to uracil for 3 h at 37 °C in a thermocycler. Deaminated DNA was then cleaned up using 100 μl (1:1 ratio) of SPRI beads and eluted in 20 μl nuclease-free water.

Splint ligation

The gDNA tube was heat-shocked for 3 min at 95 °C and immediately put on ice for 2 min. Then, 10 μl of 0.75 μM pre-annealed Splint Ligate P5 (SLP5) adapter was added. This adapter was diluted from a 12 μM stock, which contained 6 μl of 100 μΜ SLP5RC oligo, 8.4 μl of 100 μΜ SLS5ME-A-H10 oligo, 5 μl of 10× T4 RNA Ligase Buffer and 30.6 μl nuclease-free water in a PCR tube that was incubated at 95 °C for 1 min, then gradually cooled by −0.1 °C s⁻¹ to 10 °C on a thermocycler. Next, 80 μl of ligation master mixture was added to the gDNA tube at room temperature. The mixture contained 40 μl preheated 50% PEG-8000, 12.5 μl SCR buffer (666 mM Tris-HCl pH 8.0 and 132 mM MgCl₂ in nuclease-free water), 10 μl of 100 mM DTT, 10 μl of 10 mM ATP, 1.25 μl of 10,000 U ml⁻¹ T4 PNK and 6.25 μl of 400,000 U ml⁻¹ T4 ligase. The splint ligation mixture was then splinted into five 0.2-ml PCR tubes, 20 μl per tube. The tubes were shaken at 1,000 rpm for 10 s and spun down, then incubated for 45 min at 37 °C, followed by 20 min at 65 °C to inactivate the ligase. For splint ligation indexing PCR, 80 μl of the PCR reaction mixture was mixed in each splint ligated tube. The mixture contained 20 μl 5× VeraSeq GC Buffer, 4 μl 10 mM dNTPs, 3 μl VeraSeq Ultra Enzyme, 5 μl 20× EvaGreen dye, 2 μl of 10 μM N501 primer and 2 μl of 10 μM N70X-HT primer (Supplementary Table 1). The mixture was then aliquoted into a clean PCR tube with 50 μl volume and run on a thermocycler with the setting below, 98 °C for 1 min, then cycling at 98 °C for 10 s, 57 °C for 20 s and 72 °C for 30 s, for 13–19 cycles, followed by 72 °C for 10 s. The reaction was removed once the quantitative PCR (qPCR) signal began to plateau. The amplified PCR products were pooled and purified with a 0.8× volume ratio of SPRI beads (bead-to-sample ratio) and the completed DNA library was eluted in 15 μl nuclease-free water.

cDNA library generation

The separated beads containing cDNA were used for cDNA library generation. In brief, 400 μl of 1× B&W buffer with 0.05% Tween 20 was used to wash the beads twice. Then, the beads were washed once with 400 μl of 10 mM Tris-HCl pH 8.0 containing 0.1% Tween 20 for 5 min at room temperature. Streptavidin beads with bound cDNA molecules were placed on a magnetic rack and washed once with 250 μl nuclease-free water before being resuspended in a template switching oligonucleotide solution (44 μl 5× Maxima reverse transcription buffer, 44 μl of 20% Ficoll PM-400 solution, 22 μl of dNTPs, 5.5 μl of 100 mM template switch oligo, 11 μl Maxima H Minus reverse transcriptase, 5.5 μl of RNase inhibitor (Enzymatics) and 88 μl nuclease-free water). Resuspended beads were then incubated for 30 min with agitation at room temperature and for 90 min at 42 °C, with gentle agitation. After the reaction, beads were washed with 400 μl of 10 mM Tris pH 8.0 containing 0.1% Tween 20 and washed without resuspension in 250 μl nuclease-free water. Water was removed on the magnetic rack and the beads were resuspended in the PCR solution (100 μl of 2× Kappa Master mix, 8.8 μl of 10 μM primers 1 and 2 and 92.4 μl nuclease-free water). Next, the beads were mixed well and 50 μl of the PCR mixture was split into four 0.2-ml PCR tubes. The PCR programme was run as follows: 95 °C for 3 min and cycling at 98 °C for 20 s, 65 °C for 45 s and 72 °C for 3 min, for a total of five cycles, followed by 4 °C on hold. After five cycles of PCR reaction, four PCR tubes were placed on a magnetic rack and 50 μl of the clear PCR solution was transferred to four optical-grade qPCR tubes, adding 2.5 μl of 20× Evagreen dye to each tube. The sample was run on a qPCR machine with the following conditions: 95 °C for 3 min, cycling at 98 °C for 20 s, 65 °C for 20 s and 72 °C for 3 min, for 13–17 cycles, followed by 72 °C for 5 min. The reaction was removed once the qPCR signal began to plateau. The amplified PCR product was purified with a 0.8× volume ratio of SPRI beads and eluted in 20 μl nuclease-free water.

A Nextera XT DNA Library Prep Kit was used for cDNA library preparation. In brief, 2 μl (2 ng) of purified cDNA (1 ng μl⁻¹), 10 μl Tagment DNA buffer, 5 μl Amplicon Tagment mix and 3 μl nuclease-free water were mixed and incubated at 55 °C for 5 min. Then, 5 μl NT buffer was added to stop the reaction with incubation at room temperature for 5 min. PCR master mix (15 μl 2× N.P.M. Master mix, 1 μl of 10 μM P5 primer (N501) and 1 μl of 10 μM indexed P7 primer (N70X) and 8 μl nuclease-free water) was added. The PCR reaction was run 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 a total of 12 cycles. The PCR product was then purified with a 0.7× ratio of SPRI beads and eluted in 15 μl nuclease-free water to obtain the cDNA library.

Library quality check and next-generation sequencing

An Agilent Bioanalyzer D5000 ScreenTape was used to determine the size distribution and concentration of the library before sequencing. Next-generation sequencing was conducted on an Illumina NovaSeq 6000 sequencer and NovaSeq X Plus system (150 bp paired-end mode).

Data preprocessing

For RNA-sequencing data, Read 2 was processed to extract barcode A, barcode B and the UMIs. Using the STARsolo pipeline⁵⁶ (v.2.7.10b), these processed data were mapped to the mouse genome reference (mm10). This step generated a gene matrix that captures both gene-expression and spatial-positioning information, encoded through the combination of barcodes A and B. The gene matrix was then imported into R for downstream spatial transcriptomic analysis using Seurat package (v.5.1.0)⁵⁷.

For DNA-methylation data, adaptor sequences were trimmed before demultiplexing the FASTQ files using the combination of barcodes A and B. We used the BISulfite-seq CUI Toolkit (BISCUIT) (v.0.3.14)⁵⁸ to align the DNA sequences to the mouse reference genome (mm10). Methylation levels at individual CG and CH sites were stored as continuous values between 0 and 1, representing the fraction of methylated reads after quality filtering. These processed CG–CH files were then analysed independently using the MethSCAn pipeline¹⁸ to identify VMRs⁵⁹, defined as fused genome intervals with methylation-level variance in the top 2%. We used default parameter settings when running MethSCAn, and the MethSCAn filter min-sites parameter was determined from the read coverage knee plot (Extended Data Fig. 3a). The methylation levels and residuals of VMRs were then imported into R for downstream DNA-methylation analysis.

Clustering and data visualization

We mapped the exact location of pixels on the bright-field tissue image using a custom Python script (https://github.com/zhou-lab/Spatial-DMT-2024/tree/main/Data_preprocess/Image), before removing additional empty barcodes on the basis of read-count thresholds determined by the knee plot (Extended Data Fig. 3a). Clustering and data visualization were conducted using R in RStudio.

For RNA data, we used the SCTtransform function in the Seurat package (v.5.1.0), built using a regularized negative binomial model, for normalization and variance stabilization. Dimensionality reduction was performed using RunPCA function with the SCTtransformed assay. We then constructed the nearest-neighbour graph using the first 30 principal components with the FindNeighbors function and identified clusters with the default Leiden method in the FindClusters function. Finally, a UMAP embedding was computed using the same principal components with RunUMAP function.

Owing to the inherent sparsity of DNA-methylation data, it is impractical to analyse methylation status solely at the individual CpG level. Binary information at sparse loci cannot be used directly to construct a feature matrix suitable for downstream analysis. In our study, we adopted the VMR framework, which divides the genome into variable-sized tiles and calculates the average methylation level across CpGs in each tile for each pixel¹⁸. This approach results in a continuous-valued matrix, in which rows correspond to pixels and columns represent genomic tiles, with values ranging from 0 to 1. VMR methylation levels and residuals were then imputed using the iterative principal component analysis approach as suggested in the MethSCAn instructions. Initially, missing residual values were replaced with zero and missing methylation levels were replaced with the average values for that VMR interval. The principal component analysis approach was iteratively applied until updated values stabilized to a threshold. The imputed residual matrix for VMRs was then imported into the existing Seurat object as another modality. Similar to the RNA-clustering pipeline, dimensionality was reduced using the RunPCA function. The first ten principal components from the residual matrix were used for clustering and UMAP embedding.

To visualize clusters in their spatial locations, the SpatialDimPlot function was used after clustering on the basis of gene expression or VMR residuals. UMAP embedding was visualized with the DimPlot function. The FindMarkers function was applied to select genes and VMRs that were differentially expressed or methylated for each cluster. For spatial mapping of individual VMR methylation levels or gene expression, we applied the smoothScoresNN function from the FigR package⁶⁰. The SpatialFeaturePlot function was then used to visualize VMR methylation levels and gene expression across all pixels. To illustrate the relationships between clustering results from different modalities, we generated the confusion matrix and alluvial diagram using the pheatmap and ggalluvial R package⁶¹.

Integrative analysis of DNA methylation and RNA data

To integrate spatial DNA methylation and RNA data, WNN analysis in Seurat was applied using the FindMultiModalNeighbors function¹⁹. On the basis of the WNN graph, clustering, UMAP embedding and spatial mapping of identified clusters were performed for integrated visualization.

For the integration of spatial transcriptomics data of E11 and E13 mouse embryos, the top 3,000 integration features were selected, followed by the use of PrepSCTIntegration and IntegrateData functions to generate an integrated dataset. Similarly, to integrate with public single-cell transcriptomic data^25,44, we first identified anchors using the FindIntegrationAnchors function in Seurat, followed by data integration using the IntegrateData function. To integrate DNA-methylation data, common VMRs between both developmental stages were obtained and the integrated CCA method from the IntegrateLayers function was used to join the methylation data from the two developmental stages. A Wilcoxon signed-rank test was performed to compare the methylation levels and gene-expression differences between the two time points.

TF motif enrichment

To perform TF motif enrichment, we first used the MethSCAn diff function on distinct groups of cells to identify differentially methylated VMRs on the basis of the clustering assignment. The HOMER⁶² findMotifsGenome function was then applied to analyse the enrichment of known TF motifs using its default database. We followed the same parameter settings used in MethSCAn, with motif lengths of 5, 6, 7, 8, 9, 10, 11 and 12.

CpGs enrichment analysis

Enrichment analysis of individual CpGs in the differential regions (Fig. 5a,b) was performed using knowYourCG (https://github.com/zhou-lab/knowYourCG), which provides a comprehensive annotation database for each CpG, including chromatin states, TF binding sites, motif occurrences, PMD annotations and more. To avoid inflated odds ratios for high-coverage data, genomic uniformity was quantified using fold enrichment, defined as the ratio of observed overlaps to expected overlaps. The expected number of overlaps was calculated as: (number of CpGs sequenced × number of CpGs in the chromatin state feature)/total number of CpGs in the genome.

Correlation and GO enrichment analysis

Correlation analysis was performed for different clusters. We first used the findOverlaps function in GenomicRanges package (v.4.4)⁶³ to map VMRs to overlapped genes. Then, the Pearson correlation test was applied to obtain the correlation between mapped genes and corresponding VMRs. The Benjamini–Hochberg procedure was used to adjust all P values.

GO enrichment analysis was conducted using the enrichGO function from clusterProfiler package⁶⁴ (v.4.2). For GO enrichment in the comparative analysis of E11 and E13 mouse embryos, the FindMarkers function in Seurat package was used to find differential genes and VMRs in the same cluster from integrated data across two developmental stages. Differentially upregulated genes (false discovery rate ≤ 0.05) with demethylated VMRs were used for the GO analysis.

Reporting summary

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