Multiomic single-cell profiling identifies critical regulators of postnatal brain

Abstract

Human brain development spans from embryogenesis to adulthood, with dynamic gene expression controlled by cell-type-specific cis-regulatory element activity and three-dimensional genome organization. To advance our understanding of postnatal brain development, we simultaneously profiled gene expression and chromatin accessibility in 101,924 single nuclei from four brain regions across ten donors, covering five key postnatal stages from infancy to late adulthood. Using this dataset and chromosome conformation capture data, we constructed enhancer-based gene regulatory networks to identify cell-type-specific regulators of brain development and interpret genome-wide association study loci for ten main brain disorders. Our analysis connected 2,318 cell-specific loci to 1,149 unique genes, representing 41% of loci linked to the investigated traits, and highlighted 55 genes influencing several disease phenotypes. Pseudotime analysis revealed distinct stages of postnatal oligodendrogenesis and their regulatory programs. These findings provide a comprehensive dataset of cell-type-specific gene regulation at critical timepoints in postnatal brain development.

Extended Data Fig. 1 |. Representative FANS sorting for HiC and genotype control for multiome.

a, Sample identity check based on the pairwise comparison of genotypes called from RNA-seq (left) and ATAC-seq (right) reads showing a clear distinction between samples originating from the same donor (each donor was profiled in 4 brain regions) versus different donors. The Kinship scores were calculated by King software (Methods). b, Gating with forward and side scatter was used to identify DAPI+ singlets. c, DAPI+ singlets were gated and assessed for binding of both NeuN and Sox10. This allowed for the identification and subsequent collection of oligodendrocyte (DAPI+ NeuN− SOX10+) and microglia/astrocyte (DAPI+ NeuN− SOX10−) nuclei. In parallel, the NeuN+ population was gated and assessed for binding of Sox6, enabling the isolation of nuclei from both GABAergic (DAPI+ NeuN+ SOX6+) and glutamatergic neurons (DAPI+ NeuN+ SOX6−). d, Representative yield of nuclei from different cell types using the sorting strategy described in Methods.

Extended Data Fig. 2 |. Annotation of human brain single-nucleus multiome dataset.

a, RNA (left) and ATAC (right) representation of the full dataset colored by main cell types. b-c, Cosine similarity correlation on pseudobulk level of RNA between cell types defined in our multiome dataset and cell types from published datasets^,,. d, Geneset enrichment analysis (GSEA) of top 500 genes within highest contribution to gene expression variability based on age, brain region and sex, colored by P-value (one-sided Fisher’s Exact test) and size corresponding to number of genes within Gene Ontology (GO) term category.

Extended Data Fig. 3 |. Disease heritability enrichment patterns for cell type - TF regulator pairs.

Enrichment of common risk variants associated with selected disease traits in regulatory elements (peaks) linked to predicted regulatory targets (genes) of 146 TFs. For each cell type, co-localization was estimated between common risk variants and the peaks linked by ABC approach using LD-score partitioned heritability (Methods). For each cell-type and disease trait, up to 5 of the most specific FDR-significant TF regulators are labeled with their names. SCZ = schizophrenia, PD = Parkinson’s disease, MS = multiple sclerosis, MDD = major depressive disorder, BD = bipolar disorder, ASD = autism spectrum disorder, AN = anorexia nervosa, ALS = amyotrophic lateral sclerosis, AD = Alzheimer’s disease, ADHD = attention deficit hyperactivity disorder). ‘#’: indicates significance for enrichment in MAGMA LDSR calculation of P-value after FDR (Benjamini-Hochberg) correction of multiple testing across all tests in plot (Methods); ‘.’: Nominally significant for enrichment.

Extended Data Fig. 4 |. Comparison of GWAS trait enrichments across cell types between our study and other studies.

The heatmaps and scatter plots provide complementary visualizations of the same data. In the scatter plots, the blue line represents the best-fit regression line from a linear model. For Li et al., enrichments in the original study were reported at the cell subtype level (for example, SST and PVAL interneurons) rather than at the broader cell type level (for example, GABAergic neurons). To facilitate comparison, we ran the comparison against the cell subtype with the lowest P-value. SCZ = schizophrenia, PD = Parkinson’s disease, MS = multiple sclerosis, MDD = major depressive disorder, BD = bipolar disorder, ASD = autism spectrum disorder, AN = anorexia nervosa, ALS = amyotrophic lateral sclerosis, AD = Alzheimer’s disease, ADHD = attention deficit hyperactivity disorder. ‘#’: indicates significance for enrichment at BH-adjusted P-value < 0.05; ‘·’: Nominally significant for enrichment.

Extended Data Fig. 5 |. Pseudotime analysis OPC and OL population.

a, OPC and OL subtype-specific marker gene expression projected on RNA-UMAP. b, Percentage of nuclei corresponding to each brain region across distinct OPC and OL subtypes. c, UMAP projection of RNA-seq and ATAC-seq for the OPC/OL population, colored by OPC and OL subtypes. d, RNA-UMAP colored by normalized pseudotime values from monocle3, palantir and PC1 approach (upper row). Normalized expression patterns of selected OPC and OL subtype markers across pseudotime progression from each method. e, per-nuclei Pearson correlation of pseudotime values across selected methods. f, distribution of pseudotime values from monocle3 across distinct brain regions; ACC (n = 11,585 nuclei), DLPFC (n = 12,226 nuclei), CN (n = 12,933 nuclei) and Hipp (n = 10,489 nuclei). The center line represents the median, the box captures the interquartile range, and the whiskers show the maximum and minimum values within 1.5 times the interquartile range.

Extended Data Fig. 6 |. RNAscope validation of transcriptionally distinct OPC and oligodendrocyte populations.

a-b, Representative 20x stitched images of RNAscope in situ hybridization of dentate gyrus region (n = 2 samples), using probes targeting PDGFRA and ATRNL1 (marking OPC Type I), BMPER (marking OPC Type II) (a), GPR37 (marking OL Type I), and RBFOX1 and ACTN2 (marking OL Type II) (b), and DAPI for nuclear staining, with the scale bar corresponding to 200 μm. c, PDGFRA, ATRNL1 and BMPER (for OPC Type I and II, respectively) gene expression, visualized on the UMAP of OPC and OL populations. d-e, zoom in images showing representative OPC Type I cells co-expressing PDGFRA and ATRNL1 (d), and OPC Type II cells expressing BMPER (e), with scale bars corresponding to 20 μm. f, GPR37 and RBFOX1, ACTN2 (for OL Type I and II, respectively) gene expression, visualized on the UMAP of OPC and OL populations. g-h, zoom in images showing representative OL Type I cells expressing GPR37 (g), and OL Type II cells co-expressing RBFOX1 and ACTN2 (h), with scale bars corresponding to 20 μm.

Fig. 1 |. Multimodal profiling and multivariate exploration of the postnatal human brain.

a, Schematic of study design and computational analyses. 3DG, three-dimensional genome. b, WNN representation of snRNA- and snATAC-seq data projected onto a UMAP with main cell types. c, Accompanying pie charts show the proportion of nuclei represented within main covariates—sex, brain bank, brain region, age group and cell type. MSSM, Mount Sinai. d, Top two gene markers (based on highest average expression) visualized for their gene expression across main cell types. e, Cosine similarity between gene expression and corresponding gene activity (derived from chromatin accessibility) for main cell types. f, Variance explained by covariates on gene expression (left, n = 16,661 genes), chromatin accessibility (middle, n = 181,432 peaks) and cell-type composition (right, n = 10 cell types) using variance partition approach (Supplementary Data 2). Within each violin, the center line indicates the median, the box shows the interquartile range (IQR) and the whiskers indicate the highest/lowest values within 1.5× the IQR. Points outside 1.5 times the IQR from the first and third quartiles are considered outliers.

Fig. 2 |. Cell-type level eGRN analysis in postnatal human brain.

a, Schematic of eRegulon detection approach. b, Heatmap/dotplot showing the three to five most specific eRegulons for each cell type, based on their RSS specificity score. Scaled eRegulon expression is shown as color and RSS is shown as size scale. c, Heritability enrichment of common risk variants associated with selected disease traits in genes controlled by cell-specific TF regulators from b (upper left triangle; calculated by MAGMA) and in enhancers of those genes (bottom right triangle; calculated by S-LDSC) MAGMA enrichment analysis for genes in each cell-type regulator from b for selected traits. Hash sign, significance for enrichment in MAGMA at P value 0.05 after FDR (Benjamini–Hochberg) correction of multiple testing across all tests in plot (Methods); small black squares, significant for enrichment at P value 0.05. NPD, neuropsychiatric diseases; NDD, neurodegenerative diseases. d, Conservation of shared target regions between pairs of selected regulators from b. Proportion of shared peaks in a pairwise manner is shown, with yellow corresponding to the highest value. e, Conceptualization for calculation of footprinting score (top) using the TOBIAS package. Comparison of footprinting scores calculated for microglial TF motifs SPI1, RUNX1, NFATC2, IRF8 and IKZF1 located at the peaks recognized by eGRN analysis as their shared cotargets (orange violin, n = 68) versus peaks recognized by eGRN as targets for just one of them (purple violin, n = 6,175) (bottom). The difference between those groups is significant by the one-sided Kolmogorov–Smirnov test at P value < 10⁻⁸. Within each violin, the center line indicates the median, the box shows IQR and the whiskers indicate the highest/lowest values within 1.5× IQR. Schematic of footprinting was adopted from the original publication with permission from the authors.

Fig. 3 |. Summary of cell-type-specific disease heritability and disease regulome landscape.

a, Heritability enrichment of brain cell types in brain-related diseases calculated from cell-specific marker genes (RNA-seq) and peaks (ATAC-seq). Vertical dashed line denotes threshold for nominal significance (RNA-seq MAGMA/ATAC-seq LDSR regression P value at 0.05). Size of dots denotes statistical significance after applying Benjamini–Hochberg procedure to control the FDR at 0.05. NS, nonsigificant. b, Mean of mutational constraint Z score for peaks grouped by the number of cell types for which the peaks have at least one E–P ABC link. The Z scores were derived from a genome-wide mutational constraint map. c, Relative proportion (left y axis, bars) and absolute number (right y axis, solid line) of explained loci per each GWAS trait across all investigated cell types. d, Histogram of the E–P_ABC-MAX distances between the GWAS risk variant and TSS of the putatively regulated gene, with red and blue dashed lines representing the median and mean values, respectively. e, Histogram of the number of genes ‘skipped’ by E–P_ABC-MAX to reach the putatively regulated gene, with red and blue dashed lines representing the median and mean values, respectively. f, Number of explained loci per cell type across all investigated GWAS traits.

Fig. 4 |. Annotation of GWAS prioritized genes.

a, Overlap between genesets (assessed using a one-sided Fisher’s exact test) representing biological pathways and genes prioritized for each cell type in the MS GWAS. The top ten intersections are displayed, with a vertical dashed dotted line indicating nominal significance at P value 0.05 (one-sided Fisher’s exact test), and dot size reflecting FDR (Benjamini–Hochberg) significance. Only cell types with FDR-significant enrichments are included. b, Localization of prioritized genes within the synaptic hierarchy of the SynGO database. The sunburst plot features the synapse at the center, pre- and postsynaptic locations in the first ring and child terms in subsequent rings. The color scheme in the legend indicates the number of genes in each term. Selected functional terms and genes are labeled. c, Set of 16 genes prioritized by our approach for ALS GWAS. The indication of ‘previous evidence’ is derived from the ALS GWAS. The red scale reflects the ABC score of E–P_ABC-MAX link connecting the gene with GWAS loci. The blue score depicts the log₁₀(E–P_ABC-MAX distance between the enhancer and promoter). d, Normalized snATAC-seq–derived pseudobulk tracks demonstrating the complex cell-specific regulation of the C9orf72 gene with highlighted microglial ALS GWAS-related E–P link. The y axis of each pseudobulk track is normalized and corresponds to signal intensity (with each maximum at 120). Chr, chromosome.

Fig. 5 |. Characterization of oligodendrogenesis by integrative pseudotime analysis.

a, Joint ATAC-seq and RNA-seq (WNN) UMAP representation of OPC and OL population colored by subtype annotation. b, Stacked barplot representation of OPC and OL subtype proportions by age group along with subtype level composition variation across age using crumblr (Methods). The x axis represents the effect size, and the color and dot size denote log₁₀(FDR); plus symbol, statistically significant results; error bars, 95% confidence intervals, calculated from s.e. c, Pseudotime values from monocle3 projected onto RNA-UMAP. d, Boxplots for pseudotime distribution across OPC and OL subtypes (left, n_OPC Type I = 12,445 nuclei; n_OPC Type II = 1,371 nuclei; n_OL Type I = 29,560 nuclei; n_OL Type II = 3,691 nuclei and age groups (right; n_infant = 5,153 nuclei, n_childhood = 14,844 nuclei, n_adolescence = 11,873 nuclei, n_early_adulthood = 3,164 nuclei, n_late_adulthood = 11,749 nuclei)). The center line represents the median, the box captures the IQR, and the whiskers show the maximum and minimum values within 1.5 times the IQR. e, Hotspot gene modules classified into three main gene expression trends along pseudotime based on hierarchical clustering (Supplementary Fig. 17b,c). Up (red), down (green) and constant (gray) trendlines. f, Gene set enrichment analysis (GSEA) for each trendline depicting top enriched Gene Ontology (GO) terms of biological processes, with size of the dot representing gene counts and color representing adjusted P value (FDR, Benjamini–Hochberg correction) (Supplementary Data 15). ER, endoplasmic reticulum; SRP, signal recognition particle.

Fig. 6 |. Identification of key regulatory programs in oligodendrogenesis.

a, Top 20 eRegulons correlated with early and late pseudotime, with barplot showing corresponding gene counts within each eRegulon (left). Heatmap of eRegulon AUCell scores and eRegulon TF-expression values along oligodendrogenesis pseudotime from Fig. 5c, with side margins indicating relative expression levels (cutoff −2,2) (middle). Distribution of AUCell scores of OLa-eRegulons in OPC and OL. b, Relative pairwise overlap between up- and downregulated genes in OPC and OL subtypes with early (light gray) and late (dark gray) stage correlated eRegulons. c, MAGMA enrichment analysis of early- and late-stage OLa-eRegulons for ten selected traits relevant to neuropsychiatric and neurodegenerative diseases. Hash sign, significance for enrichment in MAGMA at P value < 0.05 after FDR (Benjamini–Hochberg) correction of multiple testing across all tests in plot; small black squares, significant for enrichment at P value 0.05.