Inferring and perturbing cell fate regulomes in human brain organoids

Abstract

Self-organizing neural organoids grown from pluripotent stem cells^1-3 combined with single-cell genomic technologies provide opportunities to examine gene regulatory networks underlying human brain development. Here we acquire single-cell transcriptome and accessible chromatin data over a dense time course in human organoids covering neuroepithelial formation, patterning, brain regionalization and neurogenesis, and identify temporally dynamic and brain-region-specific regulatory regions. We developed Pando-a flexible framework that incorporates multi-omic data and predictions of transcription-factor-binding sites to infer a global gene regulatory network describing organoid development. We use pooled genetic perturbation with single-cell transcriptome readout to assess transcription factor requirement for cell fate and state regulation in organoids. We find that certain factors regulate the abundance of cell fates, whereas other factors affect neuronal cell states after differentiation. We show that the transcription factor GLI3 is required for cortical fate establishment in humans, recapitulating previous research performed in mammalian model systems. We measure transcriptome and chromatin accessibility in normal or GLI3-perturbed cells and identify two distinct GLI3 regulomes that are central to telencephalic fate decisions: one regulating dorsoventral patterning with HES4/5 as direct GLI3 targets, and one controlling ganglionic eminence diversification later in development. Together, we provide a framework for how human model systems and single-cell technologies can be leveraged to reconstruct human developmental biology.

Fig. 1. Multi-omic atlas of brain organoid development reveals developmental hierarchies and critical stages of fate decision.

a, Schematic of the experimental design and UMAP embedding of integrated multi-omic metacells. Organoids from three iPS cell lines and one ES cell line were dissociated for paired scRNA-seq and scATAC-seq at time points spanning 4 days to 2 months of development. The two modalities were integrated to form metacells with RNA and ATAC components. EB, embryoid body; IPs, intermediate progenitors; N.ect., neuroectoderm; N.epi., neuroepithelium; PSCs, pluripotent stem cells. b, Examples of loci with differential accessibility during organoid development from pluripotency. c, Schematic of the branch-inference strategy. High-resolution clusters were assigned to branches on the basis of terminal fate transition probabilities calculated based on RNA velocity. d, Branch visualization in a force-directed layout. The circles represent high-resolution clusters of metacells coloured by assignment (neuroepithelium (grey); non-telencephalon progenitors (teal); telencephalon progenitors (plum); dorsal telencephalon (orange); ventral telencephalon (purple)). e, Graph representation of regional branches coloured by mean expression (log[transcript counts per 10,000 + 1]) (top) and gene activity (log[transcript counts per 10,000 + 1]) (bottom) of marker genes. The range of values is indicated for each plot. Norm., normalized. f, Stage- and branch-specific gene expression and motif enrichment z-score (Methods). Values are minimum–maximum (min–max) scaled across rows. N.t., non-telencephalon; t., telencephalon.

Fig. 2. Pando leverages multimodal measurements to infer a multiphasic GRN underlying human brain organoid development.

a, Schematic of the Pando GRN-inference framework. Candidate regions are identified through intersection of accessible peaks with CREs or conserved elements. Predicted TFs are selected for each candidate region through binding-motif matching. The relationship between TF–binding-site pairs and the expression of target genes is then fitted with a regression model. E, expression; A, accessibility; g1, target gene 1; tf1,2, transcription factors; p1,4, peaks; GLM, generalized linear model; reg., regularized. b, Signal tracks showing normalized accessibility at the transcription start site of EMX1 in the different branches and inferred regulatory regions for various transcription factors. The line colour represents the sign of the interaction and the box colour (greyscale) represents the false-discovery rate (FDR) of the most significant interaction for this region. c, UMAP embedding of the inferred TF network based on co-expression and inferred interaction strength between TFs. Colour and size represent the expression-weighted pseudotime and PageRank centrality of each TF, respectively. d, UMAP embedding shaded by module features. e, Target specificity for branch-specific TFs. f, UMAP embedding of branch-specific TF networks highlighting TFs with branch-specific targets and interactions with branch-specific accessibility. g, Groups of TFs with differential activity between the dorsal (red) and ventral (purple) telencephalon branch. TF activity is indicated by a coloured dot for each branch, connected by a line, and was calculated by multiplying the mean regulatory coefficient (coef.) with the average expression (expr.) in the branch. The sign of the activity indicates whether the regulation is mainly activating (+) or repressing (−).

Fig. 3. TF perturbations in mosaic organoids reveal critical regulators of neurodevelopmental fate decisions.

a, Schematic of the single-cell TF perturbation experiment using the CRISPR droplet sequencing (CROP-seq) method. b, The minimum–maximum-scaled average expression (log[transcript counts per 10,000 + 1]) of targeted genes in NPCs, IPs and neurons of the primary and organoid cortex. c, The proportion of cells with each perturbation for each experiment. d, UMAP embedding with cells coloured by detected gRNA (left) and branch assignment (right). e, Regional enrichment of gRNAs. The sidebar shows the number of gRNAs that were consistent and the circles represent consistent effects between experiments and statistically significant (FDR < 0.01) effects on composition. The arrows indicate the predominant observed effect. f, UMAP embedding coloured by consistent gRNAs for selected genes that had a strong effect on fate regulation. g, The Spearman correlation of HES1-target (top, n = 18 genes) and GLI3-target (bottom, n = 42 genes) genes to transition probabilities into the dorsal branch. The GRN was subsetted to retain connections that are accessible at the branchpoint (>5% detection rate). The centre line represents the median, the box limits show the 25–75% interquartile range and the whiskers indicate 1.5× the interquartile range. h, Schematic of the GLI3 loss-of-function experiment using an inducible CRISPR–Cas9 nickase system. i, UMAP embedding of scRNA-seq data from 6-week-old WT and GLI3-KO brain organoids showing the trajectories from NPCs to neurons coloured by different clusters assigned to regional branches. The inset is coloured by genetic condition. j, Stacked bar plots showing the distribution of cluster (colour) assignment per organoid for each condition. k, Differential expression (DE) in ventral telencephalic neurons for GLI3-KO data and CROP-seq data containing a GLI3 gRNA. The x and y axes indicate the coefficients of the linear model. Colours indicate significance (FDR < 10⁻⁴) in CROP-seq, the KO cell line or both.

Fig. 4. Single-cell multiome view of GLI3 loss of function reveals distinct regulomes and effectors of dorsoventral telencephalon specification.

a, Schematic of the experiment measuring the transcriptome and chromatin accessibility in the same cell at 3 weeks of brain organoid development. b, UMAP embedding coloured by cluster and labelled by projected cell fate. Inset: UMAP coloured by genetic state. c, The number of DEGs of control (WT) versus GLI3-KO cells in the different clusters. d, Differential expression in telencephalic progenitors (clusters 0 and 2) after GLI3 KO. e, DEGs after GLI3 KO for early telencephalic progenitors (week 3), ventral telencephalic progenitors (week 6) and neurons (week 6), and differential accessibility after GLI3 KO in early telencephalic progenitors (week 3). Genes are coloured according to the associated signalling pathway (if applicable) and molecular function. f,g, GRN subgraph for early telencephalic (f) and ventral telencephalon (g) progenitors, showing first- and second-order GLI3 targets. The circles represent genes for which all TFs are labelled. The edges are coloured on the basis of TF regulatory interaction. h, The GLI3-binding score (the sum of CUT&Tag signal intensity for the gene body + 2 kb) in WT organoids versus log-transformed fold change in differential expression in early telencephalic progenitors (week 3). Genes with differentially accessible (DA) CREs are coloured black. Signal tracks of GLI3 binding matched with differential accessibility peaks of HES4 and HES5 in early telencephalic progenitors. i, The z-scored mean correlation between module gene expression and branch probabilities (branch activation score) for differentially expressed TFs. j, The log-transformed fold change of genes after treatment with SHH versus GLI3 KO. GO terms are shown for common DEGs, SHH-treatment-specific and GLI3-specific DEGs. k, Schematic summarizing the results from the GLI3 and SHH perturbations.

Extended Data Fig. 1. Supplemental analysis of brain organoid developmental multiome data.

a, Phase contrast (until day 15) and bright field (day 31–60) showing examples of different stages of organoid development for four different stem cell lines. Images are representative for 96 organoids per line. Scale bar is 200 µm. b, Schematic of the experimental design and data integration strategy. c, Histogram of scRNA-seq and scATAC-seq quality control metrics. d, Histograms showing assignment log likelihoods for demultiplexing based on single nucleotide variants. e, Bar plot of number of cells for each time point (top) and stacked barplot showing proportion of cell lines (bottom) at different time points. f, Distribution of iPS cell (iPSC) lines on the UMAP embedding. g, Bar plots showing number matched and unmatched cells during MCMF bipartite matching. h, Histogram showing the number of cells per metacell for each cell line. i, Box plots showing correlation between gene expression and gene activity metrics for two multiome experiments and the integrated metacells (n = 477 genes). j, Box plots showing correlation split by stage (n = 3527 genes). Genes >95% confidence intervals of correlation to permuted background are coloured in yellow. Box center represents the median, boxes indicate 25%–75% interquantile range and whiskers 1.5 * interquantile range. k, Immunohistochemical staining for progenitor cells (SOX2, orange and GLI3, purple) and neurons (TUJ1, green) for 2 month old organoids of four cell lines. DAPI is shown in cyan. Scale bar: 200 μm. l, UMAP embedding coloured by marker gene expression (log(transcript counts per 10k+1)). The range of values is indicated for each plot. m, Hierarchical clustering of pseudotemporal bins. Top bars show stage and proportion of time points per bin. Heatmap shows min-max scaled mean accessibility (tf-idf normalized fragment counts) of stage-specific peak clusters for each pseudotime bin. Representative GREAT enrichments are shown for each stage.

Extended Data Fig. 2. Heterogeneity analysis in different stages of organoid development.

a,d, UMAP embedding of a subset of the organoid trajectory surrounding neuroectoderm cells (a) and the branching window (d) coloured by time point, velocity pseudotime, cell line, branch prediction and lovain clusters. b,e, Heatmap showing mean min-max scaled expression (log(transcript counts per 10k + 1)) of cluster markers. c, UMAP embedding coloured by cluster identities, expression patterns of cluster markers. Volcano plot shows differentially expressed (DE) genes of cluster 5 relative to other clusters. f, UMAP embedding coloured by rank-transformed CellRank transition probability to non-telencephalon, ventral telencephalon and dorsal telencephalon and coloured by expression of selected transcription factors. g, Whole-mount HCR in situ hybridizations of day 18 organoids and UMAP embedding coloured by expression of targets. Stainings were performed on 2-3 organoids per cell line and representative images were shown. Scale bar: 100 μm. The range of expression values is indicated for each feature plot.

Extended Data Fig. 3. Signalling transcriptome and regulatory element landscape of the organoid neuroepithelium from 9 stem cell lines.

a, Schematic of the experimental setup. Multiome quantification was performed on organoids in the neuroepithelial stage (~3 weeks) from a total of 9 stem cell lines. The data was combined with the data from the same stage in the early time course. b, UMAP embedding coloured by cell line, louvain clusters and anterior-posterior axis (forebrain versus non-forebrain) classification score. c, Bar plots (top) showing fraction of cells per cell line in each cluster. Dotplot (left) showing min-max scaled expression (log(transcript counts per 10k + 1)) (colour) and proportion of expressing cells (dot size) for transcription factors (TFs) and genes from different signalling pathways in clusters of 3 week old organoid data set split by cell line. All genes are annotated as TF, receptor, ligand, or TF target and if applicable, coloured by the related signalling pathway. Dotplot (right) showing expression (colour) and proportion of expressing cells (dot size) for the same genes of Extended Data Fig. 3d in mouse developing brain organizer cells of different brain regions. d, Dot plot showing cluster-specific cis regulatory elements (CREs) linked to patterning genes split by different cell lines. Colour and size indicate peak accessibility (if-idf normalized fragment counts) and proportion of expressing cells, respectively.

Extended Data Fig. 4. Trajectory reconstruction in the multiomic developmental atlas.

a, Time course UMAP embedding coloured by neuron types. b, Time course UMAP embedding coloured by RNA velocity pseudotime. c, VoxHunt plots showing expression similarity of neuron subtypes in brain organoids to voxels in five example sections of the developing mouse brain (embryonic day 13.5), as well as the structural annotation of the sections. d, UMAP embedding coloured by ranked transition probabilities. e, Scatter plot showing mean transition probabilities as computed by CellRank versus velocity pseudotime. Each dot represents one high-resolution cluster. f, UMAP embedding of the integrated time course and graph embedding coloured by gene expression (log(transcript counts per 10k +1)) (top) and gene activity (log(fragment counts per 10k +1)) (bottom) for selected marker genes. g, UMAP and graph representation coloured by transcription factor motif enrichment z-score calculated with chromVAR for selected motifs. The range of values is indicated for each feature plot.

Extended Data Fig. 5. Gene regulatory network features of brain organoid development.

a, Numbers of chromatin access peaks and percentage of H3K27ac-marked peaks accessible at day 18–23 (>5% detection) intersecting with non-protein coding conserved regions (Cons.), candidate cis regulatory regions (CRE), or exons (left). b, Representative loci showing chromatin access (top) overlaying peak, CRE, conserved elements, and exon coordinates. c, Barplot showing the number of motifs used in GRN construction from two curated databases (JASPAR, CIS-BP), and motifs assigned through amino acid sequence similarity. d, Examples of 3 TFs with no motif annotation that were assigned motifs based on sequence similarity. e, Loci for two exemplary genes (FOXG1, WLS) showing average chromatin access signal tracks, accessible peaks, CREs, conserved elements, exons and H3K27ac CUT&Tag peaks. f, Scatter plot and histograms show explained variance (x) versus number of variables (y) of models for GRN construction. g, Violin plots show the distribution of peaks (left, n = 2535 target genes) and TFs per gene (middle, n = 2535 target genes), and number of genes per TF (right, n = 720 TFs). h, Representative loci showing average chromatin access signal tracks at different developmental branches overlaying inferred transcription factor binding sites within regulatory regions. i, UMAP representation of time course coloured by gene expression (log(transcript counts per 10k + 1)), gene module activity (module score calculated with Seurat) (rows), and regulatory module enrichment z-score (calculated with chromVAR) for representative TFs (columns). The range of values is indicated for each plot. j, Variation of module activity explained by branch, pseudotime, or branch and pseudotime (n = 720 TF modules). Box plot centre lines represent the median, boxes indicate 25%–75% interquantile range and whiskers 1.5 * interquantile range. k, Branch and pseudotime specific TF modules. Colours represent the branch with highest average module activity. TFs without experimentally validated motif are shown in grey.

Extended Data Fig. 6. Target selection and experimental details for the single-cell in organoid perturbation experiment.

a, Min-max scaled mean expression (log(transcript counts per 10k + 1)) of genes targeted in the single-cell genomic perturbation experiment in neuronal progenitors (NP), intermediate progenitors (IP) and neurons in the primary human and organoid developing cortex, as well as in iPS cells, the embryoid body (EB), ventral telencephalic NPCs, inhibitory neurons of the medial ganglionic eminence (MGE in.), lateral ganglionic eminence (LGE in.), non-telencephalic NPCs, diencephalic excitatory neurons (Dien. ex.) and inhibitory neurons (Dien. in.) and mesencephalic excitatory neurons (Mesen. ex.) and inhibitory neurons (Mesen. in.). b, UMAP embedding coloured by the expression of all targeted genes. The range of expression values is indicated for each feature plot. c, Exemplary Fluorescence-activated cell sorting plots of the sorting scheme used to isolate CROP-seq vector positive iPS cells. d, Phase contrast and CROP-seq vector positive (GFP) imaging during brain organoid development. Images are representative for 48 imaged organoids. Scale Bar is 500 µm.

Extended Data Fig. 7. Guide detection and cell type annotation in the single-cell perturbation experiment in organoids.

a, Barplot showing number of cells with detected guide RNA (gRNA) for each targeted gene and stacked barplot showing the distribution of the different gRNAs targeting the same gene. b, Histogram showing the distribution of read counts for gRNA UMIs after amplicon sequencing for one organoid. UMIs marked in red were selected for downstream analyses. c, Density histograms showing the distribution of inferred KO probabilities for gRNAs of 3 different target genes. d, Barplot showing cell number and proportion of gRNAs for all target genes. e, Barplot showing the number of guides detected in sequenced cells. f, UMAP embedding with cells coloured based on experiment. g, UMAP embedding coloured by annotated neuron subtypes. h, VoxHunt plots showing expression similarity of neuron subtypes in brain organoids to voxels in five example sections of the developing mouse brain (embryonic day 13.5), as well as the structural annotation of the sections (left). i, UMAP embedding coloured by expression (log(transcript counts per 10k + 1)) of non-telencephalic (top), ventral (middle) and dorsal (bottom) neuron markers. The range of expression values is indicated for each feature plot.

Extended Data Fig. 8. Composition and expression changes after CRISPR-Cas9 perturbations in mosaic brain organoids.

a, Hierarchical clustering of Louvain clusters based on the composition of gRNAs targeting different genes. Cell type and branch annotations are shown as side bars. Compositions of organoids and composition of cells with gRNAs targeting different genes are shown below as stacked bar plots. b, Lollipop plot showing the impact of each gRNA on cell type abundance in dorsal and ventral telencephalic neurons. c, Lollipop plots showing number of differentially expressed genes (DEG) for targeted genes in the dorsal and ventral telencephalic neurons. d, Differential gene expression analysis was performed to identify potential effects on cell state. Plot shows the effect of cell composition change and the number of differentially expressed genes (DEGs). P-values were derived using an F-test based ANCOVA. e, Scatter plot shows expression changes between neurons with E2F2 targeting gRNAs and other neurons in dorsal (x-axis) and ventral (y-axis) telencephalic neurons, with each dot representing one gene. Colours of dots represent the neuron types where differential expression is detected. Lines show the correlation of expression changes in the two neuron types, with DE genes in both types and DE genes in only one type shown separately. f, Examples of functional enrichment for E2F2 DEGs in dorsal and ventral neurons with DAVID. Grey bars show enriched terms of all E2F2 DEGs, and dark bars show enriched terms of DEGs with E2F2-specific effects.

Extended Data Fig. 9. Characterization of GLI3 knock-out organoids.

a, Quantification of editing frequency as determined by the percentage and number of reads showing unmodified and modified alleles for the control and both KO cell lines. b, Frequency of frameshift of coding sequence reads as a result of the modifications seen in both KO lines. c, Western blot showing expression of Gli3-repressor (83kDA) in the control cell line. Catenin beta-1 and Ponceau were used as loading control. For western blot source data, see Supplementary Fig. 1 d, Sequences of the coding strand of the different indels of the different KO lines. The reference sequence is corresponding with the control line. The position of the gRNAs with the protospacer adjacent motif (PAM)-sequence is depicted above and underneath the sequence. Reference protein sequence with the protein sequences of each KO line of the altered protein sequences caused by the frame-shift. e, Brightfield images of brain organoid development with control and both KO cell lines. Images are representative for 16 imaged organoids per line. Scale bar is 2 mm. f, UMAP embedding showing trajectories from neural progenitor cells (NPCs) to neurons coloured by different clusters assigned to branches (dorsal, ventral, and non-telencephalon), with inset coloured by genetic condition and feature plots coloured by expression (log(transcript counts per 10k + 1)) of cell type markers. g, UMAP embedding of ventral telencephalic GLI3 KO neurons showing medial ganglionic eminence (MGE) and lateral/caudal ganglionic eminence (LGE/CGE) neuronal populations (top). Feature plots show selected marker gene expression on the UMAP embedding. The range of expression values is indicated for each feature plot. h, Volcano plot showing differential expression analysis in LGE neurons for GLI3 WT versus KO cells. i, Schematic of observed effect of GLI3 loss of function on dorsoventral telencephalic fate decisions.

Extended Data Fig. 10. GLI3 KO induced changes in telencephalic progenitors in brain organoids.

a, Heatmap showing min-max scaled expression (log(transcript counts per 10k +1)) of marker genes for unbiased Louvain clusters. b, UMAP coloured by the expression of selected marker genes. The range of values is indicated for each plot. c, UMAP embedding coloured by branch labels predicted by label transfer from the organoid time course. d, Heatmap showing DE associated with signalling pathways. e, Scatter plot showing DE in neural progenitor cells (NPCs) upon HES1 perturbation in the mosaic perturbation experiment. f, Signal tracks showing differentially accessible (DA) peaks in cluster 0 and 2. g, GREAT enrichment analysis of DA peaks in cluster 0 and 2, with box area proportional to FDR. Representative genes are shown. h, Enrichment of DE genes in the neighbourhood of GLI3 in the GRN. i, Accuracy of GRN predicted directionality of GLI3 effect at different false discovery rate (FDR) thresholds. j, Barplot showing the number of all DE genes in the GRN (DEG in GRN), all DE genes reachable from GLI3 in the graph (Reachable), DE genes where the GRN was consistent with the DE result (Overall) and DE genes for which all subpaths from GLI3 were consistent with the DE result (Full path). k, Barplot showing the fraction of DE genes directly and indirectly regulated by GLI3. l, Boxplot showing the Spearman correlation of directly (n = 39) and indirectly (n = 126) regulated DE genes with transition probabilities ventral and dorsal branched. The centre line represents the median, boxes indicate 25%–75% interquantile range and whiskers 1.5 * interquantile range. m, Barplot showing the enrichment of gene sets (HES1/5 target genes, NOTCH components) among telencephalic DE genes. n, UMAP embedding showing annotation of multiome SHH experiment. o, UMAP coloured by the expression of selected marker genes.