Spatiotemporal gene expression trajectories reveal developmental hierarchies of the human cortex

See allHide authors and affiliations

Science  08 Dec 2017:
Vol. 358, Issue 6368, pp. 1318-1323
DOI: 10.1126/science.aap8809

Building a brain

The human brain is built in an inside-out manner as a series of layers. Although progenitor cells spin off new neurons in a seemingly organized fashion, the devil is in the details. Nowakowski et al. analyzed the transcriptomes of single cells from the developing brain to elucidate the hidden complexity of brain construction. For each cell, its position within the brain matters, as well as what type of neuron is being made at what point during overall development. These individual expression patterns result in organized diversity in the brain's cortex.

Science, this issue p. 1318


Systematic analyses of spatiotemporal gene expression trajectories during organogenesis have been challenging because diverse cell types at different stages of maturation and differentiation coexist in the emerging tissues. We identified discrete cell types as well as temporally and spatially restricted trajectories of radial glia maturation and neurogenesis in developing human telencephalon. These lineage-specific trajectories reveal the expression of neurogenic transcription factors in early radial glia and enriched activation of mammalian target of rapamycin signaling in outer radial glia. Across cortical areas, modest transcriptional differences among radial glia cascade into robust typological distinctions among maturing neurons. Together, our results support a mixed model of topographical, typological, and temporal hierarchies governing cell-type diversity in the developing human telencephalon, including distinct excitatory lineages emerging in rostral and caudal cerebral cortex.

The human cerebral cortex consists of billions of neurons that emerge from an initially uniform neuroepithelium. Radial glia generate neurons during a restricted window of neurogenesis, with temporal and spatial factors influencing the fate of daughter neurons. Excitatory neurons migrate to the cortical plate, forming six cortical layers with stereotyped connectivity patterns that result in functional circuits. Across the cortical sheet, dozens of functionally distinct areas share similar cell types and lamination patterns (13). Nonetheless, differences between functional areas have long been appreciated (4) and are hypothesized to emerge from a molecular protomap along the germinal zone during neurogenesis (57). However, the extent to which early transcriptomic differences establish distinct laminar and areal cell types has been difficult to determine in bulk analysis of heterogeneous tissue samples. We used single-cell profiling to identify temporally and spatially restricted trajectories of neuronal differentiation. Our results suggest that rudimentary distinctions between laminar subtypes emerge early in development and that areal identity strongly influences distinctions among excitatory neurons, supporting a cell type– rather than lamina-centric view of cortical development.

We performed single-cell mRNA sequencing (scRNA-seq) in primary cortical and medial ganglionic eminence (MGE) samples across stages of peak neurogenesis (tables S1 to S3). To capture cells along the span of neuronal differentiation from progenitor cells to postmitotic neurons, we included samples of microdissected germinal zone and cortical plate. In addition, to compare cells of the same type across distinct cortical areas, we captured cells from the prefrontal cortex (PFC) and primary visual cortex (V1), including 13 paired specimens (Fig. 1A and figs. S1 and S2). Using unbiased clustering, we identified transcriptionally distinct cell clusters corresponding to 11 broad categories: astrocytes, oligodendrocyte precursor cells, microglia, radial glia, intermediate progenitor cells, excitatory cortical neurons, ventral MGE progenitors, inhibitory cortical interneurons, choroid plexus cells, mural cells, and endothelial cells (Fig. 1, B to D; figs. S2 to S4; and tables S4 and S5). Cells clustered largely by type rather than individual or batch (fig. S2), and major groups were composed of cells from multiple individuals (Fig. 1, E and F). To resolve finer distinctions among cells and to relate gene expression variation to cell type, developmental stage, and anatomical source, we performed weighted gene coexpression network analysis (WGCNA) (fig. S3 and tables S6 and S7) (8). We reasoned that although individual genes may not be reliably detected by using scRNA-seq because of high dropout rate, gene coexpression networks would be more robustly detected and could more accurately reflect the molecular variation related to distinct biological processes.

Fig. 1 Diverse cell types in human telencephalon development.

(A) Schematic illustrating sample collection over time, region, and lamina. (B to D) Scatterplot of 4261 cells after principal components analysis and t-stochastic neighbor embedding (tSNE), colored by (B) cluster, (C) cortical lamina source, and (D) maker gene expression. (E) Similarity matrix representing pairwise Pearson correlations between clusters across principal components. (F) Metadata describing each cluster. VZ, ventricular zone; OSVZ, outer subventricular zone; GZ, germinal zone; CP, cortical plate; V1, primary visual cortex; PFC, prefrontal cortex; MGE, medial ganglionic eminence; LGE, lateral ganglionic eminence; M1, primary motor cortex.

We first focused on variation among radial glia, which act as neural stem cells. During early development, the telencephalic neuroepithelium becomes regionalized by a small number of patterning molecules produced at local signaling centers (9), which induce topographically distinct patterns of transcription factor expression in dorsal telencephalic (excitatory cortical neuron producing) and ventral telencephalic (inhibitory cortical interneuron producing) radial glia (Fig. 2A). Unbiased clustering highlighted a robust distinction between dorsal and ventral telencephalic radial glia, including reciprocal expression of patterning transcription factors EMX1 and NKX2-1 (Fig. 2B). By contrast, dorsal telencephalic radial glia sampled from PFC and V1 intermingled across multiple clusters (fig. S5). Thus, based on transcriptomic signatures, radial glia follow a topographical hierarchy to distinguish dorsal and ventral telencephalic lineages, but within the dorsal telencephalon, radial glia appear to follow a typological hierarchy.

Fig. 2 Maturation of radial glia across germinal zones.

(A) Schematic of early radial glia patterning. (B) Heatmap of gene expression in radial glia sampled from human telencephalon, organized according to cluster. (C) WGCNA reveals age-correlated modules combined into a pseudoage score. pcw, post-conception weeks. (D) PFC and V1 radial glia show a systematic difference in maturation. (E) Trajectories of gene expression across cortical radial glia pseudoage. (F) Landscapes of inferred gene expression along pseudolamina from VZ to SVZ (y axis) and pseudoage (x axis) for genes with conserved and zone-specific trajectories. (G) Immunohistochemistry for NEUROD1 in primary tissue sections. (H) Landscape plots of mTOR modulators. (I) Detection of mTOR effector pS6 in oRG cells. (J) Model for bifurcation of radial glia subtypes during maturation.

Within each sample, dorsal telencephalic radial glia may include cells at different stages of maturation. For example, the human PFC matures several weeks earlier than V1 (10, 11). To deconvolve this axis of biological variation by using transcriptomic profiles, we combined age-correlated gene coexpression networks to define a maturation score, or “pseudoage.” Radial glia pseudoage scores correlated strongly with the sample age (Fig. 2C and figs. S5 and S6). Remarkably, pseudoage trajectories between PFC and V1 radial glia displayed a systematic divergence of ~3 weeks (Fig. 2D), revealing accelerated maturation of PFC radial glia along a common molecular program. Analysis of gene expression trajectories across pseudoage further revealed that early human cortical radial glia share a conserved gene network related to Wnt signaling (12), whereas genes involved in gliogenesis are up-regulated later in development (Fig. 2, E and F). Unexpectedly, early human radial glia also show enriched expression of proneural transcription factors, including NEUROG2, NEUROD6, and NEUROD1 (Fig. 2G and fig. S5).

Ventricular radial glia (vRG) of the cortex have bipolar morphology, with fibers spanning the thickness of the cerebrum. These early vRG cells differentiate to generate pia-contacting outer radial glia (oRG) that delaminate from the ventricular zone (VZ) and translocate to the outer subventricular zone (OSVZ) and truncated radial glia (tRG) that remain in the VZ (13). To compare radial glia in these laminae across maturation stages, we combined gene coexpression networks correlated with radial glia sampled from microdissected VZ and OSVZ to define a lamina score, or “pseudolamina” (Fig. 2F and fig. S6). Gene coexpression networks correlated with pseudolamina and pseudoage not only supported this classification of radial glia subtypes (14, 15) but also highlighted selective activation of signaling pathways. For example, DAB1, which mediates reelin signaling (16), is enriched in vRG and oRG cells, whose processes contact reelin-expressing cells in the pia, but depleted in tRG cells, whose processes are truncated (Fig. 2F). In addition, we found GLI2 enrichment in young oRG cells, which may relate to the proposed role of sonic hedgehog signaling in oRG cell specification (17), and we observed NFAT2C enrichment in ventricle-contacting cells, which is consistent with the reported roles of calcium signaling among gap junction–coupled vRG cells (18). Previous studies modulating mammalian target of rapamycin (mTOR) signaling in mouse developing cortex reported subtle phenotypes related to proliferation and neurogenesis (19), yet mutations in genes related to this pathway have strong neurodevelopmental phenotypes in humans, including autism and megancephaly (20). Surprisingly, we found that oRG cells showed enriched expression of several regulators of the mTOR signaling pathway and increased phosphorylation of the S6 ribosomal protein, a canonical readout of activated mTOR signaling, compared with the closely related vRG, tRG, and intermediate progenitor cells (Fig. 2, H and I, and fig. S7). This finding suggests that during peak stages of neurogenesis, oRG cells may be especially vulnerable to mTOR pathway gene mutations. Along with our previous observations (21), our results indicate that multiple signaling pathways display heterogeneous activation during radial glia diversification (Fig. 2J).

In addition to maturation, major gene expression changes also occur during neuronal differentiation (22). Distinct neuronal cell types have been proposed to arise through divergence of serial sister lineages from ancestral differentiation programs (23). Single-cell sequencing enables the disambiguation of these cell type–specific differentiation trajectories across related lineages. To identify a common telencephalic differentiation program, we combined gene coexpression networks across cells from the cortical and MGE neuronal lineages to generate a single consensus neuronal differentiation network (Fig. 3A and table S8). We found that gene expression trajectories are not only highly conserved across cortical PFC and V1 areas (fig. S8) but also retain striking similarity between more divergent dorsal and ventral lineages (Fig. 3B and fig. S8). Nonetheless, known and previously unidentified markers also showed lineage-specific trajectories (Fig. 3C) that may guide specialized aspects of lineage-specific neurogenesis and could underlie divergence among serial sister lineages.

Fig. 3 Landscapes of neuronal differentiation.

(A) tSNE plots colored by module eigengene expression for topographically conserved neuronal differentiation modules. (B) Scatterplot of all genes for correlation with conserved differentiation network across MGE lineage (y axis) and cortical lineage (x axis) and expression of lineage-specific modules. (C) Scatterplot of excitatory lineage cells along axes of differentiation and maturation. Labeling cells by metadata reveals undifferentiated cell enriched for GZ origin, differentiated cells enriched for CP origin, and the expected age distribution. (D) Landscapes of inferred gene expression and schematic illustrate temporally conserved and temporally restricted trajectories of neuronal differentiation.

Current methods for developmental lineage reconstruction based on gene expression reorder cells along a single pseudotime axis that may conflate signatures of differentiation with signatures of temporal maturation. However, temporal cues during cortical development are crucial for generating cellular diversity (24). Therefore, we sought to infer gene expression trajectories across both biologically relevant axes of variation by scoring cells separately along pseudodifferentiation and pseudoage axes (Fig. 3C and fig. S8D). These transformations recovered known gene expression trajectories conserved throughout neurogenesis (for example, GLI3 to EOMES to STMN2 in cortex), as well as parallel trajectories of differentiation that segregate by maturation stage (Fig. 3E and fig. S8). Together, these methods for lineage reconstruction enable the inference of gene expression trajectories from heterogeneous developmental tissue in which parallel branches of differentiation are intermingled.

Recent studies suggest that adult cortex consists of at least 8 to 19 distinct excitatory neuron subtypes distributed across six layers (2527). We sought to investigate the extent to which discrete subtypes emerge during neurogenesis before the onset of sensory experience. During peak cortical neurogenesis, the timing of neuronal cell birth predicts its terminal laminar position, although the expression patterns of cortical layer markers sometimes follow complex developmental trajectories (28). Clusters of newborn neurons and laminar markers were sequentially enriched across developmental time points, which is consistent with the temporal hierarchy of neuronal subtype generation (Fig. 4, A to D). This finding suggests that rudimentary molecular signatures that distinguish neurons generated at different time points manifest at the earliest stages of neuronal differentiation. These differences may reflect prespecification of neurons inherited from parent progenitors, or nonautonomous cues mediated through the local environment, such as molecularly distinct radial glia fibers supporting their radial migration (13).

Fig. 4 Emergence of area-specific neuronal cell clusters.

(A to C) tSNE plot highlighting excitatory neurons, colored by (A) cluster, (B) area of origin, and (C) sample age. (D) Frequency of cells from each cluster and expressing common layer markers plotted over pseudoage. (E) Collection from PFC and V1, including 13 paired samples, and example genes with area-specific expression. Cells are colored by area of origin. (F) Validation of area specificity. (G and H) Examples and overall number of differentially expressed genes across cell types. (I) Genes plotted by correlation with neuronal differentiation versus area. (J) Area-specific subtypes emerge as neurons mature. (K) Density plot of cells colored based on coexpression of SATB2 and BCL11B in PFC, but not V1. (L) Narrow transition zone marks shift to rostral coexpression.

We next evaluated whether further refinement of neuronal subtypes occurs as neurons mature. Among maturing neurons, we observed clusters containing cells captured mostly from young samples (<15.5 weeks after conception) that showed strong enrichment for markers of subplate and deep cortical layers. In addition, a second set of clusters contained a mixture of cells from young and old samples, suggesting their translaminar distribution, and showed enrichment for upper-layer markers (Fig. 4, A to C). However, genes distinguishing these maturing neuron clusters showed limited correspondence to adult cortical neuron type or layer markers (fig. S9 and tables S9 and S10) (2729), suggesting that the coarse laminar distinctions among newborn neurons (Fig. 4D) are not resolved into refined transcriptional states until later in development.

We further examined whether additional factors beyond laminar subtype shape variation between maturing excitatory neurons. Surprisingly, distinct clusters of maturing neurons reproducibly emerge according to the topographical location from which they were sampled (Fig. 4E). These area-specific clusters were enriched for previously validated markers, including KCNJ6 in PFC and TENM2 in V1 (6, 30), but also revealed previously unidentified candidates restricted to distinct subtypes of PFC and V1 neurons (figs. S10 and S11). For example, CPNE8 expression is restricted to the upper layers of PFC, whereas HCRTR2 is restricted to the subplate of PFC (Fig. 4F). Some of these patterns that distinguish individual clusters are not reflected among markers of neuronal types or layers of the adult cortex, suggesting that transcriptional states present transiently during development may influence area-specific features of neuronal identity.

According to the classic protomap hypothesis (5), molecular differences among progenitor cells parcellate the emerging cortical sheet into distinct areas. We found significant gene expression differences between PFC and V1 across cells of every type, but the number of differences increased along the axis of excitatory neuron differentiation, with the greatest number of differences among maturing excitatory neurons (Fig. 4, G and H; fig. S10; and table S11). Together with results from laminar dissections (6), these findings suggest that a limited number of area-specific gene expression differences among radial glia along a protomap cascade during neuronal differentiation into molecular differences that establish distinct neuronal subtypes at early stages of neuronal maturation (Fig. 4, I and J).

Our study suggests a new model of neuronal-type hierarchy during development, in which typological distinctions between maturing neurons appear to strongly relate to topographical hierarchy in the cortical sheet during peak neurogenesis. To further examine the notion of typological distinction between PFC and V1 enriched clusters, we compared distinct coexpression profiles of genes involved in the specification of projection patterns during development. For example, maturing PFC neurons expressing BCL11B, a marker of subcortical projections, coexpressed SATB2, a marker of intracortically projecting neurons at all stages examined (Fig. 4, K and L, and fig. S11). These cells may represent a neuronal population with a dual projection pattern (31). Only at later stages of development did neurons that express SATB2 but not BCL11B arrive at the cortical plate in the PFC. In contrast, neurons in V1 segregated the expression of these genes as early as 13.3 weeks after conception (fig. S11). Surprisingly, we observed a narrow transition zone between the caudal pattern of discrete expression and the rostral pattern of marker overlap approximately midway through the rostro-caudal extent of the cortical sheet (Fig. 4L). This narrow transition zone suggests that cell-type inferences from scRNA-seq profiling of frontal and occipital lobes do not simply reflect the opposite ends of a developmental gradient but instead represent distinct excitatory cortical neuron lineages established in rostral and caudal telencephalon.

Early studies of functional connectivity and neuroanatomy(1, 2) suggested that the cerebral cortical “column” consists of similar cell types and lamination patterns irrespective of topographical location in the cortical sheet. This perspective shaped the highly influential view of cortical uniformity. Although gene expression differences across cortical areas have been revealed through a number of studies [reviewed in (3)], the extent to which these differences establish distinct neuronal cell types has been difficult to determine in bulk analysis of heterogeneous tissue samples. We used scRNA-seq to survey cells from primary samples of the human telencephalon during early stages of area specification and layer formation in the cerebral cortex. By relating gene expression to distinct biological axes of variation—maturation, differentiation, and area (fig. S1)—we identified conserved gene networks involved in neurogenesis, as well as lineage, area, and developmental stage-specific trajectories that contribute to the diversification of cortical neurons. These gene expression signatures may highlight patterns of selective vulnerability for neurodevelopmental disorders, including autism (fig. S12), and are available via an interactive browser ( (fig. S13). Transcriptomic profiles extend the current understanding of developmental hierarchies in the developing forebrain: A topographical hierarchy segregates dorsal and ventral telencephalic lineages of developing cortical neurons; within the dorsal telencephalon, temporal and typological differences define progenitors across cortical areas, while topographical distinctions predominate across maturing neurons.

Supplementary Materials

Materials and Methods

Figs. S1 to S13

References (3238)

Tables S1 to S11

References and Notes

Acknowledgments: The authors thank S. Wang, O. Meyerson, R. N. Delgado, A. Diaz, J. H. Lui, M. Paredes, E. Huang, W. Walantus, and members of the A.R.K. laboratory for providing research resources, technical help, and helpful comments during manuscript preparation. RNA-sequencing data has been deposited at dbGaP: phs000989.v3. This study was supported by the Damon Runyon Cancer Research Foundation (grant DRG-2166-13 to A.A.P.), NIH awards U01 MH105989 and R35NS097305 to A.R.K., R01NS091544 to D.A.L., and by a gift from B. Osher. M.H. and J.K. are funded by California Institute for Regenerative Medicine grants GC1R-06673-C and GC1R-06673-C, respectively. J.A.W. is an inventor on patent application 15/055,252, 20150337295, 20130302884, 20130302883, 20130302807, 20130296196, 20130295602, 20160025761 held/submitted by Fluidigm that covers single-cell nucleic acids for high-throughput sudies, methods and devices for analysis of defined multicellular combinations, and methods, systems, and devices for multiple single-cell capturing and processing using microfluidics. The supplementary materials contain additional data. T.J.N. and A.A.P. designed the study. Single-cell experiments were performed by T.J.N., A.A.P., B.A., J.S., A.A.L., C.S.E., and E.D.L., with input from D.A.L. and J.A.W. Validation experiments and immunohistochemistry were performed by T.J.N., M.A.M.R., J.R.O., and E.D.L. Data management and read processing was performed by A.B., with initial contribution from D.V. Data analysis was performed by A.B., T.J.N., C.S.E., S.J.L., A.A.P., and X.W. Data viewer was designed by M.H. and W.J.K. The study was supervised by T.J.N., A.A.P., A.B., and A.R.K. This manuscript was prepared by T.J.N., A.B., and A.A.P., with feedback from all authors.

Stay Connected to Science

Navigate This Article