Spatial Single-cell Analysis Decodes Cortical Layer and Area Specification

Abstract

The human cerebral cortex, pivotal for advanced cognitive functions, is composed of six distinct layers and dozens of functionally specialized areas^1,2. The layers and areas are distinguished both molecularly, by diverse neuronal and glial cell subtypes, and structurally, through intricate spatial organization^3,4. While single-cell transcriptomics studies have advanced molecular characterization of human cortical development, a critical gap exists due to the loss of spatial context during cell dissociation^5,6,7,8. Here, we utilized multiplexed error-robust fluorescence in situ hybridization (MERFISH)⁹, augmented with deep-learning-based cell segmentation, to examine the molecular, cellular, and cytoarchitectural development of human fetal cortex with spatially resolved single-cell resolution. Our extensive spatial atlas, encompassing 16 million single cells, spans eight cortical areas across four time points in the second and third trimesters. We uncovered an early establishment of the six-layer structure, identifiable in the laminar distribution of excitatory neuronal subtypes by mid-gestation, long before the emergence of cytoarchitectural layers. Notably, while anterior-posterior gradients of neuronal subtypes were generally observed in most cortical areas, a striking exception was the sharp molecular border between primary (V1) and secondary visual cortices (V2) at gestational week 20. Here we discovered an abrupt binary shift in neuronal subtype specification at the earliest stages, challenging the notion that continuous morphogen gradients dictate mid-gestation cortical arealization^6,10. Moreover, integrating single-nuclei RNA-sequencing and in situ whole transcriptomics revealed an early upregulation of synaptogenesis in V1-specific Layer 4 neurons, suggesting a role of synaptogenesis in this discrete border formation. Collectively, our findings underscore the crucial role of spatial relationships in determining the molecular specification of cortical layers and areas. This work not only provides a valuable resource for the field, but also establishes a spatially resolved single-cell analysis paradigm that paves the way for a comprehensive developmental atlas of the human brain.

Fig. 1:. A spatially resolved single-cell atlas of human fetal cortical development.

a, Schematics of sampling and workflow. PFC; prefrontal cortex; PMC, premotor cortex; M1, primary motor cortex; S1, primary somatosensory cortex; Par, parietal cortex; Occi, occipital cortex; Temp, temporal cortex; Cing, cingulate cortex; Hippo, hippocampus; V1, primary visual cortex; V2, secondary visual cortex. b, c, Uniform Manifold Approximation and Projection (UMAP) of single cells analyzed by MERFISH, colored by H1 cell classes (b), and by H2 cell types (c). RG, radial glia; IPC, intermediate progenitor cells; EN-Mig, migrating excitatory neurons; EN-IT, intratelencephalic excitatory neurons; EN-ET, extratelencephalic excitatory neurons; IN, inhibitory neurons; EC, endothelial cells; oRG, outer radial glia; tRC, truncated radial glia; vRG, ventricular radial glia; INP, inhibitory neuron progenitors; CGE, caudal ganglionic eminence; MGE, medial ganglionic eminence; OPC, oligodendrocyte precursor cells. d, Spatial maps of H2 cell types from major cortical areas across gestational week (GW) 15 to 34. MZ, marginal zone; CP, cortical plate; SP, subplate; IZ, intermediate zone; oSVZ, outer subventricular zone; iSVZ, inner subventricular zone; VZ, ventricular zone. Scale bars, 500 µm. e, Dot plot showing expression of marker genes for H2 cell types. f, Cell type correspondence heatmap shows the fraction of cells from MERFISH H2 clusters that associate to clusters from published mid-gestation scRNA-seq dataset. g, Schematics illustrating the annotated fan-shaped cortical area used for relative height (RH) and cortical depth (CD) calculation. h, Violin plots showing the laminar distribution of H2 cell types from the apical to basal surface quantified by the RH. The width of the violin for each cell type is normalized to the maximum value. i, Violin plots showing CD distribution of H3 EN subtypes within the CP. Dash lines represent the borders between cortical layers (L2-L6) calculated based on the CD distribution of layer-defining clusters. Clusters with fewer than 50 cells within the annotated region are represented by dots for individual cells instead of a violin.

Fig. 2:. Progressive formation and specification of cortical layers.

a, Spatial map (left) shows the cortical layers 2 to 6 are defined by the organization of selected EN subtypes in the prefrontal cortex (PFC) at GW34; and ridgeline plot (right) shows the cortical depth (CD) distribution of layer-defining EN subtypes. The height of the ridgeline represents cell density. Dash lines represent the border between layers, calculated based on the CD distribution of layer-defining EN subtypes. b, Ridgeline plots reveal further spatial complexity among EN-ETs, deep layer (Layers 5&6) EN-IT, and upper layer (Layers 2–4) EN-ITs in the PFC at GW34. c, Spatial maps and ridgeline plots show the six-layer structure can be visualized and quantitatively defined by the distribution of layer-defining EN subtypes at GW22, despite lack of cytoarchitectural differences between layers. d, e, Spatial maps and ridgeline plots showing the cortical layers in the PFC at GW15 and GW20. f, Histogram showing the relative proportion of each cortical layer at GW15 to 34 in the PFC, parietal cortex (Par), and occipital cortex (Occi). Average is taken for replicate experiments; error bars represent standard deviation when applicable. g, h, i, Spatial maps (left) and ridgeline plots (right) showing the distribution of EN-ETs (g), deep layer EN-ITs (h), and upper layer EN-ITs (i) in the PFC and V2 at GW22. While many clusters exhibited area-dependent abundance, their laminar localization is highly consistent between PFC and V2. j, Violin plots showing the laminar distribution of H3 inhibitory neuron (IN) subtypes within the CP of PFC at GW22 and GW34. Scale bars, 500 µm.

Fig. 3:. Laminar and areal dynamic gene expression underlies neuronal specification.

a, Table compares the laminar expression patterns of validated layer marker genes between the adult human PFC and mid-gestation PFC. Genes that show conserved laminar expression pattern are colored in green, otherwise in red. b, Violin plots show the laminar expression patterns of layer-dependent genes within the CP of the PFC and V2 at GW22. Width of the violins represent the cumulative normalized expression in the cells located at a cortical depth (CD). Genes that are enriched in different layers between PFC and V2 are highlighted in red. c, Summary expression heatmap and z-score spatial maps showing the spatial-temporal expression pattern of CBLN2. The summary heatmap shows different cortical layers and laminar structures from top to bottom as rows, and cortical areas and gestational ages as columns, organized from anterior to posterior and from young to old. Average is taken for replicate samples from the same area and GW. T, temporal cortex. Scale bar, 500 µm. d, Summary expression heatmap and z-score spatial maps for CYP26A1. e, Summary expression heatmaps and z-score spatial maps of genes identified with Layers 2&3 enrichment in the PFC. f, UMAP plot of EN subtype composition in the CP shows an anterior-posterior spectrum for GW20-22. g, Dot plots and histograms of fold change for top differentially expressed genes (DEGs) between pairs of anterior- and posterior-enriched EN subtypes for each cortical layer at GW22. Genes that appear repeatedly are highlighted in red (for anterior-enriched) or in blue (for posterior-enriched). h, Spatial graphs showing the laminar and areal distribution of pairs of anterior- and posterior-enriched EN subtypes. i, Dot plot showing the expression of top areally-enriched genes at GW20 and 22 in all post-migratory EN (EN-IT and EN-ET cell classes). j, Schematic summary of the combinations of 5 marker genes that enable the approximation of cortical layers of different cortical areas at GW20 to 22.

Fig. 4:. Diversification of EN subtypes over time show continued post-mitotic fate specification.

a, b, Sankey diagram showing the correspondence between EN-ET (a) and EN-IT (b) subclusters from different gestational time points. Clustering was performed separately for each gestational age and the number of subclusters was determined single cell significant hierarchical clustering (scSHC). Nodes represent EN subtypes from this alternative clustering strategy. To differentiate these clusters from those identified through integrated analysis, an apostrophe (‘) prefixes the names of scSHC clusters. Thickness of edges represents the fraction of cells showing with correspondence, and edges with <0.1 fraction were hidden. Nodes and edges are colored based on the most enriched layer identity of the EN subtype. c, Spatial graphs for an inferred lineage (bolded in a) of Layer 6/SP EN-ETs across GW15 to GW34 show a conserved stream from GW15 to 22 before dramatic diversification at GW34. Scale bars, 500 µm. d, Spatial graphs showing a GW22 cluster EN-ITs spamming both Layers 4 and 5 specifying into seven GW34 clusters with refined layer specificity in either Layer 4 or 5 (bolded in b). Scale bars, 500 µm. e, Dot plots showing the top genes that are down- or up-regulated with gestational age among EN subtypes from each layer-based group. The genes that are downregulated over time in some categories but upregulated in others are highlighted in red. f, Curve plots show the change of relative expression over time in different groups for selected genes. g, Heatmap showing the genes with high expression variance among EN subtypes within the same group and gestational age. h, i, Dot plots showing the expression of high variance genes among EN subtypes in EN-ET-SP/L6b group at GW34 and EN-ET-L5/6 group at GW22.

Fig. 5:. Sharp molecular border between V1 and V2 at GW20 reveals early V1 specification.

a, The primary (V1) and secondary visual cortices (V2) do not exhibit morphological differences at GW20-21. Schematics and Nissl staining are taken from the Allen Reference Atlas for GW21 human fetal brain. DAPI staining image co-captured with MERFISH; H&E staining image co-captured with Visium. b, c Spatial maps for selected EN subtypes show distinct border between V1 and V2 marked by sharp transition at GW20. Scale bars 500 µm. d, Ridgeline plots showing the horizontal cell density profile for V1- and V2-enriched EN subtypes in the CP near the border. e, Ridgeline plots showing the laminar distribution for V1- and V2-enriched EN subtypes. Dash lines represent the borders between cortical layers. f, Histograms for areal distribution show V2-enriched subtypes distribute broadly in other cortical areas, while V1-enriched subtypes are exclusive to the occipital cortex. g, Z-score spatial maps of V1- and V2-enriched genes. h, Venn diagrams for the overlap between V2-enriched (left) and V1-(right) DEGs across the four pairs of subtypes in c. Only strong DEGs with Log fold change >0.5 and adjusted p-value <0.0001 were considered. See Supplementary Table 7. i, Spatial graph (left) and UMAP (right) from Visium analysis show clear V1-V2 border across all cortical layers and two parallel UMAP trajectories for V1 and V2 that bifurcated from the SP. Capture area, 6.5 x 6.5 mm. j, Spatial graphs showing additional genes identified by Visium exhibiting clear transition at the V1-V2 border. k, Imputed expression patterns for additional V1-V2 border genes matches with Visium results. l, Constellation plots of cell types in different cortical areas show V1 and V2 share a developmental lineage that is distinct with PFC. Matching nodes between V1 and V2 have connection fraction > 0.2 except for EN-IT-L3 and EN-IT-L4. See Supplementary Table 8. m, Gene ontology (GO) analysis reveals synapse and cell adhesion-associated biological processes are upregulated in V1-specific Layer 4 neurons. n, Heatmap for incoming signaling pathways shows neurexin (NRXN) signaling is specifically enriched in EN-IT-L4-V1 cluster comparing to other upper layer (UL) EN-IT clusters. o, Network heatmap shows EN-IT-L4-V1 is unique among upper layer EN-ITs to receive NRXN signaling from both IN and EN sources, leading to significant increase in overall interaction strength. The bar graphs on the top and right side are sum of interaction strength as incoming and outgoing signals, respectively.