Single-cell epigenomic reconstruction of developmental trajectories from pluripotency in human neural organoid systems

Cell fate progression of pluripotent progenitors is strictly regulated, resulting in high human cell diversity. Epigenetic modifications also orchestrate cell fate restriction. Unveiling the epigenetic mechanisms underlying human cell diversity has been difficult. In this study, we use human brain and retina organoid models and present single-cell profiling of H3K27ac, H3K27me3 and H3K4me3 histone modifications from progenitor to differentiated neural fates to reconstruct the epigenomic trajectories regulating cell identity acquisition. We capture transitions from pluripotency through neuroepithelium to retinal and brain region and cell type specification. Switching of repressive and activating epigenetic modifications can precede and predict cell fate decisions at each stage, providing a temporal census of gene regulatory elements and transcription factors. Removing H3K27me3 at the neuroectoderm stage disrupts fate restriction, resulting in aberrant cell identity acquisition. Our single-cell epigenome-wide map of human neural organoid development serves as a blueprint to explore human cell fate determination.


Regionalization within the glial lineage
We investigated regional variance also at the chromatin level when we associate region-specific peaks of the active marks H3K27ac and H3K4me3 to the closest gene (Supplementary Fig. 2j).Interestingly, H3K27me3 at regional marker genes (LHX1, FOXG1, LHX2) showed higher gene repression variance in NPCs compared to astrocytes indicating that gene repression in astrocytes is less predictive of regional identity.
Altogether, our analysis reveals that the interplay of chromatin modifications induces and stabilizes brain region diversification events.

Differential regulation of regional diversification and neuronal differentiation
We analyzed each gene across all modalities to determine how histone modifications might be differentially involved in brain region diversification versus neuronal differentiation.(Supplementary Fig. 7a-d, Supplementary Table 16, see Methods Analysis of pseudotemporal and regional variance).Overall, RNA expression and histone modifications (Supplementary Fig. 7e) showed high concordance.Variance in gene expression and H3K27ac enrichment was mostly explained by regional diversification (Supplementary Fig. 7f-g), whereas H3K27me3 showed similar variance across regional branches and over pseudotime.H3K4me3 showed lower variance across regional branches than over pseudotime (Supplementary Fig. 7g).We characterized the genome-wide distribution of branch-specific peaks.This revealed that the majority (~75%) of H3K4me3 peaks was located at promoter regions, while for H3K27ac and H3K27me3, a lower percentage (60%) of branch-specific peaks was found at promoter regions and instead, 40% of these peaks were found in distal regions (Supplementary Fig. 7h).g-i, UMAP embedding colored by expression of genes that show a brain region specific expression pattern in NPCs (e), are specifically expressed in forebrain astrocytes and NPCs (f) are specifically expressed in non-forebrain astrocytes and NPCs (g).j, Scatter plot comparing the regional variance (R 2 ) of histone modification enrichment (for H3K27ac, H3K4me3, H3K27me3) between astrocytes and NPCs.This suggests that gene repression (H3K27me3) is less predictive of regional identity in astrocytes compared to NPCs (see Methods Analysis of pseudotemporal and regional variance for details).

Supplementary Fig. 3 Pseudotemporal reconstruction of dorsal telencephalic neuron differentiation using diffusion maps
a, Schematic of the analysis.We subset the telencephalon branch and isolated only cells of the 120 days time point.The cells were again filtered for quality and pseudotime for all modalities was calculated using diffusion maps (see Methods Reconstruction of the neurogenesis trajectory for the 4 months time point for details).b, Heatmap showing scaled and smoothed gene activity scores (log(fragment counts per 10k + 1) on gene body +2kb promoter region) for H3K27me3, H3K27ac and H3K4me3 as well as RNA expression (log(transcript counts per 10k + 1)) over the telencephalic neuron differentiation trajectory from pluripotency.Pseudotime was binned and genes were Kmeans clustered based on average expression/activity of all marks and RNA in all bins (Supplementary Table 11 contains all clusters).Using this strictly refined analysis we could again identify a neuronal cluster that clearly showed priming with active chromatin showing detection rate of scATAC fragments and scRNA expression over pseudotime measured from the same cells of a primary human developing embryonic brain at gw 19.
Neuronal genes (NEUROD2, STMN2, BCL11B) show priming of chromatin.i, Boxplot (median +/-Q1/Q3) showing the shift in pseudotime between chromatin accessibility and RNA expression in the fetal cortex on the same set of neuronal genes as in b (n=357 genes).More details on the analysis can be found in the method section on Detection of inflection point for neuronal genes.

Fig. 1
Quality metrics of the scCUT&Tag dataset a, Histograms showing fragment length distribution for all samples (different timepoints including replicates) and modalities (H3K27ac, H3K4me3, H3K27me3), showing a clear nucleosome size pattern b, Violin plots showing the distribution of detected fragments per cell for each time point and chromatin modality (H3K27ac: n=33533, H3K27me3: n=34357, H3K2me3: n=42057, 8 independent timepoints each).The median value of each boxplot is indicated at the top (median +/-Q1/Q3).c, Table comparing median read per cell as a measure of data quality between different scCUT&Tag publications and the presented study (Zenk & Fleck et al. 2023).d, Histograms comparing the distribution of peak fragments per cell for each histone modification in all cell lines in the presented dataset.The red line demarcates the median.Supplementary Fig. 2 Astrocyte regionalization signatures in human brain organoids a-c, UMAP embedding of NPCs and astrocytes from the 30-240 day time points colored by cell state (a), brain region (b), and time point (c).d, UMAP embedding colored by gene expression (log(transcript counts per 10k + 1)) of general astrocyte marker genes.e, UMAP embedding colored by gene expression of general progenitor marker genes.f, UMAP embedding colored by gene expression of general outer radial glia marker genes.
modifications.c, Jitter plot showing even UMI count and fragment counts versus pseudotime for all modalities.d, Barplot showing fragment detection rate in pseudotime bins at selected example genes for all modalities.The line indicates smoothing with generalized additive models.Supplementary Fig. 4 Global analysis of neuronal genes in organoids and pseudotime reconstruction of neurogenesis in the primary developing cortex reveal widespread priming.a, Transcription factor motif enrichment in peaks on primed neuronal genes (gene body +2kb promoter region, see Methods Transcription factor motif enrichment)(cluster 6, Fig. 4a) plotted against their expression fold change in cortical neurons.Colored dots indicate a significant motif enrichment by Fisher exact test FDR<0.01.b, Ridge plots showing the distribution of pseudotime lag between active histone modifications (H3K27ac and H3K4me3) and RNA expression.Positive values indicate histone modifications preceding RNA, negative values the opposite.(see Methods Detection of inflection point for neuronal genes).Neuronal genes were selected based on having maximal expression in neurons over NPCs and detection of H3K27ac in more than 5% of cells in any neuronal high-resolution cluster.Primed genes show an enrichment of GO-terms related to neuron projection and axonogenesis, while co-transcriptionally marked genes are enriched for terms related to synapse organization.c, Boxplot (median +/-Q1/Q3) showing the shift in pseudotime between the establishment of active histone modifications, chromatin accessibility and RNA expression on the same set of neuronal genes as b (n=357 genes).d, Schematic of the multiome experiment in the primary developing brain.e, UMAP embedding of the cortical neuron trajectory from the gw 19 primary developing brain multiome data colored by cell state.f, UMAP embedding colored by pseudotime.g, UMAP embedding as is (b) colored by expression (log(normalized transcript counts +1)) of different neuronal and progenitor marker genes (top) and gene activity based on chromatin accessibility (log(fragment counts per 10k + 1) on gene body +2kb promoter region) (bottom).h, Barplots with smoothed lines