A Single-Cell Transcriptomic Atlas of Human Neocortical Development during Mid-gestation

Abstract

We performed RNA sequencing on 40,000 cells to create a high-resolution single-cell gene expression atlas of developing human cortex, providing the first single-cell characterization of previously uncharacterized cell types, including human subplate neurons, comparisons with bulk tissue, and systematic analyses of technical factors. These data permit deconvolution of regulatory networks connecting regulatory elements and transcriptional drivers to single-cell gene expression programs, significantly extending our understanding of human neurogenesis, cortical evolution, and the cellular basis of neuropsychiatric disease. We tie cell-cycle progression with early cell fate decisions during neurogenesis, demonstrating that differentiation occurs on a transcriptomic continuum; rather than only expressing a few transcription factors that drive cell fates, differentiating cells express broad, mixed cell-type transcriptomes before telophase. By mapping neuropsychiatric disease genes to cell types, we implicate dysregulation of specific cell types in ASD, ID, and epilepsy. We developed CoDEx, an online portal to facilitate data access and browsing.

Keywords: autism; cortical development; differentiation; epilepsy; evolution; human; intellectual disability; neurogenesis; schizophrenia; subplate.

Figure 1.. A catalog of cell types in developing human neocortex.

(A) Schematic illustrating experimental design and anatomical dissections. VZ: ventricular zone; iSVZ: inner subventricular zone; oSVZ: outer subventricular zone; IZ: intermediate zone; SP: subplate; CPi: inner cortical plate; CPo: outer cortical plate; RG: radial glia; IP: intermediate progenitor; MN: newborn migrating excitatory neuron; EN: excitatory neuron; IN: interneuron; O: oligodendrocyte precursor; E: endothelial cell; P: pericyte; M: microglia. (B) Scatter plot visualization of cells after principal components analysis and t-stochastic neighbor embedding (tSNE), colored by Seurat clustering, and annotated by major cell types. (C) Heatmap of gene expression for each cell. Cells are grouped by Seurat clustering, and the mean expression profile of enriched genes for each cluster was used to hierarchically cluster the Seurat clusters. The top 20 most enriched genes are shown per cluster, and anatomical marker genes in the top 20 are noted. Color bar matches Seurat clusters in B. (D and E) tSNE of cells colored by anatomical source (D), or mean expression of groups of canonical marker genes of major cell types (E). (F) Heatmap of expression profiles of canonical cell type marker genes. Cells are grouped by Seurat clustering. Color bar matches Seurat clusters in B. (G) Cluster metrics. Ratio of cells derived from GZ or CP. Percent of total cell population. Percent of cells derived from each donor. Bar colors indicate grouping of cells by major cell type, e.g. CGE and MGE derived interneurons are both blue. MP: mitotic progenitor. (H) Pseudo-time analysis using Monocle 2.0 of cells expected to be part of the neurogenesis-differentiation axis, colored by Monocle state or pseudo-time. Each point represents a cell. Pseudo-time represents an ordering of cells based upon the inferred trajectory, predicting the lineage trajectory. (I) Pseudo-time trajectory colored by Seurat clusters.

Figure 2.. Sub-clustering analysis identifies progenitor states and subtypes of excitatory and inhibitory cells.

(A) Diagram of sub-clustering analysis workflow. An iterative approach was used, cells from each initial cluster were re-processed, clustered, and analyzed from the raw counts matrix using Seurat. tSNE is colored by Seurat clustering, and annotated by major cell types. (B) Sub-clustering of progenitors. Progenitors separate by cell type and cell cycle state. (C) Sub-clustering of interneurons. InMGE sub-clusters by maturity and cell subtype. InMGE-7 displays enrichment of TAC1, a marker of PV interneurons, and does not express SST (Pfeffer et al., 2013). InMGE-6 shows strong enrichment of NPY and SST. InCGE sub-clusters by maturity. All clusters are CALB2+, with differing levels of expression likely reflecting maturity. (D) Sub-clustering of excitatory neurons. New born (ExN) and maturing excitatory neurons (ExM) sub-cluster by maturity. ExM begin to display separation of laminae markers. The excitatory upper layer enriched cluster (ExM-U) shows enrichment of laminae markers for different sub-clusters, and expression of the callosal marker LMO4 (Molyneaux et al., 2007). The deep layer cluster (ExDp1) separates by layer. ExDp1–2 is enriched for the subplate marker NR4A2 (Hoerder-Suabedissen and Molnar, 2015), ExDp1–0 is enriched for lower L5 and L6 markers (CRYM, TBR1, FOXP2) (Molyneaux et al., 2007), and ExDp1–1 and 3 are enriched for L4 and upper L5 markers (RORB, FOXP1, ETV1) (Ferland et al., 2003; Molyneaux et al., 2007). Heatmaps: Heatmaps of expression profiles by sub-cluster of groups or individual marker genes (Y-axis). The laminae bar indicates the percent of cells derived from the CP. Purple: 100% of cells derive from the CP, 0% GZ; Green: 0% of cells derive from the CP, 100% GZ. Upper layer and deep layer gene groups are the top 50 most enriched genes from the excitatory upper enriched cluster and the deep layer cluster, respectively. tSNES: tSNEs of cells are colored by features of interest: sub-cluster, anatomical source, donor, or gene expression. Grey indicates cells with an undefined transcriptional signature. For heatmaps and tSNE, gene expression is plotted as a z-score for the population of cells in the plot, therefore some cell types display differences in relative expression of cell type markers between sub-clusters of the same major cell type, but all express the marker at some level (e.g. all RG express markers of RG, but some sub-clusters of RG have higher relative expression than other sub-clusters of RG). Labels “mat” and “dif” indicate inferred order of differentiation or maturation.

Figure 3.. Cell type enrichment of TFs and co-factors.

(A) Heatmap of expression of TFs, co-factors, and chromatin remodelers enriched in RG, excitatory neurons, and deep layer excitatory neurons. Cells are grouped by cluster. Red indicates factors previously unknown to be enriched in the neocortical cell types of interest. (B) Expression of factors of interest in bulk tissue LCM laminae from developing cortex. (C) RNA FISH of fetal cortex probed with the newly identified cell-enriched TF ZFHX4 (neural progenitors in the VZ and SVZ), and known markers PAX6 (RG marker) and EOMES (IP marker). Insets show higher magnification of the VZ and SVZ. (E) Quantification of the percentage of PAX6+ or EOMES+ cells co-expressing ZFHX4. ZFHX4 is expressed in both RG and IPs. (F and G) RNA FISH of fetal cortex probed with the newly identified cell-enriched TFs CARHSP1 (neural progenitors in the VZ and SVZ), and CSRP2 (glutamatergic neurons in the CP). (C, D, F, G) Quantification of normalized fluorescence intensity per layer for each set of probes (see materials and methods). Scale bar = 250μm (left) or 100μm (inset). (H) Schematic of cell-type specific expression of factors of interest. Color indicates -log₁₀ p-value from Fisher’s test.

Figure 4.. Characterization of subplate neuron expression profiles.

(A) tSNE colored by Seurat clustering, and annotated by major cell types. (B) tSNE of cells colored by mean expression of groups of marker genes or expression of specific genes. (C) Expression of SP markers in bulk tissue LCM laminae from developing cortex. SP markers were derived from literature sources (left), or by differential expression of the SP versus the VZ, SZ, CP, and MZ and visual confirmation of SP specificity (right). (D) Sub-clustering of the deep layer excitatory cluster 1. tSNE for the full dataset colored by sub-clustering (left). tSNE of cells belonging to the deep layer excitatory cluster (right), colored by sub-clustering, mean expression of groups of marker genes, or expression of specific genes. (E to G) Expression of SP cluster enriched genes (F), ST18 co-expressed genes (G), and the intersection of both F and G (E) in bulk tissue LCM laminae from developing cortex. Genes are ordered left to right by enrichment or correlation (highest left). Light blue text indicates SP markers previously identified. (H) Eigengene of intersected ST18 co-expressed and SP cluster enriched genes (E) plotted in bulk tissue LCM laminae from developing cortex. P-values: * <0.05, ** <0.01, *** <0.001, **** <0.0001. (I) RNA FISH of fetal cortex probed with the newly identified subplate enriched TF ST18. Quantification of normalized fluorescence intensity per layer for each set of probes (see materials and methods). Scale bar = 250μm (left) or 100μm (inset).

Figure 5.. Transcriptional network discovery.

(A to H) Regulatory elements for cell-type specific genes. (A) Enhancer size by cell type. Enhancers are assigned to cell types by cell type enriched genes. (B) Density plot of enhancer sizes that are assigned to specific cell types. (C) Distance (base pairs) from enhancer to promoter by cell type. (D) Density plot of distance (base pairs) from enhancer end to promoter start that are assigned to specific cell types. (E) Enhancers per gene by cell type. (F) Histogram of enhancers per gene that are assigned to specific cell types. (G) Number of enhancers per gene versus CDS length of the gene by cell type. (H) Number of enhancers per gene versus GC content of the gene by cell type. (I) Schematic showing the computational approach used for transcriptional network discovery with the SCENIC pipeline (see materials and methods). 1) Co-expression modules between transcription factors and candidate genes are constructed. 2) Genes in co-expression modules are then pruned to genes which are inferred to be direct targets of the transcription factor, making a regulon. Direct targets are determined by the presence of the transcription factors binding motif in the regulatory elements associated with that gene. 3) The activity of each regulon is then assessed in each cell. (J) Cell type enrichment of regulon activity. Each regulon was scored as active or inactive for each cell, and cluster enrichment was then determined by Fisher’s test. Color indicates FDR-corrected -log₁₀ p-value. (K) SCENIC regulon activity in each cell (AUCell) for the indicated TF plotted on tSNE. (L) TFs with previously uncharacterized cell type or cell subtype specific activity. Regulon activity in each cell (AUCell) for the indicated TF (top panels) or expression of the TF plotted on tSNE (bottom panels).

Figure 6.. Dissecting the acquisition of a neuronal program.

(A) tSNE colored by mean expression of cell cycle phase markers. (B) tSNE colored by co-expression of groups of canonical cell type markers. Yellow indicates co-expression. (C) Percent of cells in each Seurat cluster displaying co-expression of major cell type markers. (D) Mixed transcriptomic signatures of mixed marker cells in S-phase corresponding to the expression of markers from multiple cell types. For the RG to IP comparison, RG and IP eigengenes were derived from differentially expressed genes between RG and IP cells, and similarly for the RG to Neuron comparison and IP to Neuron comparison. Boxplots: box indicates first and third quartiles; the whiskers extend from the box to the highest or lowest value that is within 1.5 * inter-quartile range of the box; and the line is the median. (E) Shared gene signatures between major cell types and mixed cell types. Overlap of gene signatures from major cell types (y-axis), and genes differentially expressed between major cell types and mixed marker cells (labeled on the grey bar). X-axis: percentage of genes differentially expressed between major cell types that are also differentially expressed between the corresponding major cell type and mixed marker cells. For example, ~85% of IP signature genes are more highly expressed in RG+IP+ cells than RG+ cells. (F) RNA FISH of fetal cortex probed with the S-phase marker PCNA (green), the RG marker PAX6 (red), the neuron marker STMN2 (magenta), and stained with DAPI (blue). Panels on the right show high magnification single-plane confocal images of individual cells expressing all three markers. Scale bar = 100μm (left) or 10μm (right). (G) Quantification of the percentage of cells co-expressing the S-phase marker PCNA, the RG marker PAX6 and the neuron marker STMN2. (H) Quantification of relative amounts of mitotic RG and relative amounts of IPs undergoing different differentiation events. (I) Diagram of mixed cell type transcriptomic states that is characteristic of neurogenic differentiation trajectories in human neocortex. P-values: * <0.05, ** <0.01, *** <0.001, **** <0.0001.

Figure 7.. Cellular determinants of disease.

(A to C) Cell type expression of ASD, epilepsy, or ID risk genes respectively. Expression of ASD risk genes is enriched in fetal glutamatergic neurons with some genes specifically expressed in other cell types. Red: gene is discussed in text. Cells are ordered by cluster. (D) Cell type enrichment of ASD, epilepsy, or ID risk genes. Numbers indicate log₂ odds ratio, the red line indicates FDR-significance threshold (p-value 0.05).

Figure 8.. Partitioned heritability analysis demonstrates enrichment of heritability in specific brain traits and neuropsychiatric diseases in diverse cell types.

(A) Schematic showing the approach to identify regulatory elements (RE) for specific cell types and assess enrichment for specific brain traits. REs of genes enriched in specific cell types are identified by chromatin accessibility correlation between the promoter of the gene and other accessible peaks within 1Mb. The set of promoter and distal RE peaks are then tested for enrichment in SNPs associated with brain traits and neuropsychiatric disease using partitioned heritability by LD score regression. (B) Heatmap showing significant partitioned heritability enrichment for specific brain traits and neuropsychiatric disorders in different cell populations. Color indicates the partitioned heritability enrichment. Numbers are the FDR-corrected p-values. References for each GWAS are in Table S8. We did not observe enrichment of IBD or finger whorl variants in the regulatory elements of any of the cortical derived cell types, supporting the cell type specificity of gene regulation. (C) For selected GWAS, barplots indicate the FDR-corrected significance or the enrichment (right) of partitioned heritability. Red vertical line indicates FDR-significance threshold (p-value 0.05). Error bars represent standard error. N: GWAS sample size.