Molecular and cellular dynamics of the developing human neocortex

Abstract

The development of the human neocortex is highly dynamic, involving complex cellular trajectories controlled by gene regulation¹. Here we collected paired single-nucleus chromatin accessibility and transcriptome data from 38 human neocortical samples encompassing both the prefrontal cortex and the 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 catalogue cell-type-specific, age-specific 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 neurogenesis-to-gliogenesis transition. We identified a tripotential intermediate progenitor subtype-tripotential intermediate progenitor cells (Tri-IPCs)-that is responsible for the local production of GABAergic neurons, oligodendrocyte precursor cells and astrocytes. Notably, 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 genome-wide association study data, we created a disease-risk map highlighting enriched risk associated with autism spectrum disorder in second-trimester intratelencephalic neurons. Our study sheds light on the molecular and cellular dynamics of the developing human neocortex.

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

a, Description of the samples used in this study. The diagram was created using BioRender. 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, The proportion of individual cell types across developmental stages and cortical regions. Bars are colour coded by cell types, the legend for which is shown in b. e, The expression of signature TFs in individual cell types (left). Middle, aggregated chromatin accessibility profiles at the promoter of signature TFs across cell types. The blue arrows represent each TF’s transcriptional start site and direction. Right, 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 colour coded by types or the niches to which they belong. b, The proportion of different cell types in individual niches. Niche numbers correspond to the legend in a. c, The neighbourhood enrichment scores of the PFC sample at infancy. The row and column annotations are colour coded by cell types as defined a. d, The percentage of significant intercellular communication determined by NCEM identified across all datasets. The row and column annotations are colour coded by cell types, as defined in a. e, The direction of cellular interactions mediated by neuregulin signalling (left). Right, the communication probability of example ligand–receptor pairs in the neuregulin signalling pathway from EN-IT-immature to other cell types. Empty space means that the communication probability is zero. P values were calculated using one-sided permutation tests. f, The direction of cellular interactions mediated by somatostatin signalling (left). Right, the communication probability of example ligand–receptor pairs in the somatostatin signalling pathway from IN-MGE-SST to other cell types. Empty space means that the communication probability is zero. P values were calculated using one-sided permutation tests.

Fig. 3. GRNs that establish cell identities.

a, The minimum–maximum-normalized TF expression levels, region-based AUC scores and gene-based AUC scores of selective eRegulons across cell types. b, GRNs 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 coloured. The size and transparency of the TF nodes represent their gene expression levels in each cell type. ACC, accessibility; GEX, gene expression; R2G, region to gene. c, UMAP plots of cells belonging to EN lineages, showing the nine trajectories. Cells are colour coded by types, regions, age groups or pseudotime. d, Standardized gene-based AUC scores of six eRegulon modules along the trajectories of EN lineages. eRegulons are colour 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 indicates adjusted P > 0.05. Statistical analysis was performed using one-sided hypergeometric tests; nominal P values were adjusted using the Benjamini–Hochberg method. f, BPs during EN differentiation. g, Trajectories of four IT neuron lineages. h, Differentially expressed genes between V1-specific and common EN-L4-IT neurons. Statistical analysis was performed using likelihood ratio tests; nominal P values were adjusted using the Benjamini–Hochberg method. i, Representative eRegulons (activators) involved in trajectory determination at BPs. j, UMAP representation of representative eRegulons involved in trajectory determination at BPs.

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

a, The expression patterns of surface proteins used for progenitor isolation. b, Schematic of the sorting strategy for isolation of progenitor subtypes (left). Right, phase-contrast images of progenitor subtypes after 5 days in culture. iSVZ, inner SVZ; oSVZ, outer SVZ; CP and SP, cortical plate and subplate. The diagram was created using BioRender. Scale bar, 50 μm. c, The 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 and 22 clones across three independent experiments. Scale bar, 100 μm. e, Immunostaining of descendants of Tri-IPCs 12 weeks after transplantation into mouse cortex, demonstrating the presence of astrocytes (GFAP⁺), OPCs or oligodendrocytes (SOX10⁺), and INs (GABA⁺). n = 2 injections. HNA, human nuclear antigen. The diagram was created using BioRender. Scale bar, 10 μm. f, SingleCellNet-predicted identities of INs and astrocytes derived from Tri-IPCs. g, Graphical summary of cell lineage relationships in late second-trimester human neocortex. The diagram was created using BioRender. h, UMAP plots of malignant GBM cells colour coded by SingleCellNet-predicted cell types. i, UMAP plots of malignant GBM cells colour coded by their main cellular states. j, The proportion of predicted cell types across different cellular states in malignant GBM cells. The legend is shown in h.

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

a, Standardized per-cell SCAVENGE TRS for four cognitive functions. b, Standardized per-cell SCAVENGE TRS for five brain disorders, including ASD, MDD, BPD, ADHD and SCZ. c, The proportion of cells with enriched trait relevance across cell types. Tiles with significant TRS enrichment (two-sided hypergeometric test, *FDR-adjusted P < 0.01 and odds ratio > 1.4) are annotated by their odds ratios. d, Standardized SCAVENGE TRS of four brain disorders plotted along the IT neuron lineage pseudotime. The best-fitted smoothed lines indicate the average TRS and the 95% confidence interval in each pseudotime bin. e, Standardized gene-based AUC scores for top ten disease-relevant eRegulons (activators) ranked by Spearman’s ρ along the IT neuron lineage pseudotime. eRegulons with SFARI ASD-associated genes as core TFs are highlighted in red. For the box plots in a,b, the centre line shows the median, the box limits show the 25th and 75th percentiles, and the whiskers show the standard error. For a,b, statistical analysis was performed using two-sided hypergeometric tests; *FDR-adjusted P < 0.01 and odds ratio > 1.4.

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

a, UMAP plots showing the distribution of cell subclasses in the snMultiome data prior to data filtering. b, UMAP plots showing the distribution of age groups in the snMultiome 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 processes (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 snMultiome data.

UMAP plots of all cells showing the expression levels of cell-type-specific marker genes. The coloured circles and numbers pinpoint specific cell types where the gene is highly 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.

a–f, 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 neighbourhood enrichment z scores of each MERFISH sample. The row and column annotations are colour 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 neighbourhood 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) signalling 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 signalling 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 and SSTR2. c, Scatter plots illustrating the Pearson correlation of log fold changes in individual genes between Ostreotide and L-054,264 treatments. d, 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 minimum-maximum normalized TF expression levels, region-based AUC scores, and gene-based AUC scores of eRegulons across cell types. d, Heatmap-dotplots showing the minimum-maximum 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 RG-vRG, EN-L4-IT, IN-MGE-PV, astrocyte-protoplasmics and OPCs. TF nodes and their links to enhancers are individually coloured. The size and the transparency of the TF nodes represent their gene expression levels in each cell type. GEX, gene expression; ACC, accessibility; R2G, region-to-gene. c–f, 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 coloured by their corresponding eRegulons. Region-to-gene links are shown as arcs and colour 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 colour 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 coloured by pseudotime. g, UMAP plots of excitatory neuron lineages coloured 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 colour 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 colour 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 colour 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 (DIV1). e, Quantification of six cell types after sorting on DIV1 (n = 5, 5, 5 cultures), 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 (DIV14). j, Quantification of six cell types after sorting on DIV14 (n = 5, 5, 5 cultures), 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 colour 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. Created in BioRender. Wang, L. (2024) BioRender.com/t85o210. 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 cultures), 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. Created in BioRender. Wang, L. (2024) BioRender.com/f96w125. 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 colour 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 colour 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 colour 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 colour 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 colour 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 colour 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 centre: 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 centre: 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 odds ratios. e, Boxplots showing standardized SCAVENGE TRS for nine cognitive and disease traits across age groups. Boxplot centre: 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 odds ratios.