Molecular and cellular dynamics of the developing human neocortex at single-cell resolution

Abstract

The development of the human neocortex is a highly dynamic process and involves complex cellular trajectories controlled by cell-type-specific gene regulation¹. Here, we collected paired single-nucleus chromatin accessibility and transcriptome data from 38 human neocortical samples encompassing both the prefrontal cortex and primary visual cortex. These samples span five main developmental stages, ranging from the first trimester to adolescence. In parallel, we performed spatial transcriptomic analysis on a subset of the samples to illustrate spatial organization and intercellular communication. This atlas enables us to catalog cell type-, age-, and area-specific gene regulatory networks underlying neural differentiation. Moreover, combining single-cell profiling, progenitor purification, and lineage-tracing experiments, we have untangled the complex lineage relationships among progenitor subtypes during the transition from neurogenesis to gliogenesis in the human neocortex. We identified a tripotential intermediate progenitor subtype, termed Tri-IPC, responsible for the local production of GABAergic neurons, oligodendrocyte precursor cells, and astrocytes. Remarkably, most glioblastoma cells resemble Tri-IPCs at the transcriptomic level, suggesting that cancer cells hijack developmental processes to enhance growth and heterogeneity. Furthermore, by integrating our atlas data with large-scale GWAS data, we created a disease-risk map highlighting enriched ASD risk in second-trimester intratelencephalic projection neurons. Our study sheds light on the gene regulatory landscape and cellular dynamics of the developing human neocortex.

Extended Data Fig. 1 |. Filtering of the snMultiome data.

a, UMAP plots showing the distribution of cell subclasses in the single-nucleus multiome data prior to data filtering. b, UMAP plots showing the distribution of age groups in the single-nucleus multiome data prior to data filtering. c, UMAP plots showing the distribution of cells removed during data filtering. d, UMAP plots showing the expression levels of genes identified in the striatum (ISL1 and SIX3), diencephalon (OTX2 and GBX2), neuronal dendrites (NRGN), and oligodendrocyte processes (MBP). e, Classes, subclasses, and types identified from the snMultiome data.

Extended Data Fig. 2 |. Quality control of the snMultiome data.

a, Violin plots, box plots, barplots, and UMAP plots of several quality control metrics for evaluating the quality of individual samples in the snMultiome data, including numbers of unique molecular identifiers (# UMIs), numbers of identified genes (# genes), number of fragments in ATAC peaks, transcription start site (TSS) enrichment scores, and proportion of individual cell types in each sample. The legend for cell types can be found in panel b. b UMAP plots of cells from individual snMultiome datasets separated by age groups. c, UMAP plots generated based on RNA or ATAC data only. The legend can be found in panel a.

Extended Data Fig. 3 |. Expression patterns of marker genes in the single-nucleus multiome data.

UMAP plots of all cells showing the expression levels of cell-type-specific marker genes. The colored circles and numbers pinpoint specific cell types where the gene is expressed. The legend for these numbers can be found in Fig. 1b.

Extended Data Fig. 4 |. Quality control and annotation of MERFISH data.

a, Violin plots, box plots, barplots, and UMAP plots of several metadata of MERFISH samples, including numbers of detected transcripts (# transcript), numbers of identified genes (# genes), age groups, regions, cell types, and niches. b, PCA plots based on cell type proportions for individual snMultiome and MERFISH samples in three different age groups. c, UMAP plots of all cells in the MERFISH dataset showing the expression levels of cell-type-specific marker genes.

Extended Data Fig. 5 |. Spatial distribution of cell types in individual MERFISH samples.

Extended Data Fig. 6 |. Difference in distribution of MGE- and CGE-derived interneurons in the second-trimester neocortex.

a, Immunostaining of MGE-derived (LGX6⁺) and CGE-derived (NR2F2⁺) interneurons in the cortex of a gestational week (GW) 24 sample. MZ, marginal zone; CP, cortical plate; SP/IZ, subplate/intermediate zone; oSVZ, outer subventricular zone; iSVZ, inner subventricular zone; VZ, ventricular zone. b, Odds ratios of the number of CGE-derived interneurons in the MZ versus ventricular/subventricular zones relative to the number of MGE-derived interneurons. Data are presented as mean values with 95% confidence intervals. P values were obtained from two-sided Fisher’s exact test; ***P < 0.001, ****P < 0.0001.

Extended Data Fig. 7 |. Intercellular communication between cell types in developing human cortex.

a, Heatmaps showing neighborhood enrichment z scores of each MERFISH sample. The row and column annotations are color-coded by cell types, the legend of which can be found in Fig. 2a. When a particular cell type is not present in the dataset, the neighborhood enrichment z scores were arbitrarily set to −50. b, Circular maps showing significant intercellular communication determined by NCEM in each MERFISH sample. c, Heatmaps showing the relative strength of outgoing (left) and incoming (right) signaling pathways in individual cell types. The bar graphs on the top and right side of the heatmaps are the sum of communication probability (interaction strength) for each cell type and signaling pathway, respectively.

Extended Data Fig. 8 |. Effects of somatostatin on the transcriptome of excitatory neurons in the second-trimester human cortex.

a, Violin plots, box plots, barplots, and UMAP plots of several metadata of scRNA-seq datasets from organotypic human brain slice cultures treated with and without somatostatin receptor agonists, including numbers of unique molecular identifiers (# UMIs), numbers of identified genes (# genes), ages, treatments, and cell types. b, UMAP plots of cells in the scRNA-seq dataset showing the expression levels of cell-type-specific marker genes. c, Gene set enrichment analysis (GSEA) highlighting the effects of L-054,264 and Ostreotide on different types of excitatory neurons. Significant terms, defined by Benjamini–Hochberg adjusted P values < 0.05, were outlined by a red circle. Abs(NES), absolute values of normalized enrichment scores.

Extended Data Fig. 9 |. SCENIC+ identifies cell-type-specific eRegulons.

a, Enrichment of eRegulon-predicted TF binding sites in ChIP-seq peaks from the human dorsolateral prefrontal cortex. P values were obtained from the two-sided Fisher’s exact test and adjusted using the Benjamini and Hochberg method. b, Overlap between eRegulon-predicted enhancer-promoter interactions and PLAC-seq loops from the developing human cortex. The P value was obtained from the two-sided Fisher’s exact test. c, Heatmaps showing the min-max normalized TF expression levels, region-based AUC scores, and gene-based AUC scores of activator eRegulons across cell types. d, Heatmap-dotplots showing the min-max normalized TF expression levels, region-based AUC scores, and gene-based AUC scores of selective eRegulons across age groups in all cells, Glutamatergic neurons, and GABAergic neurons.

Extended Data Fig. 10 |. Cell-type-specific gene regulatory networks in the developing cortex.

a, A heatmap showing Jaccard similarity matrix of target regions of cell-type-specific eRegulons listed in Fig. 3a. b, Gene regulatory networks of selective eRegulons in Astrocyte-Protoplasmics and OPCs. TF nodes and their links to enhancers are individually colored. The size and the transparency of the TF nodes represent their gene expression levels in each cell type. c, Coverage plots showing aggregated ATAC profiles across cell types on four genomic loci—SOX6, PDGFRA, HOPX, and GAD2. Identified candidate cis-regulatory elements (cCREs) are colored by their corresponding eRegulons. Region-to-gene links are shown as arcs and color-scaled based on region–gene importance scores obtained from SCENIC+ analysis.

Extended Data Fig. 11 |. Differentiation trajectories of excitatory neuron lineages.

a–e, UMAP plots of cells belonging to excitatory neuron lineages with clusters connected by a minimum spanning tree showing. The green node indicates the root node, and the red nodes indicate the ending nodes. Cells are color-coded by clusters (a), types (b), age groups (c), regions (d), or pseudotime (e). f, UMAP plots of each of the nine excitatory neuron lineages colored by pseudotime. g, UMAP plots of excitatory neuron lineages colored by the five pseudotime segments used for eRegulon activity analysis at bifurcation points. h, UMAP plots highlighting representative eRegulons involved in trajectory determination at bifurcation points.

Extended Data Fig. 12 |. Markers of V1-specific EN-L4-IT subtype.

a, UMAP plots of all EN-L4-IT color-coded by regions (left) and subtypes (right). b, UMAP plots showing the expression levels of representative differentially expressed genes between V1-specific and common EN-L4-IT neurons. c, In situ hybridization (ISH) of V1-biased (CUX1 and KCNIP1), and common-biased genes in EN-L4-IT neurons in adult human V1 and V2 areas. d, UMAP plots of EN-L4-IT subtype marker genes found in adult human V1.

Extended Data Fig. 13 |. Markers of human glial cells and their isolation strategies.

a, UMAP plots of cells belonging to glial lineages color-coded by age groups (left), regions (middle), and types (right). b, UMAP plots of cells belonging to glial lineages showing the expression levels of typical marker genes of individual cell types. c, UMAP plots of GW20 to GW23 cells belonging to glial lineages color-coded by age groups (left), regions (middle), and types (right). d, UMAP plots of GW20 to GW23 cells belonging to glial lineages showing the expression levels of typical marker genes of individual cell types. e, UMAP plots of GW20 to GW23 cells belonging to glial lineages showing the expression levels of surface markers used for glial progenitor isolation. f, Schematic of the sorting strategy for glial progenitors. VZ & iSVZ, ventricular zone and inner subventricular zone; oSVZ, outer subventricular zone.

Extended Data Fig. 14 |. Immunostaining characterization of human glial progenitor differentiation.

a–d, Immunostaining of isolated glial progenitors on days in vitro 1. e, Quantification of six cell types after sorting on days in vitro 1 (n = 5, 5, 5 samples), including RG or IPC-EN (TFAP2C⁺), IPC-Glia (OLIG2⁺EGFR⁺), OPC or oligodendrocyte (OLIG2⁺EGFR⁻), astrocyte (SPARCL1⁺), EN (NeuN⁺), and IPC-IN or IN (DLX5⁺). f–i, Immunostaining of progenies of glial progenitors on days in vitro 14. j, Quantification of six cell types after sorting on days in vitro 14 (n = 5, 5, 5 samples), including RG or IPC-EN (TFAP2C⁺), IPC-Glia (OLIG2⁺EGFR⁺), OPC or oligodendrocyte (OLIG2⁺EGFR⁻), astrocyte (SPARCL1⁺), EN (NeuN⁺), and IPC-IN or IN (DLX5⁺).

Extended Data Fig. 15 |. ScRNA-seq characterization of human glial progenitor differentiation.

a–d, UMAP plots of isolated glial progenitors and their progenies during in vitro differentiation based on single-cell RNA sequencing data color-coded by datasets (a), stages (b), seeding cell types (c), and types (d). e, UMAP plots of isolated glial progenitors and their progenies showing the expression levels of typical marker genes of individual cell types. f, A Sankey plot showing the mapping of glial progenitors and their progenies to the snMultiome atlas by SingleCellNet. g, UMAP plots of isolated glial progenitors and their progenies separated by seeding cell types and stages.

Extended Data Fig. 16 |. Lineage potential of human glial progenitors.

a, Schematic of the slice transplantation assay for glial progenitors. b–d, Immunostaining of progenies after progenitor transplantation to acute cortical slices on days in vitro 8. e, Quantification of progeny types after progenitor transplantation to acute cortical slices (n = 5, 5, 5, 5 samples), including IPC-EN (EOMES⁺), EN (NeuN⁺), Tri-IPC (OLIG2⁺EGFR⁺), astrocyte (SPARCL1+), OPC or oligodendrocyte (OLIG2⁺EGFR⁻), and IPC-IN or IN (DLX5⁺). f, Schematic of the in vivo transplantation assay for glial progenitors. g, Immunostaining of progenies after progenitor in vivo transplantation into mouse cortex (n = 2 injections). White arrows indicate HNA⁺GABA⁺ inhibitory neurons. HNA, human nuclear antigen; L2–3, layer 2–3; L6, layer 6; WM, white matter; V-SVZ, ventricular-subventricular zone; OB, olfactory bulb. h, Immunostaining of progenies after progenitor in vivo transplantation into mouse cortex (n = 2 injections). White arrows indicate HNA⁺SOX10⁺ OPCs or oligodendrocytes. Yellow arrows indicate HNA⁺GFAP⁺ astrocytes. HNA, human nuclear antigen; L2–3, layer 2–3; L6, layer 6; WM, white matter; V-SVZ, ventricular-subventricular zone; OB, olfactory bulb.

Extended Data Fig. 17 |. Mapping Tri-IPC progenies to reference data.

a, UMAP plot of a reference human ganglionic eminence dataset. Cells are color-coded by types. b, UMAP plots of human ganglionic eminence cells showing the expression levels of typical marker genes of individual cell types. c, UMAP plots of Tri-IPC-derived INs projected to the human ganglionic eminence dataset. Cells are color-coded by types and the legend can be found in panel d. d, Identities of Tri-IPC-derived INs mapped by Seurat label transfer. e, UMAP plot of mouse astrocytes from a reference developing mouse cortex dataset. Cells are color-coded by lineages and the legend can be found in panel h. f, UMAP plots of the reference mouse astrocytes showing the expression levels of typical marker genes of individual astrocyte lineages. g, UMAP plots of Tri-IPC-derived astrocytes projected to the reference mouse astrocytes. Cells are color-coded by lineages and the legend can be found in panel h. h, Identities of Tri-IPC-derived astrocytes mapped by Seurat label transfer. i, UMAP plot of human astrocytes at the infancy stage. Cells are color-coded by lineages and the legend can be found in panel l. j, UMAP plots of human astrocytes showing the expression levels of typical marker genes of individual astrocyte lineages. k, UMAP plots of Tri-IPC-derived astrocytes projected to the reference human astrocytes. Cells are color-coded by lineages and the legend can be found in panel l. l, Identities of Tri-IPC-derived astrocytes predicted by SingleCellNet (top) or mapped by Seurat label transfer (bottom). m, Proportion of each SingleCellNet-predicted cell type across GBM samples.

Extended Data Fig. 18 |. Neocortical cell association with human cognition and brain disorders.

a, UMAP plot showing the standardized per-cell SCAVENGE trait relevance score (TRS) for Alzheimer’s disease. b, Top, boxplots showing the standardized SCAVENGE TRS for Alzheimer’s disease across cell types. Boxplot center: median; hinges: the 25th and 75th percentiles; whiskers: standard error. Bottom, bar plots showing the proportion of the cells with enriched trait relevance for Alzheimer’s disease across cell types. Two-sided hypergeometry test; *FDR < 0.01 & odds ratio > 1.4. c, Boxplots showing standardized SCAVENGE TRS for nine cognitive and disease traits across regions. Boxplot center: median; hinges: the 25th and 75th percentiles; whiskers: standard error. Two-sided hypergeometry test; *FDR < 0.01 & odds ratio > 1.4. d, Heatmap showing the proportion of the cells with enriched trait relevance across regions. Tiles with significant TRS enrichment (two-sided hypergeometric test, *FDR < 0.01 & odds ratio > 1.4) are annotated by their odd ratios. e, Boxplots showing standardized SCAVENGE TRS for nine cognitive and disease traits across developmental stages. Boxplot center: median; hinges: the 25th and 75th percentiles; whiskers: standard error. Two-sided hypergeometry test; *FDR < 0.01 & odds ratio > 1.4. f, Heatmap showing the proportion of the cells with enriched trait relevance across developmental stages. Tiles with significant TRS enrichment (two-sided hypergeometric test, *FDR < 0.01 & odds ratio > 1.4) are annotated by their odd ratios.

Fig. 1 |. A multi-omic survey of the developing human neocortex.

a, Description of samples used in this study. b, UMAP plots of the snMultiome data showing the distribution of 33 cell types. c, UMAP plots showing the distribution of age groups (left) and regions (right). d, Proportion of individual cell types across developmental stages and cortical regions. Bars are color-coded by cell types, the legend of which can be found in panel a. e, Left, a dotplot of the signature transcriptional factors (TFs) in individual cell types. Middle, aggregated chromatin accessibility profiles on the promoter of signature TFs across cell types. The blue arrow represents each TF’s transcriptional starting site and gene body. Right, heatmap of normalized chromVar motif activity of signature TFs across cell types.

Fig. 2 |. Cell-cell communication in the developing human neocortex.

a, Spatial transcriptomic analysis of six neocortical samples. Cells are color-coded by types or the niches to which they belong. b, Proportion of different cell types in individual niches. Niche numbers correspond to the legend in panel a. c, Heatmap showing neighborhood enrichment scores of the PFC sample at infancy. The row and column annotations are color-coded by cell types, the legend of which can be found in panel a. d, Heatmap showing the percentage of significant intercellular communication determined by NCEM identified across all datasets. The row and column annotations are color-coded by cell types, the legend of which can be found in panel a. e, Left, a circular plot showing the direction of cellular interactions mediated by neuregulin signaling. Right, a dotplot showing communication probability of example ligand-receptor pairs in the neuregulin signaling pathway from EN-IT-Immature to other cell types. Empty space means the communication probability is zero. P-values were calculated by one-sided permutation test. f, Left, a circular plot showing the direction of cellular interactions mediated by somatostatin signaling. Right, a dotplot showing communication probability of example ligand-receptor pairs in the somatostatin signaling pathway from IN-MGE-SST to other cell types. Empty space means the communication probability is zero. P-values were calculated by one-sided permutation test.

Fig. 3 |. Gene regulatory networks that establish cell identities.

a, A heatmap-dotplot showing the min-max normalized TF expression levels, region-based AUC scores, and gene-based AUC scores of selective eRegulons across cell types. b, Gene regulatory networks of selective eRegulons in three distinct cell types (RG-vRG, EN-L4-IT, and IN-MGE-PV). TF nodes and their links to enhancers are individually colored. The size and the transparency of the TF nodes represent their gene expression levels in each cell type. c, UMAP plots of cells belonging to excitatory neuron lineages showing the nine trajectories. Cells are color-coded by types, regions, age groups, or pseudotime. d, Standardized gene-based AUC scores of six eRegulon modules along the trajectories of excitatory neuron lineages. eRegulons are color-coded by neuronal trajectories. Thick, non-transparent lines represent the average AUC scores of each module in each trajectory. e, Gene ontology enrichment analysis for target genes of individual eRegulon modules. Empty space means adjusted P values > 0.05. One-sided hypergeometric test; nominal P values were adjusted by the Benjamini and Hochberg method. f, Bifurcation points during excitatory neuron differentiation. g, Trajectories of four intratelencephalic neuron lineages. h, Volcano plots highlighting differentially expressed genes between V1-specific and common EN-L4-IT neurons. Likelihood ratio test; nominal P values were adjusted by the Benjamini and Hochberg method. i, A dotplot highlighting representative eRegulons (activators) involved in trajectory determination at bifurcation points. j, UMAP plots highlighting representative eRegulons involved in trajectory determination at bifurcation points.

Fig. 4 |. Multipotent progenitors during transition from neurogenesis to gliogenesis.

a, Violin plots showing the expression patterns of surface proteins used for progenitor isolation. b, Left, schematic diagram showing the sorting strategy for isolation of progenitor subtypes. Right, phase-contrast images of progenitor subtypes after five days in culture. VZ & iSVZ, ventricular zone and inner subventricular zone; oSVZ, outer subventricular zone; CP & SP, cortical plate and subplate. c, Proportion of individual cell types across progenitor subtypes and differentiation stages during progenitor differentiation in vitro. d, Clonal analysis demonstrating multipotency of individual progenitor cells (n = 26, 29, 22 clones across three independent experiments). e, Immunostaining of progenies of Tri-IPCs 12 weeks after transplantation into mouse cortex, demonstrating presence of astrocytes (GFAP⁺), OPCs or oligodendrocytes (SOX10⁺), and IN (GABA⁺) (n = 2 injections). HNA, human nuclear antigen. f, SingleCellNet predicted identities of interneurons (INs) and astrocytes derived from Tri-IPCs. g, Graphical summary of cell lineage relationships in late second-trimester human neocortex. h, UMAP plots of malignant GBM cells color-coded by their main cellular states. i, UMAP plots of malignant GBM cells color-coded by SingleCellNet-predicted cell types. j, Proportion of predicted cell types across different cellular states in malignant GBM cells. The legend can be found in panel i.

Fig. 5 |. Cell type association with human cognition and brain disorders.

a, Standardized per-cell SCAVENGE trait relevance score (TRS) for four cognitive functions. Boxplot center: median; hinges: the 25th and 75th percentiles; whiskers: standard error. b, Standardized per-cell SCAVENGE TRS for five brain disorders, including autism spectrum disorder (ASD), major depressive disorder (MDD), bipolar disorder (BPD), attention-deficit/hyperactivity disorder (ADHD), and schizophrenia (SCZ). Boxplot center: median; hinges: the 25th and 75th percentiles; whiskers: standard error. Two-sided hypergeometric test, *FDR < 0.01 & odds ratio > 1.4. c, Heatmap showing the proportion of the cells with enriched trait relevance across cell types. Tiles with significant TRS enrichment (two-sided hypergeometric test, *FDR < 0.01 & odds ratio > 1.4) are annotated by their odd ratios. d, Standardized SCAVENGE TRS of four brain disorders plotted along the intratelencephalic (IT) neuron lineage pseudotime. The best-fitted smoothed lines indicate the average TRS and the 95% confidence interval in each pseudo-time bin. e, Heatmaps of standardized gene-based AUC scores for top ten disease-relevant eRegulons ranked by Spearman’s ρ along the IT neuron lineage psuedotime. eRegulons with SFARI ASD-associated genes as core TFs are highlighted in red.