Single-cell multi-cohort dissection of the schizophrenia transcriptome

Abstract

The complexity and heterogeneity of schizophrenia have hindered mechanistic elucidation and the development of more effective therapies. Here, we performed single-cell dissection of schizophrenia-associated transcriptomic changes in the human prefrontal cortex across 140 individuals in two independent cohorts. Excitatory neurons were the most affected cell group, with transcriptional changes converging on neurodevelopment and synapse-related molecular pathways. Transcriptional alterations included known genetic risk factors, suggesting convergence of rare and common genomic variants on neuronal population-specific alterations in schizophrenia. Based on the magnitude of schizophrenia-associated transcriptional change, we identified two populations of individuals with schizophrenia marked by expression of specific excitatory and inhibitory neuronal cell states. This single-cell atlas links transcriptomic changes to etiological genetic risk factors, contextualizing established knowledge within the human cortical cytoarchitecture and facilitating mechanistic understanding of schizophrenia pathophysiology and heterogeneity.

Fig. 1.. Multiresolution dissection of cellular subpopulations.

(A) Overview of study design and data analysis strategies. Nuclei were pooled after sample-specific hashtag labeling, allowing removal of intersample doublets. Cell annotation was performed on the combined dataset. DE analysis was performed within each dataset independently and results merged through meta-analysis. Downstream analyses of biological pathways and cis- and trans-regulatory factors related to cell type–specific schizophrenia DEGs (right). CON, control; SZ, schizophrenia; McL, McLean cohort; MSSM, Mount Sinai School of Medicine cohort. (B) ACTIONet plot of putative cell types. Green and red clusters represent excitatory and inhibitory subtypes of neurons, respectively, with darker shades indicating an association with deeper cortical layers. (C) Annotation of cell types using curated markers (31) from previous snRNA-seq (26, 27, 79) and spatial transcriptomics (38) studies. (D) Projection of known marker genes verifies cell-type annotations and cortical layer associations (85).

Fig. 2.. Cell type–specific DEGs and associated biological pathways in schizophrenia.

(A) Number of down-regulated (blue, left of 0) and up-regulated (red, right of 0) genes with FDR < 0.05 and abs(log₂(fold change)) > 0.1 in meta-analysis combining results of both cohorts. Shaded regions of each bar indicate the number of genes also significantly different with the same direction of change in a bulk tissue RNA-seq study of schizophrenia PFC (20). (B) Significance of overlap between cell type–specific DEG sets identified in each cohort and in the meta-analysis, computed with the R package Super-ExactTest (48). Columns are scaled independently with darker color indicating greater significance of overlap, and cross hatches indicating the specific level of significance of each tested overlap. (C) Concordance of meta-analysis DEGs and DEGs identified by prior bulk cortex RNA-seq (20) with up- and down-regulated DEGs considered separately. Columns are scaled and highlighted as in (B). (D) Signed significance of DE for the 10 most broadly and divergently dysregulated genes. FC, fold change. (E) Biological pathways overrepresented within cell type–specific DEGs determined by gene set enrichment and semantic clustering analysis of gene ontology (GO) terms. A total of 119 enriched GO terms revealed 14 distinct biological themes listed in panel (F) and datafile S7. Numbers indicate significance rank. Themes are named by the most significant term they contain. (F) Aggregate enrichment of GO terms within each biological theme across cell types by down-regulated (blue) and up-regulated (red) DEGs, excluding themes implicated by only one term and nonsignificant cell types. (G) Overrepresentation of synaptic compartment genes [SynGO database (45)] within cell type–specific DEGs.

Fig. 3.. Schizophrenia DEGs implicate a coherently expressed TF module.

(A) Prioritization of TF coexpression modules by overrepresentation of annotated TF targets within DEGs performed using ChEA3 (86). Red bars depict ChEA3 enrichment scores for up-regulated DEGs and blue bars, for down-regulated DEGs. (B) TFs within module 10 form a highly connected protein-protein association network supported by multiple lines of evidence [STRING database (47)]. TFs named in red are fine-mapping prioritized schizophrenia risk genes reported by the Psychiatric Genomics Consortium (12). (C) Genomic loci bound by MEF2C, SATB2, or TCF4 in postmortem human PFC neurons are significantly associated with schizophrenia DEGs across multiple neuronal and no non-neuronal populations.

Fig. 4.. Association of DEGs with genetic risk variants across cell types and disorders.

(A) Enrichment of GWAS signal (gene-level scores computed with H-MAGMA; standard MAGMA results presented for comparison in fig. S16) within cell type–specific schizophrenia DEGs (red) and equal numbers of genes preferentially expressed in each cell type (PEG, blue) assessed by permutation test against randomly selected genes. Asterisks indicate significant deviation from random expectation (Bonferroni-corrected P < 0.01). PEGs do not significantly deviate from random expectation in any cell type. (B) Genome-wide association between cell type–specific schizophrenia differential expression and genomic variants implicated by GWAS across multiple neuropsychiatric disorders computed with H-MAGMA (51), and by exome sequencing for schizophrenia (13) (SZ exome) computed by gene set enrichment analysis (GSEA). Pie charts and fill color both indicate −log₁₀(FDR). Empty tiles indicate no significant association (FDR > 0.01). (C) Visualization of genes most strongly associated (FDR < 0.05) with common schizophrenia risk variants (GWAS) computed with H-MAGMA (y axis, left), or rare schizophrenia risk variants (exome sequencing, y axis, right), with significant expression changes in schizophrenia (x axis) within excitatory neuronal populations showing the strongest associations in panel (B).

Fig. 5.. Cell type–specific differential expression of high-confidence schizophrenia risk genes.

Cell type–specific and bulk PFC differential expression of 114 genes within the Psychiatric Genomics Consortium’s broad fine-mapped set (12). Shading indicates significance and direction of change (red, up-regulated; blue, down-regulated). Significance of association of each gene with schizophrenia common risk loci computed with H-MAGMA (Gene Link Score). Significance of the index SNP linked to each gene by the PGC3’s statistical fine mapping (SNP Score). PGC prioritized indicates that the gene is present within the Psychiatric Genomics Consortium’s set of 120 prioritized schizophrenia risk genes. Index SNPs in red are linked to multiple DEGs.

Fig. 6.. Heterogeneity of transcriptional changes across individual subjects.

(A) Individuals plotted by TPS computed within the space of schizophrenia-associated transcriptional change observed in the McL cohort (x axis) or the MSSM cohort (y axis). Values are scaled, with larger values indicating a greater association with schizophrenia. Four subgroups of individuals are identified: SZ, red points in the upper right quadrant; CON, blue points in the lower left quadrant; SZ_CON-like, red points outlined in blue in the lower left quadrant; and CON_SZ-like, blue points outlined in red in the upper right quadrant. Off-diagonal individuals (gray points) are inconsistently classified between analyses. R, Pearson correlation coefficient. Regression line computed using all data, including gray points. (B) TPS within individuals clustered into the four groups identified in panel (A). Values are scaled, with darker red indicating a value more characteristic of the schizophrenia group, and darker blue indicating a value more characteristic of the control group. Correlation, cell type–specific correlation between each individual’s transcriptional signature and the transcriptional signature characteristic of the SZ group; TPS, Transcriptional Pathology Score averaged across neuronal cell types for each individual; Ex_SZCS, In_SZCS, Ex_SZTR, scaled expression of each transcriptional signature within each individual; PRS, schizophrenia polygenic risk score. Individuals with inconsistent classification in panel (A) were omitted from panel (B). (C) Enrichment of the In_SZCS transcriptional signature (top) within all individuals across the four groups identified in panel (A), showing a high association with diagnostic subgroups. PRS (bottom) does not show this association. (D) Correlation between cell-state transcriptional signatures and schizophrenia transcriptional change across individuals, with In_SZCS and Ex_SZCS showing positive (R = 0.56, P = 4.6e-13, and R = 0.35, P = 1.9e-5, respectively) and Ex_SZTR showing inverse association with TPS (R = −0.39, P = 1.8e-6). (E) qPCR validation of increased expression of five top In_SZCS-associated genes in SZ as compared to SZ_CON-like individuals, performed in RNA extracted from whole-cortex tissue samples of 13 SZ and five SZ_CON-like individuals within the McL cohort. Black dots indicate fold change of expression of each target gene in one SZ individual relative to average expression across all SZ_CON-like individuals normalized to 1, computed with the ΔΔct method and beta-2 microglobulin (B2M) used as the reference gene. Red dots indicate outlier values. (F) Enrichment of proteins localized to the synaptic compartment as annotated within the SynGO database (45) within the top 100 genes characterizing the merged Ex_SZCS and In_SZCS transcriptional signatures. 1, postsynaptic specialization; 2, postsynaptic density; 3, integral component of postsynaptic density membrane; 4, postsynaptic membrane; 5, extrinsic component of post-synaptic membrane; 6, integral component of postsynaptic membrane; 7, presynaptic active zone; 8, synaptic vesicle membrane; 9, integral component of synaptic active zone membrane; 10, presynaptic membrane; 11, integral component of presynaptic membrane.