Single-cell epigenomics reveals mechanisms of human cortical development

Abstract

During mammalian development, differences in chromatin state coincide with cellular differentiation and reflect changes in the gene regulatory landscape¹. In the developing brain, cell fate specification and topographic identity are important for defining cell identity² and confer selective vulnerabilities to neurodevelopmental disorders³. Here, to identify cell-type-specific chromatin accessibility patterns in the developing human brain, we used a single-cell assay for transposase accessibility by sequencing (scATAC-seq) in primary tissue samples from the human forebrain. We applied unbiased analyses to identify genomic loci that undergo extensive cell-type- and brain-region-specific changes in accessibility during neurogenesis, and an integrative analysis to predict cell-type-specific candidate regulatory elements. We found that cerebral organoids recapitulate most putative cell-type-specific enhancer accessibility patterns but lack many cell-type-specific open chromatin regions that are found in vivo. Systematic comparison of chromatin accessibility across brain regions revealed unexpected diversity among neural progenitor cells in the cerebral cortex and implicated retinoic acid signalling in the specification of neuronal lineage identity in the prefrontal cortex. Together, our results reveal the important contribution of chromatin state to the emerging patterns of cell type diversity and cell fate specification and provide a blueprint for evaluating the fidelity and robustness of cerebral organoids as a model for cortical development.

Fig. 1. Single-cell chromatin state atlas of the developing human brain.

a, Schematic cross-section of developing cortex, highlighting major cell types. b, Experimental workflow. c, Uniform manifold approximation and projection (UMAP) of primary scATAC-seq cells (n = 6 individuals, 77,354 cells) coloured by clusters. d, UMAP projection of primary scATAC-seq cells coloured by brain region. Som., somatosensory cortex. e, UMAP projections of gene activity scores for GFAP (marking glia), EOMES (IPCs), DLX1 (interneuron lineage cells), and NEUROD2 (excitatory lineage cells). f, UMAP projection of primary scATAC-seq cells coloured by broad cell type. g, Top, Sankey plot linking scATAC-seq clusters and cell-type predictions. Endo, endothelial; astro/oligo, astrocyte/oligodendrocyte precursors; nEN, newborn excitatory neurons. Bottom left, pile-ups of ATAC-seq signal for each cluster within sets of the top 1,000 enriched peaks for each cluster (Fisher’s exact, two-sided). Pile-ups are centred on peak centres and the ±10-kb flanking region is depicted. Bottom right, significantly enriched transcription factor (TF) motifs for each cluster-specific peak set (hypergeometric test, one-sided). h, Left, predicted enhancer–gene interactions for RGs (pink curves) overlayed with ATAC-seq signal tracks and peaks for RGs, IPCs, nENs, dlENs, and ulENs. Red boxes highlight predicted enhancers of ARX that overlap with validated VISTA forebrain enhancers. Right, LacZ staining marking regions of enhancer activity for the enhancer candidates hs122 and hs145 in embryonic day (E)12.5 mouse embryos, depicting activity in the forebrain^, (images reproduced with permission from VISTA Enhancer Browser; https://enhancer.lbl.gov/). At E12.5, embryos have an average crown-rump length of 8.6 mm. i, Enrichment and depletion of peaks that overlap with promoter-interacting regions, cell-type-specific peaks, and peaks that meet both criteria in copy number variant (CNV) regions enriched in children with NDD (n = 70 NDD-associated CNVs; Fisher’s exact, two-sided, P < 0.05). Asterisks, Bonferroni-corrected significance. Error bars, 95% confidence interval (CI). j, Enrichment and depletion of peaks that overlap with predicted enhancers in promoter and gene regions of genes associated with ASD and NDD, including genes enriched in de novo non-coding mutations (SFARI845; https://gene.sfari.org/database/human-gene/, DDD299, COE253) (Fisher’s exact, two-sided). Asterisks, Bonferroni-corrected significance. k, Heat map of heritability enrichment based on linkage disequilibrium (LD) score regression analysis of genome-wide association study (GWAS) summary statistics in cell-type-specific peak sets coloured by −log₁₀(P). Asterisks, FDR < 0.05. From left to right: Psychiatric Genomics Consortium (PGC) schizophrenia (SCZ) GWAS; an additional PGC schizophrenia GWAS; PGC bipolar (BP) disorder; PGC major depressive disorder (MDD) GWAS; PGC ASD GWAS.

Fig. 2. Dynamic changes in chromatin accessibility during human cortical neurogenesis.

a, Workflow for co-embedding scATAC-seq and scRNA-seq data from the same samples. Left, experimental workflow. Top middle left, UMAP projection of scATAC-seq cells from visual cortex (n = 3 individuals) coloured by leiden clusters. Bottom middle left, UMAP projection of scRNA-seq cells from visual cortex (n = 2 individuals) colored by Leiden clusters. Middle right, UMAP projection of co-embedded cells coloured by assay. Right, UMAP projection of co-embedded scATAC-seq and scRNA-seq cells coloured by Leiden clusters. b, Sankey plot depicting mappings between scATAC-seq clusters, scRNA-seq clusters, and co-embedded clusters. c, Heat map of correlations between scATAC-seq and scRNA-seq clusters based on a set of cell-type marker genes (Methods). d, Left, schematic depicting cell-type marker genes in the cortical excitatory neuronal lineage. Right, projection of log-normalized gene expression and gene activity scores in co-embedded space for SOX2 (RGs), EOMES (IPCs), SATB2 (ulENs), and CRYM (dlENs). e, UMAP projection of co-embedded cells coloured by pseudotime with principal graph overlaid. f, Heat map depicting the average proportion of cells with peaks that are differentially accessible across pseudotime (n = 25,415). Cells are binned by pseudotime into ten equally sized bins. g, Peak accessibility for four individual peaks across ten pseudotime bins with regression line overlaid. h, Predicted enhancer–gene interactions (pink curves) overlaying ATAC-seq signal tracks and peaks with each of the four enhancers in g highlighted in red. i, Heat maps depicting gene expression (left) and gene activity scores derived from open chromatin (right) for 615 cell-type marker genes. Values are averaged within 20 equally sized bins of pseudotime. j, Comparison of moving averages of normalized gene activity scores (red), gene expression (blue), and motif enrichment (green) across pseudotime for PAX6 (left), EOMES (middle), and MEF2C (right).

Fig. 3. Areal differences in chromatin state of progenitor cells foreshadow the emergence of area-specific types of excitatory neurons.

a, Differentiation trajectories for excitatory neurons from the PFC (left) and V1 (right). b, UMAP projection of PFC and V1 scATAC-seq cells (n = 3 individuals) coloured by cell type predictions. Cells from the excitatory lineage are outlined. IN-STR, striatal interneurons. Suffixes 1–3 denote subclusters from Nowakowski et al.. c, d, UMAP projections of PFC and V1 scATAC-seq excitatory lineage cells coloured by area of origin (c) and pseudotime value (d). e, f, UMAP projections of PFC and V1 scRNA-seq excitatory lineage cells (n = 2 individuals) coloured by area of origin (e) and pseudotime value (f). g, h, Left, PFC and V1 scATAC-seq (g) and scRNA-seq (h) excitatory lineage cells ordered from bottom to top by pseudotime value with PFC–V1 divergence branch point shown (Methods). Cells coloured by EOMES gene activity score (g) or expression (h), highlighting IPCs. Right, schematic illustrating the excitatory neuron differentiation trajectory based on chromatin accessibility, in which PFC–V1 divergence becomes apparent at the level of IPCs (g), or gene expression, in which PFC–V1 divergence is not apparent in IPCs (h). i, Pile-ups of PFC and V1 signal in PFC and V1 differentially accessible (DA) peak sets. Pileups are centred on peaks and show ±10-kb flanking regions. j, Transcription factor motif enrichments of RA-related transcription factors in set of 4,176 PFC-specific peaks (Fisher’s exact, two-sided, FDR < 0.05). k, UMAP projection of deviation scores of motif enrichment for TGIF1. l, Experimental design to test role of RA in organoid area identity. m, UMAP projection of scRNA-seq data from day 70 organoids (n = 11,415 cells). Cells coloured by treatment. VA, vitamin A. n, Left, schematic of expected expression patterns of BCL11B, SATB2, AUTS2, and NR2F1 in primary human cortex. Right, images of primary developing human cortex from the PFC (left) and V1 (right) immunostained for CTIP2 and SATB2 (top) or AUTS2 and NR2F1 (bottom). Representative images shown from n = 2 specimens. o, Left, UMAP projection of cells coloured by expression of NEUROD2, TBR1, SATB2, and NR2F2. Right, images of organoids cultured without (left) or with vitamin A (right) immunostained for CTIP2 and SATB2 (top) or AUTS2 and NR2F1 (bottom). Representative images shown from n = 3 lines.

Fig. 4. Cell type-specific differences in chromatin accessibility between cerebral organoids and the developing human brain.

a, Schematic depicting experimental workflow. iPS cells, induced pluripotent stem cells. b, UMAP projection of all organoid scATAC-seq cells (n = 5 organoids from 3 different lines and 3 different time points; 23,555 cells) coloured by Leiden clusters. (Cluster 16 not depicted, see Extended Data Fig. 10e.) c, UMAP projections of gene activity scores for GFAP marking RGs, EOMES marking IPCs, DLX1 marking interneurons, and NEUROD2 marking excitatory neurons. d, UMAP projection of all organoid scATAC-seq cells coloured by sample. e, Heat map of Pearson correlations between primary and organoid scATAC-seq clusters based on a common peak set. MGE IN, interneurons identified in the MGE. f, Venn diagram of overlap between the full primary peak set (red), a down-sampled primary peak set (blue), and the organoid peak set. g, UMAP projection of enrichment Z-scores of peaks that overlap between primary and organoid datasets on the primary scATAC-seq dataset. h, UMAP projections of enrichment Z-scores of RG-specific peaks (Fisher’s exact, two-sided, FDR < 0.05) in all primary scATAC-seq cells (left) and all organoid scATAC-seq cells (right). i, Left, proportion of cell-type-specific primary peaks in the organoid peak set. Right, proportion of gene-linked enhancers for each cell type in the organoid peak set.

Extended Data Fig. 1. Batch correction and quality control metrics for primary scATAC-seq data.

a, Density curve of fragment size distribution for deduplicated, uniquely mapped fragments that passed quality filters for each sample. b, Plot of promoter ratio against ln(read depth) for all cell barcodes detected. Cells were included in downstream analysis that had a promoter ration of 10–60% and ln(read depth) of 3–5. Red lines indicate upper and lower thresholds. c, Histogram of cellular coverage in 5-kb genomic bins for all cells. Cells with reads in fewer than 500 bins were removed from downstream analysis. d, Heat map of all-by-all Pearson correlations between all samples that had bulk libraries prepared in parallel. Bulk and aggregate single-cell libraries for each sample are included. Correlation was calculated in the space of the merged primary peak set (n = 459,953 peaks). e, Heat map of all-by-all Pearson correlations between all aggregate single-cell libraries. Correlation was calculated in the space of the merged primary peak set (n = 459,953 peaks). f, Heat map of all-by-all correlations between primary single-cell data aggregated by area of origin. Correlation was calculated in the space of the merged primary peak set (n = 459,953 peaks). g, h, UMAP projections of all primary scATAC-seq cells before (g) and after (h) batch correction coloured by sample (see Methods). i, UMAP projection of log(read depth) for deduplicated, uniquely mapped reads that passed quality filters on all primary scATAC-seq cells. j, UMAP projection of fraction of fragments in peaks on all primary scATAC-seq cells. k, UMAP projection of all primary scATAC-seq cells coloured by condition (fresh or frozen). l, Bar plot depicting the proportion of cells from each brain region for each Leiden cluster for all primary scATAC-seq cells. m, Heat map of all-by-all Pearson correlations between clusters. Correlation is calculated in the space of all primary peaks (n = 459,953).

Extended Data Fig. 2. Batch Correction of primary scATAC-seq samples.

a, UMAP projections of all primary scATAC-seq cells that passed quality control before batch correction with all cells from each sample coloured in red. b, UMAP projections of all primary scATAC-seq cells that passed quality control after batch correction with all cells from each sample coloured in red.

Extended Data Fig. 3. Gene activity scores correlate with cell type-specific expression of marker genes.

a, UMAP projections of all primary scATAC-seq cells coloured by gene activity score. From top left to bottom right: NKX2.1 (also known as NKX2-1) marking MGE cells, AQP4 marking glia/astrocytes, TBR1 marking excitatory neurons, FEZF2 marking deep layer excitatory neurons, HES1 marking radial glia, HOPX marking outer radial glia, SATB2 marking upper layer excitatory neurons, CCL4 marking microglia, CRYAB marking truncated radial glia, LHX6 marking MGE-derived interneurons, OLIG1 marking oligodendrocyte precursors, and SOX2 marking radial glia. b, Schematic of the CellWalker algorithm used to assign cell-type labels to scATAC-seq cells based on integration with scRNA-seq data. c, UMAP projection of cell-type labels assigned by CellWalker. d, UMAP projection of radial glia with cell-type assignments from CellWalker. e, UMAP projections of gene activity scores for PAX6 and GLI3, genes that are ubiquitously expressed in radial glia; HOPX, an oRG-specific gene; and CRYAB, a tRG-specific gene. f, UMAP projection of Z-scores of enrichment of oRG- and tRG-specific peaks (Fisher’s exact, two-sided test, P < 0.05). g, UMAP projections of Z-scores of enrichment of area-specific peaks for each area in all primary scATAC-seq cells (Fisher’s exact, two-sided test, P < 0.05).

Extended Data Fig. 4. Annotation of primary scATAC-seq peaks.

a, Distribution of primary scATAC-seq peaks in genomic features. b, Distribution of primary scATAC-seq peaks around transcription start sites. c, Bar plot of log(fold enrichment) of primary scATAC-seq peaks in chromatin states. Chromatin states defined by the 25-state model from Roadmap Epigenomics (see Methods). d, Intersection of cell-type-specific peaks (Fisher’s exact, two sided, FDR < 0.05), predicted enhancer peaks (see Methods), and peaks overlapping promoter-interacting regions identified by H3K4me3 PLAC–seq. e, UMAP projections of Z-scores of enrichment of cell-type-specific peaks (Fisher’s exact, two-sided, FDR < 0.05) for each broad cell type. f, Browser tracks highlighting cell-type-specific predicted enhancers. Left, highlighting predicted enhancers linked to SOX2 in RGs that are not present in ulENs. Right, highlighting a predicted enhancer for GRIN2B that is present in ulENs and not RGs. g, h, Bar plots of −log₁₀(P) of gene ontology biological processes that are enriched in cell-type-specific predicted enhancers of MGE-derived interneurons (g) and RGs (h) (see Methods). i, Heat map of Z-scores of transcription factor motif enrichments of key lineage-associated transcription factors in each cluster.

Extended Data Fig. 5. scATAC-seq peaks overlap with previously annotated bulk ATAC-seq peaks and validated forebrain enhancers.

a–c, Overlap of all primary peaks with the peak set from de la Torre-Ubieta et al. (a), Li et al. (b), and Markenscoff-Papadimitriou et al. (c). Left to right: Venn diagram of overlaps, UMAP projection of Z-scores of enrichment of overlapping peaks in all primary scATAC-seq cells, bar plot of proportions of cell-type-specific peaks (Fisher’s exact, two-sided, FDR < 0.05) present in overlapping set, bar plot of proportions of predicted enhancers for each cell type (see Methods) present in the overlapping set. d, Venn diagram of overlap of all predicted enhancers (see Methods) with gene-linked enhancers from Amiri et al., high confidence enhancers from Wang et al., and putative regulatory elements from Markenscoff-Papadimitriou et al.. e, Venn diagram of overlaps of VISTA forebrain enhancers with all primary scATAC-seq peaks, and all predicted enhancers (see Methods).

Extended Data Fig. 6. Enrichment and depletion of disease associated variants in scATACseq peaks.

a–c, Enrichment and depletion of cell-type-specific peaks (a), peaks that overlap with H3K4me3 PLAC–seq interactions (b) and cell-type-specific peaks that overlap with H3K4me3 PLAC–seq interactions (c) (Fisher’s exact, two-sided, FDR < 0.05) in promoter and gene regions of genes associated with ASD and NDD, including genes enriched in DNMs (SFARI845, DDD29960, COE25361). Asterisks indicate tests that pass Bonferroni significance. d, Enrichment and depletion of DNMs in predicted enhancer peaks for each cell type in ASD probands compared with unaffected siblings. DNM data from a total of 2767 probands and 1855 unaffected siblings were included in the analysis. No tests reached Bonferroni significance. Bars represent 95% CI. e, Enrichment and depletion of DNMs in cell-type-specific peaks for each cell type in ASD probands compared with unaffected siblings. DNM data from a total of 2,767 probands and 1,855 unaffected siblings were included in the analysis. No tests reached Bonferroni significance. Bars represent 95% CI. f, Enrichment and depletion of predicted enhancer peaks in copy number variant (CNV) regions enriched in paediatric cases of NDD (n = 70 NDD-associated CNVs). No tests reached Bonferroni significance. Bars represent 95% CI. g, Enrichment of cell-type-specific enhancers located in TADs with neurodevelopmental disease-associated genes (left) and schizophrenia-associated genes (right). Three distinct sets of TADs were used (top to bottom): TADs defined in the cortical plate (CP) of developing human cortex, TADs defined in the germinal zone (GZ) of developing human cortex, and TADs defined in adult cortex.

Extended Data Fig. 7. Dynamic patterns of gene expression, chromatin accessibility, and transcription factor motif enrichment across pseudotime.

a–c, Heat map of normalized expression (a), normalized gene activity scores (b) and deviation scores of motif enrichments (c) of 40 key lineage transcription factors across pseudotime. Pseudotime is binned into 20 equally sized bins and expression, activity scores and deviations are averaged across all cells in each bin. Deviation scores determined using ChromVAR (Methods). d, e, Heat maps of correlations between gene expression and gene activity scores for 40 key lineage transcription factors (d) and for 615 cell-type marker genes (e).

Extended Data Fig. 8. Chromatin state profiling reveals divergence of PFC and V1 excitatory lineages.

a–d, UMAP projections of scATAC-seq cells from PFC and V1 samples before batch correction coloured by sample (a), after batch correction coloured by sample (b), coloured by Leiden cluster (c) and coloured by gene activity score for EOMES, a marker of IPCs (d). e–h, UMAP projections of scRNA-seq cells from PFC and V1 samples before batch correction coloured by sample (e), after batch correction coloured by sample (f), coloured by Leiden cluster (g) and coloured by gene activity score for EOMES, a marker of IPCs (h). i, j, Projection of NHLH1, PPP1R17, and NEUROD4 gene activity scores on PFC and V1 scATAC-seq cells (i) and scRNA-seq cells (j) ordered by pseudotime with PFC/V1 divergence branch points displayed. k, Volcano plot of peaks that are differentially accessible between PFC and V1 IPCs. Peaks highlighted in red have log-transformed fold change (logFC) > 0.5 and FDR < 0.05 (n = 1,819). l, Volcano plot of genes that are differentially expressed between PFC and V1 IPCs. Genes highlighted in red have logFC > 0.5 and FDR < 0.05 (n = 11).

Extended Data Fig. 9. Modelling the PFC–V1 split in the developing cortex.

a, b, UMAP projections of Z-scores of enrichment of PFC-specific peaks (a; n = 4,176) and V1-specific peaks (b; n = 21,030) in all PFC and V1 scATAC-seq cells (Fisher’s exact, two-sided, FDR < 0.05). c, Top enriched transcription factor motifs in V1-specific peak set as determined by HOMER (hypergeometric test, one-sided). d, UMAP projection of ChromVAR deviation scores of motif enrichment of MEF2C in all PFC and V1 scATAC-seq cells. e, UMAP projection of scRNA-seq data from organoids (n = 3) cultured in the presence of vitamin A, without the presence of vitamin A, and in the presence of DEAB. Cells coloured by cluster. f, Schematic depiction of classifier method used to assign area identity to organoid cells on the basis of defined area-specific gene modules. g, Bar plot depicting proportion of excitatory neurons from each treatment group classified as more PFC-like or more V1-like on the basis of calculation of module eigengene values for area-specific modules (see Methods). Asterisks indicate significant differences in proportions (DEAB: PFC, 1,160/2,622; V1, 1,462/2,622; NoVA: PFC, 563/1,556; V1, 993/1,556; VA: PFC, 1,831/2,976; V1, 1,145/2,976) (χ² test, one-sided, vitamin A versus no vitamin A: P = 3.209 × 10^–59; vitamin A versus DEAB: P = 2.79 × 10⁻³⁸). h, Violin plots depicting expression levels of SATB2, NR2F1, and NR2F2 for excitatory neurons from each treatment group. i, Heat map showing gene expression of differentially expressed genes between excitatory neurons cultured with and without vitamin A. j, Images of organoids cultured with and without vitamin A stained with DAPI and immunostained for NR2F1 and AUTS2. All images taken at 10× resolution. Representative images shown from n = 3 lines.

Extended Data Fig. 10. Comparison of organoid and primary peaks reveal significant differences in chromatin landscapes.

a, Representative image of organoid (1323-4 line depicted) immunostained for FOXG1, a marker of cortical identity. Image taken at 10× resolution. Representative image shown from n = 3 lines. b, Density curves of fragment size distributions for each organoid sample. Fragments are deduplicated, uniquely mapped fragments that have passed quality filters (see Methods). c, UMAP projection of all organoid scATAC-seq cells coloured by cluster. d, UMAP projection of log(read depth). Fragments are deduplicated, uniquely mapped fragments that have passed quality filters (see Methods). e, Bar plot depicting the proportions of cells in each cluster from each organoid sample. f, Venn diagram of overlaps of all primary scATAC-seq peaks, all organoid scATAC–seq peaks, all H3K27ac peaks from Amiri et al., and all ATAC-seq peaks from Trevino et al.. g, UMAP projection of Z-scores of enrichment of peaks that are present in the organoid scATAC-seq dataset but not in the primary scATAC-seq dataset. h, UMAP projection of Z-scores of enrichment of peaks that are present in the primary scATAC-seq dataset but not in the organoid scATAC-seq dataset. i, Genome browser tracks depicted a predicted enhancer of SOX2 that is present in both the full and down-sampled primary RGs but not detected in organoid RGs. j, UMAP projection of scRNA-seq data from the same three organoid lines used for scATAC-seq analysis coloured by line (n = 19,509 cells). k, UMAP projection of scRNA-seq data coloured by cluster. l, UMAP projection of normalized gene expression of PAX6, SOX2, EOMES, and NEUROD2. Maximum value was set at the 99th quantile. m, Bar plot depicting the proportions of cells in each cluster from each organoid sample. n, UMAP projection of normalized gene expression of FOXG1 split by organoid line. Maximum value was set at the 99th quantile.