Temporally distinct 3D multi-omic dynamics in the developing human brain

Ethics statement

Paediatric tissues obtained from the University of California, San Francisco (UCSF) were collected from autopsy sources through the UCSF Pediatric Neuropathology Research Laboratory. Informed consent was obtained from the next of kin for all paediatric samples obtained from the Pediatric Neuropathology Research Laboratory. Paediatric samples collected through autopsy were de-identified before acquisition and thus exempt from Institutional Review Board (IRB) review. For tissues obtained through the gynaecology clinic, patients were asked about their interest in donating tissue to research after making the decision for termination of pregnancy. Patients who agreed signed a written consent after receiving information, both written and oral, given by a physician or midwife. They were informed that agreeing to donate would not affect their medical care and that neither the donor nor the clinical team would benefit from the donation. The use of abortion material was reviewed and approved by the UCSF Committee on Human Research. Protocols were approved by the Human Gamete, Embryo, and Stem Cell Research Committee (IRB GESCR number 10-02693; IRB number 20-31968) at UCSF.

Adult human brain samples and a post-mortem GW-35 sample (based on adjusted age) were banked by the National Institutes of Health (NIH) NeuroBioBank at the University of Maryland Brain and Tissue Bank (Supplementary Table 1). Informed consent was obtained from the patient or next of kin for all samples obtained from NIH NeuroBioBank. The tissue collection and repository is overseen by The University of Maryland IRB with IRB protocol number HM-HP-00042077, as well as The Maryland Department of Health and Mental Hygiene IRB with IRB protocol number 5-58. When an individual of any age dies, the medical examiner or coroner contacts the next of kin and asks whether they would be willing to talk to a staff member at the University of Maryland Brain and Tissue Bank about an NIH-funded tissue procurement project. If the family agrees, the medical examiner or coroner contacts the bank, and a staff member obtains a recorded telephone consent from the next of kin for donation. A written verification of the consent is then faxed to the office of the referring medical examiner or coroner. Alternatively, the recording can also be played over the telephone for confirmation. The University of California, Los Angeles IRB has determined that our study using post-mortem human tissue obtained from the NIH NeuroBioBank involves no human participant and requires no IRB review.

Brain specimens

Collections were carried out at post-mortem intervals of less than 24 h. Tissue was collected through the UCSF, and the NIH NeuroBioBank. Patient consent was obtained before collection according to institutional ethical regulations. The UCSF Committee on Human Research reviewed protocols that were approved by the Human Gamete, Embryo and Stem Cell Research Committee (IRB GESCR number 10-02693) and associated IRBs (20-31968) at UCSF. For prenatal samples, mothers gave consent for the related medical procedure before any consent for donation of tissue. In the tissue consent process, they were informed that agreeing to donate would not affect their medical care and that neither the donor nor the clinical team would benefit from the donation. Specimens collected from UCSF were evaluated by a neuropathologist as control samples collected through the UCSF Pediatric Neuropathology Research Laboratory. Additional samples were identified from the NIH NeuroBioBank, according to their IRB approval. Tissues were cut coronally, and areas of interest were sampled. Tissue blocks of 1 mm used for the snm3C-seq3 assay were flash-frozen in liquid nitrogen and stored at −80 °C. Blocks used for histological analyses were fixed with 4% paraformaldehyde (PFA) for 2 days and cryoprotected in a 30% sucrose gradient. The tissue was then frozen in OCT, and blocks were cut at 30 µm with a cryostat and mounted onto glass slides. For each sample used, we cresyl-stained three sections spanning the block to ensure our position using anatomical landmarks, such as the lateral ventricle, presence of the caudate, thalamus and HPC.

snm3C-seq3

For prenatal brain samples, snm3C-seq3 was carried out without the labelling of neuronal nuclei using anti-NeuN antibody, whereas postnatal samples were labelled by anti-NeuN antibody during the procedure to isolate nuclei. For snm3C-seq3 carried out without labelling, frozen powder of brain tissue was resuspended in 10 ml of DPBS with 2% formaldehyde and incubated at room temperature for 10 min with slow rotation. The crosslinking reaction was quenched with 1.17 ml of 2 M glycine for 5 min at room temperature. The crosslinked tissue sample was pelleted by centrifugation at 2,000g for 10 min at 4 °C. The same centrifugation condition was used to pellet nuclei throughout the snm3C-seq3 procedure. The pellet was resuspended in 3 ml NIBT (10 mM Tris-HCl pH 8.0, 0.25 M sucrose, 5 mM MgCl₂, 25 mM KCl, 1 mM dithiothreitol, 0.1% Triton X-100 and 1:100 protease inhibitor cocktail (Sigma number P8340)). The resuspended tissue sample was dounced with a dounce homogenizer (Sigma number D9063) 40 times with a loose pestle and 40 times with a tight pestle. For snm3C-seq3 carried out with anti-NeuN labelling, anti-NeuN antibody (PE-conjugated, clone A60, Millipore-Sigma number FCMAB317PE) was added to NIBT at a 1:250 dilution and was incubated with the tissue lysate during the homogenization steps for a total of 15 min. The lysate was mixed with 2 ml of 50% iodixanol (prepared by mixing OptiPrep density gradient medium (Sigma number D1556) with diluent (120 mM Tris-Cl pH 8.0, 150 mM KCl and 30 mM MgCl₂) with a volume ratio of 5:1). The lysate was gently layered on top of a 25% iodixanol cushion and centrifuged at 10,000g for 20 min at 4 °C using a swing rotor. The pellet of nuclei was resuspended in 1 ml of cold DPBS followed by quantification of nuclei using a Biorad TC20 Automated Cell Counter (Biorad number 1450102).

The in situ 3C reaction was carried out using an Arima Genomics Arima-HiC+ kit. Each in situ 3C reaction used 300,000 to 450,000 nuclei. Nuclei aliquots were pelleted and resuspended in 20 µl H₂O mixed with 24 µl conditioning solution and incubated at 62 °C for 10 min. After the incubation, 20 µl of stop solution 2 was added to the reaction and incubated at 37 °C for 15 min. A restriction digestion mix containing 7 µl of 10× NEB CutSmart buffer (NEB number B7204), 4.5 µl of NlaIII (NEB number R0125), 4.5 µl of MboI (NEB number R0147) and 12 µl of 1× NEB CutSmart buffer was added to the reaction followed by incubation at 37 °C for 1 h. The restriction digestion reaction was stopped by incubation at 65 °C for 20 min. A ligation mix containing 70 µl of buffer C and 12 µl of enzyme C was added and then incubated at room temperature for 15 min. The reaction was then kept at 4 °C overnight.

Before fluorescence-activated nucleus sorting, 900 µl cold DPBS supplemented with 100 µl ultrapure BSA (50 mg ml^−1, Invitrogen number AM2618) was added to the in situ 3C reaction. To fluorescently stain nuclei, 1 µl of 1 mg ml^−1 Hoechst 33342 was added before sorting. Fluorescence-activated nucleus sorting was carried out at the Broad Stem Cell Research Center Flow Cytometry core of the University of California, Los Angeles using BD FACSAria sorters. Single nuclei were sorted into 384-well plates containing 1 µl M-Digestion buffer containing proteinase K and about 0.05 pg lambda DNA isolated from dcm⁺ Escherichia coli (Promega number D1501).

Single-nucleus DNA methylome library preparation with snmC-seq3

snmC-seq3 is a modification of snmC-seq2⁴³ that provides improved throughput and reduced cost. Key differences between snmC-seq3 and snmC-seq2 include the usage of 384 instead of 8 barcoded degenerated (RP-H) primers (Supplementary Table 8) for the initiation of random-primed DNA synthesis using bisulfite-converted DNA as a template. The expanded multiplexing allows the combination of 64 single nuclei into the downstream enzymatic reactions, which provides an eightfold reduction of the usage of Adaptase and PCR reagents. In addition, the amounts of Klenow exo^−, exonuclease 1 and rSAP are reduced by tenfold compared to snmC-seq2, further reducing reagent cost. A detailed bench protocol for snm3C-seq3 is provided through protocol.io (https://doi.org/10.17504/protocols.io.kqdg3x6ezg25/v1).

Probe library design for RNA MERFISH and chromatin tracing

We selected 40-bp target sequences for DNA or RNA hybridization by considering each contiguous 40-bp subsequence of each target of interest (the mRNA of a targeted gene or the genomic locus of interest) and then filtering out off-targets to the rest of the transcriptome or genome including repetitive regions, or too high or low GC content or melting temperature. More specifically, our probe design algorithm was implemented with three steps: build a 17-base index based on reference genome hs1 assembly (DNA) or the hg38 transcriptome (RNA); quantify 17-base off-target counts for each candidate 40-bp target sequence; filter and rank target sequences on the basis of predefined selection criteria as previously described^30,44.

MERFISH gene selection

The MERFISH gene selection was carried out by first using a BICCN dataset from GW-18–19 brain and using NSForest v2⁴⁵ with default parameters to identify marker genes for the cell-type clusters in this data⁴⁶. This list of genes was supplemented with additional marker genes from the literature (that is, DCX, GFAP and so on) as well as genes with differential methylation in the snm3C-seq3 data in HPC of mid-gestational human brains. The target sequences for each gene were concatenated with one or two unique readout sequences to facilitate MERFISH or single-molecule fluorescence in situ hybridization imaging. The final list of encoding probes for RNA imaging used is shown in Supplementary Table 5.

Design of chromatin probes

We designed probes for DNA hybridization similarly to those for RNA MERFISH as described in refs. ^30,47. Briefly, we first partitioned chromosome 14 into 50-kb segments and selected for imaging a fifth of these segments uniformly spaced every 250 kb (amounting to 354 target genomic loci using genome reference hs1). After screening against off-target binding, GC content and melting temperature, about 150 unique 40-bp target sequences were selected for each 50-kb segment. We concatenated a unique readout sequence to the target probes of each segment to facilitate sequential hybridization and imaging of each locus. The final list of encoding probes for DNA imaging used is shown in Supplementary Table 4.

Primary probe synthesis

The encoding probes were synthesized from template oligonucleotide pools, following the previously described method⁴⁸. First, we amplified the oligonucleotide pools (Twist Biosciences) using a limited cycle quantitative PCR (approximately 15–20 cycles) with a concentration of 0.6 µM of each primer to create templates. These templates were converted into the corresponding RNAs using the in vitro T7 transcription reaction (New England Biolabs, E2040S) and the PCR product as the templates. The resulting RNAs were then converted to complementary single-stranded DNA using reverse transcription. During the reverse transcription step, a primer with a 5′ acrydite-modified end was used to facilitate the incorporation of the encoding probes into a protective acrylamide gel cast on the sample allowing for more than 100 rounds of hybridization without substantial probe loss. Subsequently, the single-stranded DNA oligonucleotides were purified using alkaline hydrolysis to remove RNA templates and cleaned with columns (Zymo, D4060). The resulting probes were stored at −20 °C.

Sample preparation for multiplexed imaging experiments

The samples were prepared similarly to previously described for cell culture samples with notable modifications³⁰. Briefly, fresh frozen brain tissues were sectioned into coronal sections of 18 μm thickness at −20 °C using a Leica CM3050S cryostat. Sections were collected on salinized and poly-l-lysine (Millipore, 2913997)-treated 40-mm, round number 1.5 coverslips (Bioptechs, 0420-0323-2). The tissue sections were fixed with 4% PFA (Electron Microscopy Sciences, 15710) in 1× PBS with RNase inhibitors (New England Biolabs, M0314L) for 10 min at room temperature before being permeabilized with 0.5% Triton X-100 (Sigma-Aldrich, T8787). Then coverslips were treated with 0.1 M hydrochloric acid (Thermo Scientific, 24308) for 5 min at room temperature. Tissue sections were next incubated with pre-hybridization buffer (40% (vol/vol) formamide (Ambion, AM9342) in 2× SSC (Corning, 46-020-CM)) for 10 min. Then 50 μl of encoding probe hybridization buffer (50% (vol/vol) formamide (Ambion, AM9342), 2× SSC, 10% dextran sulfate (Millipore, S4030)) containing 15 μg of RNA-encoding probes for the targeted genes and 150 μg of DNA-encoding probes for the targeted chromosome 14 loci was incubated with the sample first at 90 °C for 3 min, followed by 47 °C for 18 h. The sample was then washed with 40% (vol/vol) formamide in 2× SSC, 0.5% Tween 20 for 30 min before being embedded in thin, 4% polyacrylamide gels as described previously⁴⁹.

Imaging and adaptor hybridization protocol

MERFISH measurements were conducted on a custom microfluidics–microscope system with the configuration described previously described^30,44. Briefly, the system was built around a Nikon Ti-U microscope body with a Nikon CFI Plan Apo Lambda 60× oil immersion objective with 1.4 NA and used a Lumencor CELESTA light engine and a scientific CMOS camera (Hamamatsu FLASH4.0). The different components were synchronized and controlled using a National Instruments data acquisition card (NI PCIe-6353) and custom software³⁰.

To enable multimodal imaging, we first sequentially hybridized fluorescent readout probes and then imaged the targeted genomic loci and then the targeted mRNAs. Then we carried out a series of antibody stains and imaging. Specifically, the following protocol was used in order: 118 rounds of hybridization and chromatin tracing imaging, sequentially targeting the 354 chromosome 14 loci using three-colour imaging; 16 rounds of hybridization and MERFISH imaging, combinatorially targeting 298 genes using three-colour imaging; 3 rounds of sequential staining for 6 different antibodies using two-colour imaging.

The protocol for each hybridization included the following steps: incubate the sample with adaptor probes for 30 min for DNA imaging or 75 min for RNA imaging at room temperature; flow wash buffer and incubate for 7 min; incubate fluorescent readout probes (one for each colour) for 30 min at room temperature; flow wash buffer and incubate for 7 min; flow imaging buffer. The imaging buffer was prepared as described previously³⁰ and additionally included 2.5 μg ml^−1 4′,6-diamidino-2-phenylindole (DAPI).

Following each hybridization, the sample was imaged, and then the signal was removed by flowing 100% formamide for 20 min and then re-equilibrating to 2× SSC for 10 min.

RNA MERFISH measurement of the 298-gene panel was carried out with an encoding scheme of 48-bit binary barcode and a Hamming weight of 4. Therefore, in each hybridization round, about 75 adaptor probes were pooled together to target a unique subset of the 298 genes. The genes targeted in each round of hybridization are highlighted in Supplementary Table 9.

Immunofluorescence staining

Antibody imaging was carried out immediately after completing the DNA and RNA imaging. The sample was first stained for 4 h at room temperature using two primary antibodies of two different species (mouse and rabbit), washed in 2× SSC for 15 min and then stained for 2 h using two secondary antibodies for each target species conjugated with fluorescent dyes. Supplementary Table 10 lists all of the antibodies used in this study.

MERFISH image acquisition

We imaged approximately 650 fields of view covering the HPC. After each round of hybridization, we acquired z-stack images of each field of view in four colours: 750 nm, 647 nm, 560 nm and 405 nm. Consecutive z-sections were separated by 300 nm and covered 15 μm of the sample. Images were acquired at a rate of 20 Hz.

Processing of snm3C-seq3 data

Sequencing reads were first demultiplexed by matching the first 8 bp of R1 reads to the predefined well barcodes (https://github.com/luogenomics/demultiplexing). Demultiplexed reads were trimmed to remove sequencing adaptors using Cutadapt 1.18 with the following parameters in paired-end mode: -f fastq -q 20 -m 50 -a AGATCGGAAGAGCACACGTCTGAAC -A AGATCGGAAGAGCGTCGTGTAGGGA. Then 18 bp and 10 bp were further trimmed from the 5′- and 3′-end of the R1 reads, respectively; and 10 bp were trimmed from both the 5′- and 3′- ends of the R2 reads. snm3C-seq3 reads were mapped to the hg38 reference genome using a modified Taurus-MH package (https://github.com/luogenomics/Taurus-MH)³. Briefly, each read end (R1 or R2) was mapped separately using Bismark with Bowtie1 with read1 as complementary (always G to A converted) and read2 (always C to T converted) as the original strand. After the first alignment, unmapped reads were retained and split into three pieces by 40 bp, 42 bp and 40 bp resulting in six subreads (read1 and read2). The subreads derived from unmapped reads were mapped separately using Bismark Bowtie1. All aligned reads were merged into BAM using the Picard SortSam tool with query names sorted. For each fragment, the outermost aligned reads were chosen for the chromatin conformation map generation. The chromatin contacts with both ends mapped to the same positions were considered duplicates and removed for further analysis. Duplicated reads were removed from BAM files using the Picard MarkDuplicates tool before the generation of allc files using the Allcools bam-to-allc tool (https://lhqing.github.io/ALLCools/)²⁶.

Single-molecule fluorescence in situ hybridization

Single-molecule fluorescence in situ hybridization was carried out according to the RNAscope manual (multiplex details). Sequences of target probes, preamplifiers, amplifiers and label probes are proprietary and commercially available (Advanced Cell Diagnostics (ACD)). Typically, the probes contain 20 ZZ probe pairs (approximately 50 bp per pair) covering 1,000 bp. Here we used probes against human genes as single-plex probes, outlined below: Hs-MEF2C (452881), Hs-GAD1-C2 (404031-C3), Hs-RBFOX3-C2 (415591-C2), Hs-TLL1-C3 (439211), Hs-TRPS1 (831611-C3), Hs-PROX1 (530241), Hs-ALDH1L1-C3 (438881-C3), Hs-LRIG1-C2 (407421-C2).

Slides were dried at 60 °C for 1 h and fixed in 4% PFA for 2 h. After several washes in PBS, slides were treated with ACD hydrogen peroxide for 10 min and then washed in water twice before treatment in 1× target retrieval buffer (ACD) for 5 min (at 95–100 °C). After being washed in water and then 100% alcohol, the slides were baked at 60 °C for 30 min. After moistening samples with water, protease treatment was carried out for 15 min at 40 °C in a HybEZ oven. Hybridization of probes and amplification was carried out according to the manufacturer’s instructions. In short, tissue sections were incubated in the desired probe (2–3 drops per section) for 2 h at 40 °C in the HybEZ oven. The slides were washed twice in 1× wash buffer (ACD) for 2 min each and incubated in 5× SSC at room temperature overnight. Amplification and detection steps were carried out using the Multiplex kit (ACD, 320293) for single-plex probes. The following was carried out in repeated cycles for each probe. About four drops of AMP x-FL were added to entirely cover each section and the slide was placed in the HybEZ oven. The slide was incubated for 30 min at 40 °C. Slides were removed from the HybEZ slide rack, and excess liquid was removed before being submerging them in a Tissue-Tek staining dish filled with 1× wash buffer. Slides were washed in 1× wash buffer for 2 min at room temperature. The next AMP x-FL treatment was added, and the cycle was repeated. Slides were washed in PBST, incubated with DAPI for 30 s at room temperature, and mounted in aqua mount (Lerner). Images were taken using a 100× objective on a Leica Stellaris confocal microscope.

Single-cell bimodal data quality control and preprocessing

Cells were filtered on the basis of several metadata metrics: mCCC level < 0.03; global CG methylation level > 0.5; global non-CG methylation level < 0.2; and total 3C interactions > 100,000. Methylation features were calculated as fractions of methylcytosine over total cytosine across gene bodies ±2 kb flanking regions and 100-kb bins spanning the entire genome. Methylation features were further split into CG and CH methylation types. These features were then filtered on mean coverage of ≥10 and values with coverage of <5 were imputed as the mean feature value by sample. Principal component analysis (PCA) was then run using Scanpy⁵⁰ default parameters followed by k-nearest neighbours using only the top 20 principal components by the amount of variance explained and k = 15. Iterative clustering was then carried out with a combination of Leiden unsupervised clustering and UMAP dimensionality reduction, identifying clusters as cell types by marker gene body CH and CG hypomethylation. We observed certain batch effects in our dataset that are associated with the time the data were generated. Harmony²² was used on metadata features to mitigate batch effects occurring between samples in the principal-component feature space. The developmental trajectories of cortical and hippocampal cell types were reconstructed using shared CG methylation feature patterns at cell-type marker genes and integration of cells derived from different age groups using Harmony²². Cells were first separated by their L2 (major cell-type groups) annotation using the shared marker gene approach, and then Harmony integration by pairwise ages for all L2 groups was used to link L3 cell types across ages.

Integration of snm3C-seq3 data with datasets for single-nucleus RNA-sequencing or single-nucleus assay for transposase-accessible chromatin with sequencing

The single-nucleus RNA (snRNA) and prenatal snm3C-seq CG methylation data were co-embedded by inverting the sign of the methylation matrix (owing to an inverse correlation of gene expression to CG methylation). PCA was then applied on the combined data, and Harmony²² was used to correct for the systematic differences between the two modalities. The co-embedded UMAP was then generated from the k-nearest neighbours graph with k = 20 using the top 20 principal components. Annotation transfer was carried out for each cell in the CG methylation data by taking the cell’s top snRNA neighbour and assigning the RNA label to the CG methylation cell. The Jaccard index was then computed on the CG methylation annotation versus the snRNA liftover annotation. The correlation of features between the two modalities was calculated by taking the k = 1 nearest neighbour for all methylation cells and computing the Pearson correlation of raw CG methylation fraction to log-scaled gene expression counts for all genes across all paired cells. Generalized annotations were then made for the co-embedding by running Leiden unsupervised clustering and naming the clusters by their most representative cell type to assess the relative quantity of similar cells in each dataset. The co-embedding of the single-nucleus assay for transposase-accessible chromatin (ATAC) and snm3C-seq CG methylation data followed the same procedure as the snRNA co-embedding, but the Pearson correlation was computed on the raw ATAC gene activity scores versus raw CG methylation fractions.

Pseudotime analysis

Pseudotime analysis was run following the methods outlined in ref. ²⁴. Each pseudotime analysis had clustering preprocessing steps, PCA, k-nearest neighbours with k = 15 using 20 principal components, and Leiden, recomputed for its respective subset of the data. The computed Leiden clusters were then used to initialize a partition-based graph abstraction. This partition-based graph abstraction is used as the precomputed initialization coordinates for the visualization with force-directed graph drawing by the ForceAtlas2 package⁵¹. A root node is then set in the Leiden cluster furthest from the adult cell types, and Scanpy’s implementation of diffusion-based pseudotime was used. In multimodal pseudotimes the same cell is set as the root node in each modality. Genes are selected for display compared to the pseudotime scores by sorting by correlation and anticorrelation to the pseudotime score as well as requiring the 3CGS to have variance of >0.1 and gene length of >90 kb. For Fig. 2n and Supplementary Fig. 4j,u, gene examples were selected by highest gene body CG methylation correlation to the pseudotime and 3CGS anticorrelation. For Fig. 2o and Supplementary Fig. 4k,v, gene examples were selected by highest gene body CG methylation anticorrelation to the pseudotime and 3CGS correlation to the pseudotime. Distribution comparisons are computed by the Wilcoxon rank-sum test.

DMR and TF-binding motif analysis

All CG methylation DMRs were identified from pseudobulk allc files using Methylpy (https://github.com/yupenghe/methylpy)⁵². DMRs identified from a multi-sample comparison of all cell types were used for analyses in Fig. 5, as well as disease heritability enrichment analyses. Trajectory-DMRs were identified using pairwise comparisons of adjacent development stages of a cell-type trajectory. Branch-DMRs were identified using multi-sample comparisons, including the mother cell population from an earlier developmental stage and daughter populations from a later developmental stage. TF-binding motif enrichment analysis was carried out similarly to previously described^17,25,53. DMR regions were lifted over to the hg19 reference genome for the TF-binding motif enrichment analysis. TF-binding motif position weight matrices were obtained from the MEME motif database and scanned across the human hg19 reference genome to identify hits using FIMO (–output-pthresh 1E-5, — max-stored-scores 500,000 and –max-strand)^54,55. DMRs were extended 250 bp both upstream and downstream for overlapping with TF-binding motif hits. The overlap between TF-binding motif hits and DMRs (extended ±250 bp) was determined by requiring a minimum of 1-bp overlap. The enrichment of TF-binding motifs in DMRs was assessed using DMRs (extended 250 bp from centre) identified across adult human tissues (tissue DMRs) as the background⁵². The overlaps between TF-binding motif hits and the foreground DMR list were compared to the overlaps between TF-binding motif hits and tissue DMRs (background) using the hypergeometric test (MATLAB hygecdf).

Single-cell embedding based on chromatin contact

Single-cell contact matrices at 100 kb resolution were imputed by scHiCluster⁵⁶ with pad = 1. The imputed contacts with distance of >100 kb and <1 Mb were used as features for singular value decomposition dimension reduction. Principal components were normalized by singular values and L2 norms per cell and then used for k-nearest neighbour graph construction (k = 25) and UMAP. A total of 25 dimensions were used for the full dataset (Fig. 1 and Extended Data Fig. 2b,e), 20 dimensions were used for the RG subtypes (Extended Data Fig. 2f) and 10 dimensions were used for the MGE or CGE lineage (Extended Data Fig. 2d).

Chromatin loop, differential loop and SIP analysis

Chromatin loops were identified with scHiCluster⁵⁶ for each cell type identified in this study. To identify loops from a group of cells, single-cell contact matrices at 10 kb resolution were imputed by scHiCluster with pad = 2 for the contacts with a distance less than 5.05 Mb (result denoted as Q_cell). We carried out loop calling only between 50 kb and 5 Mb, given that increasing the distance leads to only a limited increase in the number of statistically significant loops. For each single cell, the imputed matrix of each chromosome (denoted as Q_cell) was log-transformed, Z-score-normalized at each diagonal (result denoted as E_cell) and a local background between >30 kb and <50 kb was subtracted (result denoted as T_cell), as in SnapHiC⁵⁷. We then generated pseudobulk matrices for each sample by taking the average across single cells. To compute the variance of each matrix across single cells in loop and differential loop analysis, for each pseudobulk sample, we saved both the mean and mean of squares. Specifically, six pseudobulk matrices were generated as \({Q}_{{\rm{bulk}}}=\sum {Q}_{{\rm{cell}}}/{n}_{{\rm{cell}}}\), \({Q2}_{{\rm{bulk}}}=\sum {Q}_{{\rm{cell}}}^{2}/{n}_{{\rm{cell}}}\), \({E}_{{\rm{b}}{\rm{u}}{\rm{l}}{\rm{k}}}=\sum {E}_{{\rm{c}}{\rm{e}}{\rm{l}}{\rm{l}}}/{n}_{{\rm{c}}{\rm{e}}{\rm{l}}{\rm{l}}}\), \({E2}_{{\rm{b}}{\rm{u}}{\rm{l}}{\rm{k}}}=\sum {E}_{{\rm{c}}{\rm{e}}{\rm{l}}{\rm{l}}}^{2}/{n}_{{\rm{c}}{\rm{e}}{\rm{l}}{\rm{l}}}\), \({T}_{{\rm{bulk}}}=\sum {T}_{{\rm{cell}}}/{n}_{{\rm{cell}}}\), \({T2}_{{\rm{bulk}}}=\sum {T}_{{\rm{cell}}}^{2}/{n}_{{\rm{cell}}}\). A pseudobulk-level t-statistic was computed to quantify the deviation of E and T from 0 across single cells from the cell group, with larger deviations representing higher enrichment against the global (E) or local (T) background. E_cell is also shuffled across each diagonal to generate E_shufflecell, then minus local background for T_shufflecell, and further merged into pseudobulks E_shufflebulk, E2_shufflebulk, T_shufflebulk and T2_shufflebulk, to estimate a background of the t-statistics. E2_shufflebulk is defined as \({E2}_{{\rm{shufflebulk}}}=\sum {E}_{{\rm{shufflecell}}}^{2}/{n}_{{\rm{shufflecell}}}\). T2_shufflebulk is defined as \({T2}_{{\rm{shufflebulk}}}=\sum {T}_{{\rm{shufflecell}}}^{2}/{n}_{{\rm{shufflecell}}}\). An empirical false discovery rate (FDR) can be derived by comparing the t-statistics of observed cells versus shuffled cells. We required the pixels to have average E of >0, fold change of >1.33 against doughnut and bottom left backgrounds, fold change of >1.2 against horizontal and vertical backgrounds⁵⁷, and FDR of <0.01 compared to global and local backgrounds.

Differential loops were identified between age groups in the same major lineage. The detailed analysis framework is shown at https://zhoujt1994.github.io/scHiCluster/hba/loop_majortype/intro.html. To compare the interaction strength of loops between different groups of cells, analysis of variance or a Kruskal–Wallis test can be used. This test is more generalizable, as it does not require the data to be normally distributed. However, in practice, it is very expensive computationally to enumerate through all cells and all loops to run the tests. Therefore, we adopt an analysis of variance framework to compute the F statistics for each loop identified in at least one cell group using either Q_cell (result denoted as F_Q) or T_cell (result denoted as F_T). This analysis requires only Q_bulk, Q2_bulk, T_bulk and T2_bulk for each pseudobulk sample to capture the variability across cells rather than the matrices of each single cell, which makes it feasible across thousands of cells and millions of pixels. We log-transformed and then Z-scored F_Q and F_T across all of the loops being tested and selected the ones with both F_Q and F_T > 1.036 (85th percentile of standard normal distribution) as differential loops. The threshold was decided by visually inspecting the contact maps as well as the correlation of interaction and loop anchor CG methylation. These thresholds selected the top ≈5% loops as differential for downstream analyses.

To identify SIPs, annotated transcription start sites from GENCODE v33 annotation were intersected with chromatin loops, allowing a maximum distance of 5 kb. To determine a threshold of cumulative loop score for SIPs, the cumulative loop score of promoters was modelled by a half-Gaussian distribution with the mean equal to 0 and standard deviation equal to the standard deviation of cumulative loop scores. The threshold for SIP was selected with a P value of 0.001.

Identification of domains and differential domain boundaries

Single-cell contact matrices at 25 kb resolution were imputed by scHiCluster⁵⁶ with pad = 2 for the contacts with distance less than 10.05 Mb. Domains were identified for each single cell. Insulation scores were computed in each cell group (major type or major type in a brain region) for each bin with the pseudobulk imputed matrices (average over single cells) and a window size of 10 bins. The boundary probability of a bin is defined as the proportion of cells having the bin called as a domain boundary among the total number of cells from the group.

To identify differential domain boundaries between n cell groups, we derived an n × 2 contingency table for each 25-kb bin, in which the values in each row represent the number of cells from the group that has the bin called as a boundary or not as a boundary. We computed the chi-square statistic and P value of each bin and used the peaks of the statistics across the genome as differential boundaries. The peaks are defined as the local maximum of chi-square statistics with an FDR of <1 × 10^−3 (Benjamini and Hochberg procedure). If two peaks are within less than five bins of each other, we kept only the peak with a higher chi-square statistic. We also require the peaks to have a Z-score-transformed chi-square statistic of >1.960 (97.5th percentile of standard normal distribution), fold-changes between maximum and minimum insulation score of >1.2, and differences between maximum and minimum boundary probability of >0.05.

3CGS

3CGS is defined as the sum of off-diagonal values from the row and column of the transcription start site bin to the row and column of the transcription end site bin in the imputed 10-kb contact matrices.

Correlation between loop strength and CG methylation

Differential loops for cell-type trajectories were used in these analyses (Supplementary Table 6), with each trajectory analysed separately. Single cells were grouped into meta-cells ordered by the developmental pseudotime scores to boost the power of correlation analyses. Each meta-cell is composed of 20 single cells. Specifically, the joint embedding of CG methylation and chromosome conformation for each trajectory was used to generate a k-nearest neighbour (k = 20, self-included) graph of single cells. Each cell was merged with its other 19 nearest neighbours to generate a meta-cell. To avoid meta-cells that are highly similar to each other, we first computed the number of shared cells between each pair of meta-cells and removed highly similar meta-cells. We repeated this process until no pairs of meta-cells shared more than five cells. To further alleviate the bias towards large homogeneous cell populations, we downsampled the number of meta-cells that originated from RG-1 to half of the number.

The methylation level of a meta-cell at each 10-kb bin was computed by the sum of methylated basecalls divided by the sum of total basecalls over its 20 composed single cells. The average methylation levels of the two anchors were used to correlate with the loop interaction strength. The loop interaction strength of a meta-cell was the average of imputed interactions over the 20 composing single cells. The diffusion pseudotime of a meta-cell was also the mean across its composing single cells. The meta-cells were ordered according to the pseudotime, and the cross-correlation was used to measure the temporary discrepancy between chromatin conformation and DNA methylation changes during development. Intuitively, the cross-correlation measures the correlation between loop strength and CG methylation after shifting the DNA methylation values along the developmental axis for a certain distance. We used the argument of the minima of the cross-correlation for the loops having a negative correlation between interaction and CG methylation to evaluate their timing differences. If the loop interaction is changing before DNA methylation, the DNA methylation values need to be moved backward to maximize its (absolute) correlation with genome structure, therefore corresponding to a negative shift in our measurement. A left-skewed distribution indicates that interaction changes earlier than methylation and vice versa.

Distribution of the distance between interacting loci analysis

To count the number of cis (intra-chromosomal) contacts in each cell and bulk Hi-C data²⁸, we divided the contacts into 143 logarithmic bins, the first of which was for contacts that were separated by less than 1 kb. Each subsequent bin covered an exponent step of 0.125, using base 2. Contacts in bins 1–37 were determined to be noisy and were eliminated, leaving bins 38–141 as valid bins.

The following metrics were used for the following analysis: percentage median, the percentage of contacts in bins 38–89 out of all valid bins; percentage long, the percentage of contacts in bins 90–141 out of all valid bins.

Cells were clustered by the distribution of their distance between interacting loci (k-means, k = 10) and reordered by the average value of log₂[percentage median/percentage long] of each cluster (Fig. 3a–d and Extended Data Fig. 4a–c).

Each cell was assigned to a group by the following criteria (Fig. 3e,f and Extended Data Fig. 4d–k): SE, log₂[percentage median/percentage long) > 0.4; INT, −0.4 > log₂[percentage median/percentage long) > 0.4; LE, log₂[percentage median/percentage long) < −0.4.

To find the clusters enriched in each cell type, we first calculated the percentage belonging to each cluster by cell type (Extended Data Fig. 4d). The enrichment score was obtained by normalizing the fraction of each cell type by the relative cluster sizes. For each cell type, an average of log₂[percentage near/percentage long] scores (representing the ratio of median to long-range interactions) for all individual cells of that type was computed to compare with compartmentalization metrics.

We have computed empirical P values to determine the significance of k-means clusters showing distinct distributions of chromatin contact distances (Fig. 3a and Extended Data Fig. 4c). For each pair of clusters, we randomly selected one cell from each cluster. This process was repeated 1,000 times for each cluster pair. Cell pairs whose chromatin contract distribution showed a Pearson’s correlation coefficient greater than 0.8 were considered similar. For each pair of clusters, we counted the number of times (out of 1,000) that the correlation coefficient exceeded the threshold. This gave an estimate of the similarity in cell patterns between the two clusters.

Chromatin compartment detection and analysis

Pseudobulk .cool files for all cell types were generated and balanced with cooler v0.8.3 at 1 Mb resolution⁵⁸. For each cell type, compartments were assigned with cooltools v0.5.1, through eigenvector decomposition of each chromosome’s cis-interaction matrix⁵⁹. Each 1-mb genomic bin (excluding chromosomes X, Y and M) was assigned to the A or B compartment by the sign of its chromosome’s eigenvector that has the highest Pearson correlation with GC content. A positive sign indicates bin membership in the A compartment; a negative sign indicates bin membership in the B compartment. We use the magnitude of the eigenvector value as the strength of compartment assignment. Using cooltools, we generated saddle plots, which visualize the distribution of observed/expected (O/E) contact frequency between genomic bins stratified by their eigenvector value. For pseudobulk files with more than 30 million contacts, we subset each matrix’s bins whose assignment strengths are in the top 20th percentile for their compartment. Then, we find the sum of O/E values of AA, BB, AB and BA interactions between these bins. For computing AA or BB O/E interaction dominance, we find the fraction of O/E signal explained by these AA or BB interactions, respectively, out of the total O/E signal for the pseudobulk matrix (Extended Data Fig. 6b). Similarly, the formula for compartmentalization strength score is: (sum(AA O/E) + sum(BB O/E))/(sum(BA O/E) + sum(AB O/E)) (Extended Data Fig. 6a,e). When computed on pseudobulk files with more than 30 million contacts, compartmentalization strength scores had no significant correlation with total pseudobulk contacts (Pearson correlation = 0.18, P = 0.10). Two-sided Mann–Whitney U-tests were used to compare distributions of these metrics between groups of cell types (Extended Data Fig. 6b,e).

Differential compartments were identified across all age groups for each major (L2) lineage using dcHiC v2.1 at 100 kb resolution (adjusted P < 0.01)⁶⁰. These results were used to identify 100-kb genome bins that transitioned from the A compartment to the B compartment (AB transition) or vice versa (BA transition) between the earliest and latest ages in each lineage (Extended Data Fig. 6h). For each lineage, the transitioning bins’ CG methylation levels were computed at each age and normalized by subtracting the CG methylation level at the earliest age. The distribution of CG methylation levels for AB versus BA transitions at each age and lineage was compared with two-sided Mann–Whitney U-tests (Extended Data Fig. 6j–l).

Polygenic heritability enrichment analysis

Polygenic heritability enrichment of DMRs and/or chromatin loops was analysed using a stratified linkage disequilibrium score regression (S-LDSC)-based partitioned heritability approach⁶¹. The genome-wide association study (GWAS) summary statistics included in this study were as follows: schizophrenia⁴⁰, bipolar disorder⁶², major depressive disorder⁶³, attention deficit hyperactivity disorder⁶⁴, autism spectrum disorder⁶⁵, Alzheimer’s disease⁶⁶ and height from the UK Biobank⁶⁷ (downloaded from https://alkesgroup.broadinstitute.org/sumstats_formatted/). For each cell type, binary annotations were created using DMR and/or chromatin loop. We considered two types of genomic region—DMR: including all DMRs for a given cell type; loop-connected DMR: including the subset of DMRs that overlap with any of the chromatin loop-called in the matching cell types. To create binary annotations, SNPs in these genomics regions were assigned as 1 and otherwise 0. Then we assessed the heritability enrichment of each of these annotations conditional on the ‘baseline model’³⁹. We reported heritability enrichment and proportion of heritability using Enrichment, Enrichment_std_error, Prop._h2, Prop._h2_std_error columns in S-LDSC results. To assess statistical significance for heritability enrichment differences across annotations (for example, differences between cell types in a developmental trajectory), we used a t-test to test the differences of heritability enrichment of two cell types with d.f. = 200 + 200 − 2, in which 200 corresponds to the number of jackknife samples in the S-LDSC block jackknife procedure.

Overlap between fine-mapped variants and DMR and/or chromatin loop for schizophrenia

We used statistical fine-mapping results that were previously performed in the latest PGC schizophrenia study⁴⁰. We filtered for autosomal high-confidence putative causal SNPs with posterior inclusion probability of >10%, and retained 190 independent association loci (containing 569 SNPs in total), with each loci containing a credible set with 3.0 SNPs on average. We used Fisher’s exact test to assess the overlap between these 569 fine-mapped SNPs and DMR and/or chromatin loop annotations using all SNPs in the GWAS summary statistics as the background (see above for constructing DMR and/or chromatin loop annotations). We reported odds ratios of the overlap. We also assessed the overlap between 190 schizophrenia fine-mapped loci (as aggregates of 569 putative causal SNPs) and DMR and/or chromatin loop annotations (Fig. 5i). We define the overlap between fine-mapped loci and DMR and/or chromatin loop annotations on the basis of whether any high-confidence putative causal SNP in the fine-mapped loci is located in the annotation. Furthermore, we overlapped putative causal SNP and DMR and/or chromatin loop annotations to Genotype-Tissue Expression high-confidence fine-mapped cis-expression quantitative trait locus (eQTL) data (downloaded from https://www.gtexportal.org/home/downloads/adult-gtex/qtl): we first identified SNP–gene pairs such that the putative causal SNP is located in DMRs and connected to the transcription start site of any gene through chromatin loops, and then we overlapped these SNP–gene pairs with cis-eQTL–eGene pairs.

Chromatin tracing analysis

Localization of fluorescent spots

To calculate fluorescent spot localizations for chromatin tracing data, we carried out the following computational steps: we computed a point spread function for our microscope and a median image across all fields of view for each colour channel based on the first round of imaging to be used for homogenizing the illumination across the field of view (called flat-field correction); to identify fluorescent spots, the images were flat-field-corrected, deconvoluted with the custom point spread function, and then local maxima were computed on the resulting images. A flat-field correction was carried out for each colour channel separately.

Image registration and selection of chromatin traces

Imaging registration was carried out by aligning the DAPI channel of each image from the same field of view across imaging rounds. First, the local maxima and local minima of the flat-field-corrected and deconvolved DAPI signal were calculated. Next, a rigid translation was calculated using a fast Fourier transform to best align the local maxima or minima between imaging rounds.

Nuclear segmentation was carried out on the DAPI signal of the first round of imaging using the Cellpose algorithm⁶⁸ with the ‘nuclei’ neural network model. Following image registration, chromatin traces were computed from the drift-corrected local maxima of each imaged locus as previously described³⁰.

RNA MERFISH analysis

The MERFISH decoding followed a similar strategy to that of the MERlin algorithm but operated on spots identified in the images rather than individual pixels. Briefly, the drift- and chromatic-aberration-corrected local maxima (spots) were grouped into clusters, with each cluster containing all spots from all imaging rounds in a 2-pixel radius of an anchor spot. Clusters were generated for every possible anchor spot. Any cluster containing spots from at least four images was then assigned a gene identity by best matching the MERFISH codebook. Each cluster was ranked by the average brightness and the interdistance between the contained spots. These measures were used to filter the decoded cluster and best separate the more confident spots from the less confident.

Protein density quantification

The antibody images were flat-field-corrected, deconvolved and then registered to the chromatin traces using the DAPI signal as described previously. For each chromatin trace, the fluorescent signal of each antibody was sampled at the 3D location corresponding to each genomic locus.

Refinement of the nuclear segmentation

Refined 3D nuclear segmentation was carried out using the 3D Cellpose ‘nuclei’ model based on the NUP98 fluorescent stain.

Integration of snmC-seq3 and RNA MERFISH

RNA profiles generated by MERFISH were unbiasedly clustered using the Leiden method and annotated on the basis of known marker genes shown in Supplementary Fig. 5a. The MERFISH transcriptomic data were then integrated with the data from mid-gestation HPC snmC-seq3 methylation samples by subsampling a set of 220 shared genes, inverting the methylation matrix and using Harmony²² to correct the systematic differences between the two modalities. RG-1 and RG-2 annotations were generated by label transfer from the methylation data using k = 9 nearest-neighbour majority voting in the batch-corrected UMAP space. Label transfer was blocked in mature cell types (dentate gyrus, CA1 and CA2–3) and those not derived from hippocampal RG cells (ependyma and choroid plexus).

MERFISH data availability

Raw imaging data will be provided on request owing to the extraordinary file sizes. Processed data are available at the Gene Expression Omnibus under the accession number GSE213950 as a scanpy.h5ad file. The main ‘.X’ matrix of the object contains log-normalized counts. The full contents of the scanpy object are described below. For brevity, standard contents added by scanpy (for example, connectivities and distances added by sc.pp.neighbors) are not listed.

• obs

  • volm: total pixel volume of the cell based on DAPI segmentation;

  • x_um_abs, y_um_abs: global x and y coordinates of the cell in micrometres;

  • zc, xc, yc: pixel coordinates of the cell centre relative to the field of view;

  • Leiden: unsupervised Leiden clustering;

  • L1: excitatory versus inhibitory;

  • dpt_pseudotime: pseudotime calculated from RG-1;

  • Final_anno_v3: annotation used in figures;

  • Hpc_regional: spatial subset of cells restricted to the HPC;

  • hpcRG: RG-1 and RG-2 annotation in this zone;

  • Fimbria_regional: spatial subset of cells restricted to the HPC;

  • fimbriaRG: RG-1 and RG-2 annotation in this zone;

  • Ventricular_regional: spatial subset of cells restricted to the ventricular zone;

  • ventricularRG: RG-1 and RG-2 annotation in this zone;

  • Refined_volume: Recalculated cell volume based on Nup98 antibodies.

• obsm

  • X_fov: the field-of-view identifier each cell was imaged in;

  • X_raw: raw count matrix;

  • X_spatial: the spatial coordinates of the cells;

  • blank: the count of each blank barcode per cell;

  • X_h_score: a csr sparse matrix containing chromatin trace results. The matrix should be reshaped to 50,374 ×4 × 354 × 5, representing the number of cells, maximum number of homologues, number of chromatin regions, and the z, x, y coordinates followed by the brightness and score of the fluorescent spot. Missing data (that is, containing fewer than four homologues or missing regions) are filled with 0 s;

  • H3 K9 trimethylation, Pol2PSer2, SRSF2, K27Ac, LAMA1, NUP98: antibody signals localized at each chromatin region. Stored as a csr sparse matrix and can be reshaped to 50,374 × 4 × 354, similar to X_h_score.

Reporting summary

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