Molecular features driving cellular complexity of human brain evolution

Abstract

Human-specific genomic changes contribute to the unique functionalities of the human brain^1-5. The cellular heterogeneity of the human brain^6,7 and the complex regulation of gene expression highlight the need to characterize human-specific molecular features at cellular resolution. Here we analysed single-nucleus RNA-sequencing and single-nucleus assay for transposase-accessible chromatin with sequencing datasets for human, chimpanzee and rhesus macaque brain tissue from posterior cingulate cortex. We show a human-specific increase of oligodendrocyte progenitor cells and a decrease of mature oligodendrocytes across cortical tissues. Human-specific regulatory changes were accelerated in oligodendrocyte progenitor cells, and we highlight key biological pathways that may be associated with the proportional changes. We also identify human-specific regulatory changes in neuronal subtypes, which reveal human-specific upregulation of FOXP2 in only two of the neuronal subtypes. We additionally identify hundreds of new human accelerated genomic regions associated with human-specific chromatin accessibility changes. Our data also reveal that FOS::JUN and FOX motifs are enriched in the human-specifically accessible chromatin regions of excitatory neuronal subtypes. Together, our results reveal several new mechanisms underlying the evolutionary innovation of human brain at cell-type resolution.

Extended Figure 1:. Annotation and quality control of single-nuclei RNA-seq and single-nuclei ATAC-seq.

(A) Distribution of sex and humanized age of samples. (B) Broadly annotated UMAP of nuclei per species. (C) Total nuclei number per sample after filtering. (D) Normalized, log (ln) transformed expression values of major cell type markers. (E) Violin plots of number of detected UMIs (log₁₀ transformed) per major cell type. (F) Percentage of cells contributed per individual per species per major cell type. (G) Broad annotation of snATAC-seq data per species. (H) Total nuclei number per sample after quality control. (I) Nucleosome band pattern per sample; each line represents one sample. First, second and third peaks represent nucleosome free, mononucleosome and dinucleosome fractions, respectively. (J) Percentage of cells contributed per individual per species per major cell type. (K) Clarity of annotation transfer from snRNA-seq to snATAC-seq as displayed by Jaccard similarity index, which is the number of nuclei with the same final annotation and prediction (intersection) divided by the total number of nuclei with a given annotation or prediction (union). y-axis represents final annotation; x-axis represents the prediction which was assigned by label transfer per nucleus. Higher values indicate more similarity between final annotation and initial prediction. (L) Fraction of reads in peaks per sample (N = 9280, 5383, 5657, 4655, 5941, 4381, 5691, 4690, 3321, 6426, 5984, 6793 left to right). Boxplots represent median and interquartile range.

Extended Figure 2:. Annotation of oligodendrocyte lineage cells.

(A) UMAP visualization of integrated and annotated oligodendrocyte lineage nuclei in snRNA-seq. Oligodendrocyte: mature oligodendrocytes, COP: committer oligodendrocyte progenitor cells, OPC: oligodendrocyte progenitor cells. (B) Percentage of nuclei per sample for each subtype in snRNA-seq. (C) Normalized and scaled (z-scored) expression values of major oligodendrocyte lineage cell type markers. (D) UMAP visualization of annotated oligodendrocyte lineage nuclei in snATAC-seq per species. (E) Clarity of annotation transfer from snRNA-seq to snATAC-seq as displayed by Jaccard similarity index (similar to Extended Figure 1K). (F) Percentage of cells contributed per individual per species per major cell type. (G-H) smFISH of PDGFRA (OPC) and MOG (MOL) in humans (G) and chimpanzees (H) (region: posterior cingulate cortex. Images span all cortical layers in both species. Scale bar is 100 μm). Similar results have been obtained for 4 bins across 2 humans and for 6 bins across 3 chimpanzees (see Extended Figure 4C-D).

Extended Figure 3:. Integration and annotation of neurons.

(A) Annotated UMAP of excitatory neurons integrated across all species in snRNA-seq. (B) Percentage of nuclei per sample for each excitatory subtype in snRNA-seq. (C) Annotated snATAC-seq per species in the UMAP space. All 14 subtypes identified in snRNA-seq are also distinctly annotated in snATAC-seq for all species. (D) Percentage of nuclei per sample for each excitatory subtype in snATAC-seq. (E) Excitatory subtype markers for validation of annotation (expressions are normalized and log transformed). Note that the individual cells are plotted for C1QL3 since the expression level is not sufficient for a violin plot. (F) Clarity of annotation transfer from snRNA-seq to snATAC-seq as displayed by Jaccard similarity index (similar to Extended Figure 1K). (G) Annotated UMAP of inhibitory neurons integrated across all species in snRNA-seq. (H) Percentage of nuclei per sample for each inhibitory subtype in snRNA-seq. (I) Known inhibitory subtype markers for validation of annotation. Expression levels are normalized and log transformed. Note that the individual cells are plotted for NMBR, PAX6, SYT6 since the expression level is not sufficient for a violin plot. (J) Annotated snATAC-seq per species in the UMAP space. All 8 subtypes identified in snRNA-seq are also distinctly annotated in snATAC-seq for all species. (K) Percentage of nuclei per sample for each inhibitory subtype in snATAC-seq. (L) Clarity of annotation transfer from snRNA-seq to snATAC-seq as displayed by Jaccard similarity index (similar to Extended Figure 1K).

Extended Figure 4:. Additional analyses on the oligodendrocyte lineage.

(A) UMAP of MOLs and OPCs in the anterior cingulate cortex (ACC). (B) Percentage of cells contributed per individual per species per cell type. (C-F) Fractions of MOLs and OPCs in smFISH experiments per section (see Figure 1). Stitched column images encompassing all layers were divided into 5 equal parts from upper (Section 1) to lower layers (Section 5) in all images from human and chimpanzee. (C-D) are data from posterior cingulate cortex (PCC), and (E-F) are data from ACC. Each data point is a bin that contains sections from all layers. C-D: 4 bins from 2 humans, 6 bins from 3 chimpanzees. E-F: 9 bins from 3 humans and 3 chimpanzees. Data are represented as mean values +/− SEM. (G) Deconvoluted proportions from OLIG2+ bulk RNA-seq dataset (reference datasets: (left) chimpanzee, (right) rhesus macaque from this study). N = 22 (human), 10 (chimpanzee), 10 (rhesus macaque) individuals. P-value: Wilcoxon rank sum test, two-sided). (H-I) Fraction of OPCs or MOLs in glia in (H) caudate nucleus and (I) dentate gyrus per species. Each dot represents a sample (p-value: Wilcoxon rank sum test, two-sided. Caudate nucleus: N = 6 per species. Dentate gyrus: N = 6 for human, 3 for rhesus macaque). Box plots represent median and interquartile range in panels G-I. (J) Number of species-specific regulatory changes (PCC snRNA-seq (top), ACC snRNA-seq (middle), and PCC snATAC-seq (bottom, log₁₀ transformed for better readability). (K) Distributions of UMIs (unique molecular identifiers) in ACC and PCC oligodendrocyte lineage nuclei (N=12 individuals both for PCC and ACC). Box plots represent median and interquartile range. (L) Enrichment results between species-specifically expressed genes in ACC (x-axis) and PCC (this study, y-axis). Blue asterisk indicates a significant overlap(FDR < 0.05, Fisher’s exact test, one-sided).

Extended Figure 5:. Additional analyses of the regulatory changes in neuronal subtypes.

(A-B) Number of regulatory changes that are human-specific, chimpanzee-specific or differential between rhesus macaque - human and rhesus macaque – chimpanzee in (A) snRNA-seq or (B) snATAC-seq (log₁₀ transformed for better readability). (C) Scatter plots of number of HS-Genes and CS-Genes per neuronal subtype. Dashed rectangles indicate the subtypes with an excess number of human-specific regulatory gene expression changes (Two-sided chi-square test, FDR < 0.05). Shaded area indicates 95% confidence interval around the best fit (R indicates Spearman’s rank correlation coefficient). (D) Same as (C) for HS-CREs and CS-CREs identified in snATAC-seq data. (E) Percentage distribution of excitatory HS-Genes that are found in only one subtype or shared among increasing number of subtypes (x-axis). Sum of all percentages equal 100. From left to right: in excitatory snRNA-seq, excitatory snATAC-seq, inhibitory snRNA-seq¸ inhibtory snATAC-seq.

Extended Figure 6:. Comparisons of neuronal expression patterns between this dataset and previous comparative bulk datasets.

(A-C) Enrichments of species-specific expression patterns between this study and previous bulk studies between excitatory neurons (left) and inhibitory neurons (right). (A) Transcriptomic changes between the Kozlenkov et al. dataset and this dataset, (B) epigenomic changes between the Kozlenkov et al. dataset and this dataset, (C) transcriptomic changes between the Berto et al. dataset and this dataset. FDR values are from a Fisher’s exact test with multiple testing correction. (D-E) Subtype-specific changes are captured less in the bulk RNA-seq datasets. (D) Comparison of excitatory HS-Genes between a previous bulk analysis and this dataset. Top: odds ratio between the bulk dataset and this dataset with increasing subtype specificity of HS- Genes (from right to left). Bottom: percentage of HS- Genes that were also found in the bulk dataset. (E) Same comparison as (D) with HS-CREs. (F-G) Subtype-specific changes are captured less when the subtypes are combined within the same dataset. (F) Same comparison as (D) with HS-Genes but this time pseudobulk data results are obtained by combining the subtypes in this study. (G) Same comparison as (F) with HS-CREs.

Extended Figure 7:. Associations between HS-Genes and HS-CREs.

(A) The specificity of association between HS-Genes and HS-CREs decreases with increasing distance from the transcription start site (TSS). Y-axis shows the odds ratio, which is defined by the ratio of HS-Genes associated with HS-CREs divided by the ratio of not significant genes (NS-Genes) associated with HS-CREs. We calculated the odds ratio for increasing the distance from the TSS in both directions for four different associations (from left to right): HS-Open-CRE & HS-Up-Gene, HS-Open-CRE & HS-Down-Gene, HS-Close-CRE & HS-Up-Gene, HS-Open-CRE & HS-Down-Gene. The value for each observation was obtained by taking the mean across all cell types. (B-D) Enrichments between HS-CRE associated genes within a 10kb window from the TSS and HS-Genes per cell type. (E-F) Putative target genes of human-specific FOXP2 upregulation in (E) L5-6_THEMIS_1 and (F) L4-6_RORB_2 cells. All genes show human-specific up / down regulation in their respective subtype and reside within 500kb of at least one human-specific chromatin accessibility change that has a FOXP2 motif. Dark blue circles indicate the genes that are not altered in the other 12 excitatory subtypes (similar to FOXP2 itself). Red loop in (A) indicates that FOXP2 itself is also identified with this analysis in the L5-6_THEMIS_1 subtype.

Extended Figure 8:. Further associations between human-specific substitutions and human-specific chromatin accessibility changes.

(A) Pie-chart distribution of published HARs overlapping CREs in this dataset. (B) Ratio of non-BA23 CREs overlapping BA23 CREs (denominator: all CREs in BA23). Each dot represents an independent library prep. Red datasets indicate cortical regions, blue datasets indicate sub cortical regions. (Sample sizes; Superior Middle Temporal Gyri: 8, Middle Frontal Gyrus: 12, Parietal Lobe: 7, Hippocampus: 16, Caudate: 32, Putamen: 11, Substantia Nigra: 14. Box plots represent median and interquartile range). (C) Overlap between cortical HARs (identified in this study) and published HARs (p-value: One-sided chi-square test). (D) Number of HS-CREs associated with a cortical HAR or a published HAR. (E-F) Examples of HS-Open-CRE associated HARs. Bottom panel shows the multi-species alignment for CELF4 HAR. Dots represent no change from the human (hg38) sequence. Human-specific changes conserved in other lineages are highlighted in shaded blue. (G) Enrichment of human-specific substitutions within the HS-CREs per major cell type. Enrichment is tested by a negative binomial regression model with CRE length and evolution of the CRE as the predictor variables (HS-CRE or not HS-CRE) and number of human-specific substitutions as a response variable(Significance: likelihood ratio test). (H) Example of an HS-Open-CRE with many human specific substitutions. (I) Overlap of substitutions that are specific to the human lineage (in comparison to chimpanzee, gorilla and gibbon) and previously identified modern human substitutions. (J) Log fold changes of substitution and HS-CRE association for substitutions on the human (blue boxplots) and modern human lineage (tile red dots) per cell type (except for excitatory cells). Human lineage-specific substitutions were randomly down sampled for 100 times to 12,161 (the number of modern human-specific variants) for comparison. Box plots represent median and interquartile range.

Extended Figure 9:. Supplementary motif enrichment results.

(A-B) Hierarchical clustering of motif enrichments (log-fold change) in HS-Open-CREs across (A) excitatory and (B) inhibitory neuronal subtypes. Transcription factors (TFs) associated with each motif enrichment are displayed in rows and the neuronal subtypes are displayed in columns. Only the motifs enriched in at least one subtype are displayed. (C-D) Accessibility of (C) FOX and (D) FOS / JUN family TFs. Accessibility is assessed by the normalized gene activity scores (calculated using Cicero) per gene per subtype. (E) Annotated UMAP of excitatory neurons in snRNA-seq of surgically resected samples (referred to as PMI0 compared to postmortem BA23 human samples that are referred to as PMI24 in this figure). (F) Percentage of nuclei per sample for each excitatory subtype in snRNA-seq. (G) Enrichments of species-specifically expressed genes when PMI0 or PMI24 datasets were used as the human dataset in the comparative analyses. (H) Pearson correlations (test for p-value is two-sided) between the log fold changes of HS-Open-CRE motif enrichments when PMI0 or PMI24 datasets were used as the human dataset in the comparative analyses. (I) Heatmap of motif FOS / JUN motif enrichments per excitatory subtype in HS-Open-CREs. Colors correspond to −log₁₀(FDR); numbers correspond to log fold change of enrichment.

Extended Figure 10:. Comparisons with external datasets

(A) Expression levels of three ambient RNA markers highly expressed in neurons (SYT1, SNAP25 and NRGN) in the Ma et al. dataset. The dot plot is generated through the interactive web tool linked to the original publication. Dashed square brackets indicate glial cell types, which show exceptionally low levels in the human dataset. Note that the smallest dot shows the presence of a transcript in 40% of the cells. (B) Same as (A) using this PCC dataset. Neuronal ambient RNA markers are detected at very low levels in glial cells across species after ambient RNA removal. (C-E) Enrichment of HS-Genes between the previous study (y-axis) and the current study (x-axis) with two alternative methods. (F) Enrichment of HS-CREs between the previous study (y-axis) and the current study (x-axis) with two alternative methods. For simplicity, we combined all HS-Genes from the subtypes of a major cell type (e.g. all excitatory neuronal subtypes were combined for the excitatory cell type comparisons). P-values were computed using a Fisher’s exact test (one-sided) and false discovery rate (FDR) was calculated per panel.

Figure 1:. The oligodendrocyte lineage is proportionally altered in human evolution.

(A) The fraction of OPC and MOL (mature oligodendrocytes) in glia are altered in humans. Each dot represents a sample (red: snATAC-seq, pink: snRNA-seq. N = 4 individuals per species per assay. P-value: likelihood ratio test, two-sided. See methods). (B-D) smFISH shows increased PDGFRA (OPCs) and decreased MOG (MOLs) signals in humans compared to chimpanzees (region: posterior cingulate cortex). (B) A representative image, scale bar is 100 μm. (C-D) Quantification of the fraction of OPCs and MOLs. Each data point is the average of all subareas in a section (2-4 subareas/5 sections/individual. Human: 10 sections. Chimpanzee:15 sections, see Methods). The p value is the main effect of species from a linear mixed model (random effect: individual, two-sided). Error bars represent SEM. (E) The fraction of OPCs or MOLs in glia based on snRNAseq from anterior cingulate cortex (n = 4 individuals per species, P-value: Wilcoxon rank sum test, two-sided). (F) smFISH similar to (B), but for anterior cingulate cortex. (G-H) Quantification of the fraction of OPCs mOLs as in (C-D) Human: 15 sections, chimpanzee:15 sections, see Methods. (I) Deconvoluted proportions of cells from OLIG2 expressing bulk RNA-seq (reference dataset: human snRNA-seq from this study, n = 22 (human), 10 (chimpanzee), 10 (rhesus macaque) individuals. P-value: Wilcoxon rank sum test, two-sided). (J-K) Fraction of OPC or MOLs in glia per species in the primary motor cortex. (J) Human – marmoset comparison, (K) human – rhesus macaque comparison. N = 5 (human), 4 (marmoset), 3 (rhesus macaque) individuals. P-value: Wilcoxon rank sum test, two-sided. Boxplots represent median and interquartile range in panels A, E, I-K.

Figure 2:. Human-specific gene regulatory changes in the oligodendrocyte lineage.

(A) The proportions of HS (human-specific) regulatory changes to HS + CS (chimpanzee-specific) regulatory changes are higher in OPCs than MOLs (left: snRNA-seq, right: snATAC-seq). P value: chi-square test, two-sided. (B) Same as in (A), except using the anterior cingulate cortex snRNA-seq dataset. (C) Gene ontology enrichment for HS-Down-Genes in OPCs highlights altered cytoskeletal function (p-value: Fisher’s exact test, one-sided). (D) Overlap of HS-Down-Genes in OPCs and primate-conserved COP (committed oligodendrocyte progenitor) markers reveal two COP markers with loss-of-function in human OPCs (p-value: Fisher’s exact test, one-sided). (E-F) Expression levels of (E) SH3RF3 and (F) KIF21A1 across cell types in human, chimpanzee, rhesus macaque. FDR corrected p-values compare the expression levels in OPCs between species (see Supplementary Table 3). (G) snATAC-seq coverage plots of the Human-DOWN-CRE near SH3RF3 in OPCs. TSS: transcription start site. Track scales are the same in all species. (H) UMAP plot of oligodendrocyte lineage cells in mouse adult frontal cortex dataset. NFOL: newly formed oligodendrocytes. (I) Expression pattern of primate-conserved COP markers across mouse oligodendrocyte lineage cell types. Only Sh3rf3 expression is decreased in COPs or NFOLs compared to OPCs. (J) Violin plots of Sh3rf3 and Kif21a expression in mouse oligodendrocyte lineage cell types.

Figure 3:. Subtype and cortical region-specific upregulation of FOXP2 in human neurons.

(A) Expression levels of FOXP2, CNTNAP2 and MET in the posterior cingulate cortex (PCC) show subtype-specific expression changes in the human brain. Human-specific expression is labeled with a red asterisk, x-axis denotes species and excitatory subtypes. (B-D) smFISH of FOXP2 and THEMIS in anterior cingulate cortex (ACC) shows greater number of FOXP2/THEMIS+ cells in human compared to chimpanzee. (B) A representative image. Solid circles show FOXP2 and THEMIS overlapping cells; dashed circles show THEMIS+ cells without FOXP2 expression. Scale bar is 50 μm. (C-D) Quantification of the FOXP2+ cells (C) and (D) FOXP2+ puncta per cell in (left: THEMIS+ cells, right: THEMIS− cells). The p value is the main effect of species from a linear mixed model (random effect: individual, two-sided) with each data point representing a subarea/image/individual. Error bars represent SEM. N = 3 individuals per species. (E-F) FOXP2 is upregulated in the PCC compared to pre-frontal cortex (PFC) and ACC in THEMIS+ neurons (E) but not among all excitatory neurons (F). Y-axis: normalized and log-transformed expression levels. (G) snATAC-seq coverage plots of HS-Open-CREs near FOXP2 in L5-6_THEMIS_1 and L5-6_FEZF2 2-3 neurons. The HS-Open-CREs shown have human-specific chromatin accessibility in L5-6_THEMIS_1 neurons but not in L5-6_FEZF2_2-3 neurons. Track scales are the same in all species.

Figure 4:. Significant association of chromatin accessibility changes with sequence divergence.

(A) Enrichment of publicly available human accelerated regions (HARs) within the HS-CREs. Enrichment is tested by a logistic regression model with CRE length and evolution of the CRE as the predictor variables (HS-CRE or not HS-CRE) and HAR as the response variable (HAR or not HAR, P-value: likelihood ratio test, two-sided). (B) Odds ratio of HAR and HS-CRE association compared to HAR and NS-CRE (non-significant CREs; all CREs that are not HS-CREs) for published HARs (dashed blue line) and HARs identified in this study per log likelihood cutoff (x-axis). Red line indicates the log likelihood cutoff that corresponds to p=0.001 (one-sided). Red dot indicates the odds ratio that corresponds to p=0.001 cutoff. (C) Cortical HAR analysis reveals stronger association with HS-CREs (same enrichment analysis as (A) between cortical HARs and HS-CREs). (D) Modern human-specific variant enrichment within the HS-CREs. Enrichment is tested by a negative binomial regression model with CRE length and evolution of the CRE as the predictor variables (HS-CRE or not HS-CRE) and the number of modern human-specific variants as a response variable (P-value: likelihood ratio test, two-sided). (E) Log fold changes of substitution and HS-CRE association for substitutions on the human (blue) and modern human lineage (tile red) per excitatory subtype. Human lineage-specific substitutions were randomly down sampled for 100 times to 12,161 (the number of modern human-specific variants) for comparison. Empirical p-value (two-sided) is reported for conserved (in blue) and accelerated (in tile red) subtypes in modern humans. Boxplots represent median and interquartile range. (F) Gene ontology enrichment of HS-Open-CREs with modern human-specific variants in L2-3_2 subtype. (G) snATAC-seq coverage plots of Human-Open-CREs near EPHB1 in L2-3_2 neurons which has a modern human variant. Track scales are the same in all species.

Figure 5:. Accessibility changes highlight subtype-specific transcription factor target evolution in the human brain .

(A-B) FOS / JUN (A) and FOX (B) family transcription factors (TFs) enrichments in HS-Open-CREs. Heatmap shows the log fold changes of TF motif enrichment within HS-Open-CREs per TF motif and per subtype (Blue asterisks: FDR < 0.05. Red rectangles: distinctly more accessible TFs within these subtypes). (C) Ratio of HS substitutions by CRE length per group of CREs. Groups from left to right: HS-Open-CREs that contain at least one motif, HS-Open-CREs that do not contain a motif, NS-CREs (non-significant CREs) that contain a motif, NS-CREs that do not contain a motif. Only the enriched subtypes per TF group and the highlighted TFs -in panel A- per TF family were used for these comparisons. (N = 1519, 2894, 38486, 109452 CREs left to right). Boxplots represent median and interquartile range. P-value: Wilcoxon rank sum test (two-sided). (D) Same CREs as (C) except using the mean log likelihood ratios per CREs computed in the HAR analysis.(E-F) Same comparison as (C-D) except using FOX motifs. (N = 1338, 3075, 38774, 109164 CREs left to right.). (G) Track plot of an HS-Open-CRE with a HAR and contains a human-specific gain of a FOS/JUN motif. Identical sequences with respect to the human sequence are shown with dots. The motif representation on the bottom is a FOS motif.