Single-cell dissection of schizophrenia reveals neurodevelopmental-synaptic axis and transcriptional resilience

Multiplexed single-cell profiling of schizophrenia

To investigate schizophrenia-associated cell-type-specific transcriptional disruption within the complex cytoarchitecture of the human brain, we used single-nucleus RNA sequencing (snRNA-seq) to profile postmortem tissue samples from the prefrontal cortex (Brodmann Area 10) (Fig. 1a). We curated a cohort of 24 schizophrenia and 24 control subjects, balanced for gender (12 male and 12 female subjects per group), age (ranging from 22 to 94 years, average 63.5 years), and postmortem interval (ranging from 6.9 to 26.3 hours, average 16.8 hours) (Supplementary Table 1). We reviewed medical records and incorporated psychiatric medication exposures into our analysis to control for potential confounding effects of pharmacotherapy. We multiplexed samples, pooling three schizophrenia and three control samples in each of eight batches, using the Multiplexing Using Lipid Tagged Indices (MULTI-Seq)²⁷. Compared to standard snRNA-seq protocols, sample multiplexing allows us to capture more cells from each individual while reducing batch effects (by profiling both schizophrenia and control samples in the same sequencing library), reducing the rate of undetected doublets (by using sample hashtags to remove cross-individual doublets), and lowering sample preparation costs (by assessing more than one sample per 10x kit channel). After doublet removal, we obtained 560,020 single-nucleus transcriptomes, including 266,431 cells from schizophrenia and 293,589 from control individuals, profiled at an average depth of 12k cells per sample and 35k reads per cell (Supplementary Table 2). We next removed genes expressed in less than ten cells in each batch, and cells with fewer than 500 identified genes or more than 10% of unique molecular identifiers (UMIs) in mitochondrial genes, resulting in 17,460 genes detected in 506,378 cells.

• Download figure
• Open in new tab

Figure 1. Multiresolution dissection of cellular subpopulations.

a. Overview of study design and data analysis strategies. TF - transcription factor. CUT&Tag - Cleavage Under Targets and Tagmentation. b. ACTIONet plot of putative cell types/states. Green and red clusters represent excitatory and inhibitory subtypes of neurons respectively, with darker shades indicating an association with deeper cortical layers. c. Projection of known marker genes verifies cell type annotations and cortical layer associations⁸². d. Annotation of transcriptional archetypes using curated markers from previous studies. e. The percentage of cells within each cell-type/state contributed by SZ and CON subjects. f. Cell-type/state decomposition of individual samples. Colors are consistent with those used in panel b.

Multiresolution dissection of cellular subpopulations

We used our recently developed multiresolution cell-state discovery platform, ACTIONet²⁸, to identify both discrete and continuous cell states (“archetypes”) using a coupled archetypal/network analysis. We first applied ACTIONet to reconstruct the geometry of the cell-state landscape and used it to detect and remove outlier cells that fall on the periphery of the constructed cell-cell similarity graph, keeping only the 386,065 highest-quality cells with reproducible expression patterns. In the second round, we used the ACTIONet on the filtered cells and identified 20 cell states in the combined population of schizophrenia and control subjects, nearly all of which form dense clusters in the cell-cell similarity network (Fig. 1b, Extended Data Fig. 1a), with two notable exceptions that we discuss below. Our 20 cell states captured all major cell types of the human prefrontal cortex, including subtypes of excitatory and inhibitory neurons and glia (Fig. 1b, d). We verified these annotations based on known marker gene expression patterns, by projecting individual marker genes on the cell similarity network (Fig. 1c). De novo markers for these cell types are consistent with our previous studies of the human prefrontal cortex²⁸ (Supplementary Table 3). Neuronal subtypes, particularly excitatory neurons, showed higher numbers of expressed genes and identified UMIs (Extended Data Fig. 2a, b).

We also captured all major subtypes of GABAergic inhibitory neurons, including calcium-binding protein parvalbumin (PV) expressing neurons, neuropeptide somatostatin (SST) neurons, and ionotropic serotonin receptor 5HT3a (5HT3aR) neurons. Within these groups, we detected two PV-expressing subtypes of inhibitory neurons (Basket and Chandelier), and two 5HT3aR-expressing subtypes (VIP+ and Reelin+), and the recently-described Rosehip population²⁹. Cell groupings observed in the ACTIONet plot are consistent with the developmental origin of cardinal interneuron subtypes (medial versus caudal ganglionic eminence), demonstrating their shared transcriptional signatures³⁰.

Excitatory neurons exhibit a strong layer-specific pattern. We observed a consistent association of superficial-layer excitatory neurons (layers II/II) with marker genes CUX2 and CBLN2. However, intermediate and deep-layer excitatory neurons are more intermixed, with layer IV/V neurons enriched for RORB, FOXP2, and RXFP1, while deep layer neurons are marked with TLE4, SEMA3E, and HTR2C genes. Within the deep cortical layers V and VI, we observe three distinct populations of excitatory neurons: cortico-fugal projection neurons (Ex-L5/6) expressing FEZF2, and two distinct populations of cortico-cortical projection neurons (Ex-L5/6-CCa, Ex-L5/6-CCb).

To further investigate differences between Ex-L5/6-CCa and Ex-L5/6-CCb, we focused on cells from control individuals associated with these cell-types and performed differential expression analysis (Wilcoxon’s rank-sum test, Supplementary Table 4). We found that Ex-L5/6-CCa is enriched for genes involved in the dopamine receptor signaling pathway, whereas Ex-L5/6-CCb is enriched for many key genes in glutamate signaling. These distinct expression profiles point to physiological differences between two deep layer neuronal populations likely related to specifics of their afferent and efferent connectivity with distinct cortical and subcortical regions.

In addition to major cell types and their corresponding subtypes, we found two transcriptional archetypes that capture continuous cell states shared across multiple subtypes of excitatory neurons. The first cross-cutting cell state (Ex-NRGN) captures both a localized cell-type characterized by the expression of the NRGN gene described in prior snRNA-seq studies of postmortem human brain²³, as well as a transcriptional signature distributed among cells of different subtypes. The Ex-NRGN cell-state does not show enrichment for layer-specific markers. In addition to NRGN, this population is marked by expression of the Brain Expressed and X-Linked gene family members (BEX1 and BEX3), calcium-binding enzymatic cofactor Calmodulin 3 (CALM3), and YWHAH, which encodes a 14-3-3 signal transduction protein implicated in the conversion to the psychosis of at-risk individuals³¹. Cells associated with the Ex-NRGN cell state highly express mitochondrial genes (Extended Data Fig. 2c), an observation reproducible in an independent study²³ (Extended Data Fig. 2d).

The second cross-cutting cell-state (Ex-SZTR) is of particular interest, as it is preferentially associated with excitatory neurons from schizophrenia patients (Fig. 1e). Ex-SZTR cells are enriched in superficial layer excitatory neurons, as evidenced by the expression of layer-specific marker genes as well as the localization of these cells within the cell network (Fig. 1b,d, Extended Data Fig. 1a). The top-ranked associated genes are enriched for pathways related to synapse organization and synaptic plasticity, as well as neuronal process morphology (Supplementary Table 5). We name this cell state Ex-SZTR for “schizophrenia transcriptional resilience” as it is preferentially found in schizophrenia individuals whose transcriptional profiles are surprisingly non-schizophrenia-like, as we discuss later in the text.

Cell-type-specific schizophrenia-associated transcriptional perturbations

Across all cell types/states, we identified a total of 1,637 up-regulated and 2,492 down-regulated differentially expressed genes (DEGs) in schizophrenia (3,742 distinct schizophrenia-perturbed genes, as 387 genes are both up- and down-regulated in distinct cell types), with the highest number of up-regulated genes in Ex-SZTR and the highest number of downregulated genes observed in Ex-L5/6-CCb (Fig. 2a, Supplementary Table 6). The majority of observed perturbations occur in neuronal populations, and ∼60% of observed perturbations are downregulation events. Among inhibitory neurons, PV-expressing subtypes (both Chandelier and Basket) which have been the focus of much prior work³², indeed show the highest number of perturbed genes. Pairwise comparison of the schizophrenia-associated transcriptional dysregulation between different cell subpopulations reveals similarity in transcriptional changes within all neuronal subpopulations. Among neuronal cell types, we observed a particularly high concordance of dysregulated genes among superficial layer excitatory neurons (Ex-L2/3, Ex-L4, and Ex-L4/5) and Ex-L5/6-CCb (Extended Data Fig. 3).

• Download figure
• Open in new tab

Figure 2. Differential gene expression analysis.

a. Total count of up- and downregulated genes in each cell-type/state and the overlap of cell-type/state-specific SZ DEGs with DEGs observed in previous studies of bulk cortical tissue²⁸.b. Cell-type-specific differential expression of selected top-ranked genes across all 20 cell-types/states as well as previous studies of bulk cortical tissue. Many selected genes demonstrate consistent up or downregulation across most cell populations, while others show distinct patterns across major categories (NRXN3, RASGRF2, BACH2, DMD), and several genes are dysregulated in opposite directions within Ex-SZTR and multiple neuronal cell-types (GRIA4, DLC1, KIAA0232, HDAC4, ATF7IP2). Blue indicates downregulation, and red indicates overexpression in SZ. Comparisons with nominal p>0.05 are not colored.c. Ranking of individuals based on an aggregate Transcriptional Pathology Score computed across all neuronal cell-types/states. Red indicates a transcriptomic signature more typical of SZ, while blue indicates a signature more typical of CON. *Individual CON7 has a first degree relative with SZ. d The first column depicts the relative pseudobulk expression across all neuronal cell-types of TCF4, CLU, SHANK2, or UNC13A on the y-axis plotted against the transcriptional pathology score ranking of each subject on the x-axis. Columns two through four show representative photomicrographs of fluorescent in situ hybridization for CUX2 in red, a marker of excitatory neurons in layers II and III of the cortex, and TCF4, CLU, SHANK2, or UNC13A in blue or green. Column five depicts quantitation of the in situ hybridization signal as detected transcripts per single CUX2 positive cell using DotDotDot⁷⁶. Points indicate individual cell counts, boxes indicate the median and the interquartile range, and whiskers indicate the largest value within 1.5 times the interquartile range above the 75th percentile. e. Functional enrichment of cell-type/state-specific SZ DEGs within neuronally relevant biological processes.

To evaluate the reproducibility of our cell-type-specific dysregulated genes, we compared them to bulk-level dysregulated genes across 559 schizophrenia and 936 control homogenized prefrontal cortex samples¹⁹. We observed a significant overlap between the 2,450 upregulated (p-value<10⁻⁸, Fisher’s exact test) and 2,371 downregulated (p-value<10⁻²⁹) bulk-level genes and the upregulated and downregulated gene sets identified here, with downregulated genes in excitatory neurons having the strongest association (Fig. 2a). However, among our cell-type-specific DEGs we identify novel genes not detected in prior tissue-level observations, including genes in the Ex-SZTR cell state and subtypes of inhibitory neurons, highlighting the importance of single-cell resolution analysis.

Among our DEGs (Fig. 2b), the transcriptional regulator TCF4¹⁴, which lies in a schizophrenia-associated genetic locus, was the most widely-perturbed gene, upregulated in 14 of 20 cell types. Indeed, TCF4 was also detected as upregulated in bulk tissue RNA-Seq¹⁸. TCF4 plays critical roles in both neurodevelopment and in regulating the excitability of prefrontal cortex neurons³³, and was previously predicted as a “master regulator” of schizophrenia gene network perturbations³⁴. Additionally, the molecular chaperone CLU, which also lies in a schizophrenia-associated genetic locus¹⁴ and is most prominently expressed in astrocytes, was broadly overexpressed across astrocytes, all excitatory neuron cell types, and most inhibitory neuron cell types. Confirming our findings, CLU was previously recognized to be hypomethylated in schizophrenia bulk postmortem brain samples³⁵, and overexpressed in laser-capture-microdissected excitatory neurons and in parvalbumin-expressing inhibitory neurons^36,37. We also observed two gene families, neurexins (NRXN1, NRXN2, NRXN3) and SHANK genes, to be commonly dysregulated across multiple neuronal subpopulations. Neurexin genes serve as presynaptic cell-adhesion molecules and interact with neuroligins within the postsynaptic membrane to mediate multiple aspects of synapse formation, structure, and function³⁸. In contrast, SHANK genes encode postsynaptic scaffolding proteins essential for the organization of glutamatergic synapses³⁹.

Multi-gene transcriptional pathology score

We next calculated a “transcriptional pathology score” (TPS) for each individual, based on the correlation of their expression deviation (relative to the mean of all individuals) in each cell type, with the vector of schizophrenia-associated differential expression in that cell type (Fig. 2c, Extended Data Fig. 4, Supplementary Table 7). As expected, schizophrenia individuals showed significantly higher transcriptional pathology scores than control individuals lying at opposite extremes of the TPS distribution, indicating agreement between our predicted dysregulation scores and known phenotypic states. Individual CON7 was an exception to this trend, showing schizophrenia-like expression profiles but lacking a schizophrenia diagnosis; however, CON7 had a first-degree relative (his son) diagnosed with schizophrenia, suggesting a higher genetic predisposition than expected possibly driven by a familial strong-effect genetic alteration. Surprisingly, we found that TPS ranking was consistent across all neuronal cell types, rather than confined to the deep layer excitatory and PV-expressing basket cell-types that showed the highest number of DEGs, indicating that global pan-transcriptome gene-expression signatures of schizophrenia are robust across all neuronal cell-types (glial cell types did not show such consistency).

Strikingly, schizophrenic individuals with an abundant subset of excitatory neurons in the Ex-SZTR cell state showed among the lowest transcriptional pathology scores, even lower than most control individuals. In fact, the two control individuals marked by the Ex-SZTR state show the two lowest schizophrenia transcriptional pathology scores. This suggests that the presence of the Ex-SZTR cell state is correlated with “transcriptional resilience”, whereby even schizophrenia cases show transcriptional profiles consistent with control individuals, and control individuals show even more non-schizophrenia-like transcriptional signatures. We do not find the Ex-SZTR cell-state to be associated with exposure to specific medications or levels of exposure to pharmacologic classes. Instead, we interpret it as a potential mechanism of transcriptional compensation or resilience to schizophrenia at the molecular level worthy of further investigation.

Experimental validation of differentially expressed genes

We next used fluorescence in situ hybridization (RNAscope) to experimentally validate both the cell-type-specific expression and the differential expression of four differentially-expressed genes across six individuals. We selected two upregulated genes (TCF4, CLU) and two down-regulated genes (SHANK2, UNC13A), which showed strong and consistent changes across the majority of neuronal cell-types, and validated them in two schizophrenia, two control, and two transcriptionally-resilient schizophrenia individuals. The resulting images (Fig. 2d, Extended Data Fig. 5) show clear localization of these transcripts in excitatory neurons of the superficial cortical layers (CUX2 staining). Quantification of transcript abundance (using dotdotdot⁴⁰) confirmed the highest expression for schizophrenia individuals, a lower expression for control individuals, and the lowest expression for transcriptionally-resilient individuals for TCF4 and CLU, and the opposite trend for SHANK2 and UNC13A, strongly confirming the cell-type-specificity, differential expression, and directionality of our findings, and also confirming the surprising behavior of transcriptionally-resilient individuals.

Dysregulated genes converge on synaptic plasticity, organization, and development

Observed gene expression changes suggest that schizophrenia predominantly impacts the transcriptional state of neuronal populations. To investigate whether neuronal processes are consistently dysregulated, or, alternatively, diverse neuronal functions are perturbed across cell subpopulations more specifically, we designed and performed extensive pathway enrichment analyses. We curated and organized relevant biological pathways into 14 neuronally-associated functional categories within three major themes with direct relevance to the etiology of schizophrenia: synaptic organization⁴¹, synaptic plasticity², and neurodevelopment⁴² (Supplementary Table 8a). We observe a pan-neuronal overrepresentation of pathways related to post-synaptic processes within both up- and downregulated genes, consistent with genetic and proteomic evidence that schizophrenia perturbations converge on the postsynaptic density^43,44. Neurodevelopmental pathways show higher enrichment among upregulated genes, whereas glutamate signaling and synaptic plasticity are dominantly downregulated in schizophrenia. Among subpopulations of inhibitory neurons, PV- and SST-expressing neurons show the most robust functional enrichment with up- and downregulated genes, respectively (Fig. 2e, Supplementary Table 8b,c).

In addition to these major categories, we identified the disruption of cytoskeletal processes to be specifically enriched in the Ex-SZTR cell-state. This includes differential expression of genes involved in lamellipodium organization and assembly, migration (CAPZB, CARMIL1, SLIT2), morphologic regulation of axons, dendrites, and dendritic spines (CARMIL1, GOLPH3, SRGAP2, VAV2⁴⁵, VAV3, WASF3⁴⁶), and axon guidance (ABLIM1, NCK1⁴⁷, PTPRO, SLIT2). These findings in the Ex-SZTR cell-state, which is predominantly associated with markers of the superficial cortical layers, are consistent with previous studies describing aberrant pyramidal cell density and morphology predominantly within layers II and III^48–50.

Correlated patterns of genetic and cell-type-specific expression perturbations

Seeking insights into the association between transcriptional alterations and schizophrenia manifestation, we investigated the relationship between dysregulated genes and genetic risk factors. By assessing the cell-type-specific differential expression of genes within 145 genomic loci significantly associated with schizophrenia¹⁴, we discovered significant differential expression events in at least one cell type within 68 of these loci (Fig. 3a, Supplementary Table 9). These novel associations identify putative causal mechanisms for more schizophrenia risk loci than previously possible using bulk tissue gene expression data, explaining the functional relevance of these loci with candidate genes, cell-type of action, and direction of effect.

• Download figure
• Open in new tab

Figure 3. Correlation of differentially expressed genes with GWAS-associated genesets.

a. Schizophrenia DEGs suggest mechanisms of action for 68 of 145 GWAS implicated loci. Shown are the 68 GWAS loci¹⁴ containing significant SZ DEGs, along with the most perturbed gene, the cell-type in which the largest DE event occurs, and a heatmap depicting all DE events for that gene within neuronal populations (red - upregulated; blue - downregulated, circles mark the maximum differential expression event for each gene). b. Correlation plots depicting the overlap between cell-type/state-specific schizophrenia DEGs and genes implicated by Genome-Wide Association Studies of six distinct neuropsychiatric disorders. Within each pairwise association, the ratio of purple clock-face indicates the relative amount of overlap, with the intensity of coloring indicating the significance of the association. c. Visualization of the relationship between the level of significance of genome-wide association with schizophrenia (y-axis) and the magnitude of differential expression in schizophrenia (x-axis) for individual genes within three neuronal populations.

The majority of schizophrenia DEGs providing explanation for these 68 genomic loci were maximally perturbed in excitatory neuronal populations, while 13 were most robustly altered in inhibitory neurons, and two genes in astrocytes and two in oligodendrocytes. Perturbation of explanatory genes within these regions was split nearly equally between up and downregulation (32 up vs. 36 down). Many explanatory genes were dysregulated in a small number or only a single cell-type, or dysregulated in opposing directions across multiple cell-types, resulting in their not being identified as schizophrenia DEGs in prior studies of bulk tissue. For example, the CALB2 gene encodes the calcium binding protein calretinin which is expressed specifically in In-VIP and In-SST populations and contributes to long-term potentiation through regulation of neuronal excitability⁵¹, and CALB2 was specifically upregulated in only the In-VIP cell-type. ALMS1 encodes a microtubule organizing protein, and this gene was downregulated specifically in PV-expressing interneurons. Similarly, we found multiple GWAS explanatory genes that were specifically dysregulated in the Ex-SZTR cell-type (CPEB1, GPR135, ZNF804A), or upregulated in Ex-SZTR but downregulated across other populations (BCL11B, RALGAPA2).

In addition to schizophrenia¹⁴, we further considered four psychiatric disorders that are known to share genetic risk factors⁵²: major depressive disorder (MDD)⁵³, bipolar disorder (BD)⁵⁴, autism spectrum disorder (ASD)⁵⁵, and attention-deficit/hyperactivity disorder (ADHD)⁵⁶. This approach allows us to distinguish schizophrenia from general psychiatric illness-related associations. As a point of contrast, we included Alzheimer’s disease (AD)⁵⁷, a neurodegenerative disorder not expected to be genetically related to schizophrenia. For each gene/trait pair, we computed an aggregate genetic perturbation score using H-MAGMA⁵², a tool that associates GWAS risk variants to genes based on proximity to gene bodies, promoter regions, or regions linked to genes through chromatin looping. We only considered evidence of distal regulatory links occurring in the adult dorsolateral prefrontal cortex²⁶. As a measure of transcriptional perturbation for schizophrenia, we used the absolute value of the t-statistics computed in the cell-type-specific differential expression analysis.

Correlation analysis between genetic and transcriptional perturbation scores for schizophrenia uncovered a strong association within neuronal subpopulations, suggesting a substantial causal effect for observed transcriptional alterations (Fig. 3b). We found that many of the cell-types whose transcriptional perturbations show a strong association with schizophrenia genetic risk variants are also associated with risk variants for bipolar and major depressive disorders, which is consistent with their previously-reported strong genetic relationships both at the level of genetic correlations and gene-level overlaps⁵². We found weaker relationships between schizophrenia transcriptional perturbations and ASD risk^52,58, consistent with the known lower correlation of genetic risk between these disorders⁵⁸. Indeed, schizophrenia transcriptional perturbations showed overlap with genetic risk loci for all neuropsychiatric disorders we assessed, while as expected we found no overlap with genetic risk for the neurodegenerative disorder AD.

Across excitatory neuronal subpopulations, the relation between transcriptional and genetic perturbation was strongest for the Ex-SZTR cell state, followed by deep-layer corticocortical projection neurons. The Ex-SZTR cell state is predominantly associated with superficial-layer excitatory neurons (layers II/III), whereas cortico-cortical excitatory neurons are enriched for layer V. The observed layer-specificity of perturbations is consistent with the enrichment of schizophrenia genetic variants proximal to genes preferentially expressed in layer II and layer V⁵⁹. Among inhibitory neurons, we detected the strongest association between schizophrenia DEGs and genomic variants in In-PV-Basket cells, consistent with the known dysregulation of PV-expressing GABAergic neurons in schizophrenia hypothesized to contribute to the disruption of synchronous neural activity in schizophrenia^60,61. Unlike neuronal populations, we did not observe any significant association between patterns of differential expression and genetic association in either glial or endothelial cells, suggesting non-neuronal cell types do not play a primary role in schizophrenia pathologies mediated by transcriptional dysregulation.

We next examined genes that show both genetic and transcriptional perturbations, focusing on the cell-types most enriched for schizophrenia DEGs (Ex-SZTR, deep-layer cortico-cortical excitatory neurons, and In-PV-Basket cells, Fig. 3c, Extended Data Fig. 6). Among these genes we found TCF4, a regulator of pyramidal cell excitability whose target genes are also schizophrenia-associated and cluster in neurodevelopmental pathways⁶²; CLIP1, a linker protein involved in microtubule trafficking of cargo within axons and dendrites⁶³; CACNA1C, an ion channel with key roles in synaptic plasticity and neurodevelopment; ZKSCAN3, a transcription factor implicated in autophagy; and PGBD1, a transposase with brain-specific expression and unknown functionality.

Transcriptional dissection of trans-acting factors reveals common regulators

While individual DEGs may be affected by cis-acting genetic variants, coordinated expression changes are often driven by common upstream transcriptional regulators⁶⁴. To prioritize such potential factors, we ranked a total of 1,632 transcription factors (TFs), based on the degree to which their putative target genes overlap with schizophrenia DEGs⁶⁵ (Supplementary Table 10). From this data, we identified, for each cell type, a list of candidate regulators and tested whether their functional annotations converge to similar pathways as those observed for DEGs (Fig. 2c). Surprisingly, TFs were strongly enriched only in neurodevelopmental pathways (Fig. 4a, Supplementary Table 11), and not in synaptic function or signaling pathways.

• Download figure
• Open in new tab

Figure 4. Schizophrenia-associated transcriptional regulators.

a. Functional enrichment of the top-ranked transcriptional regulators identified by ChEA3 analysis. b. Co-expression heatmap of the most perturbed neurodevelopmentally-driven TFs across pseudobulk profiles from all cell-types/states and all individuals identifies five distinct TF modules, with high enrichment of GWAS-associated TFs in the MEF2C/SATB2 module (boxed top left). c. Overlap of TFs within the MEF2C/SATB2 module and the ME37 module identified in the PyschENCODE brain development dataset.

To assess the association of schizophrenia genetic risk loci with transcriptional regulators in addition to downstream transcriptional dysregulation directly, we focused on a subset of GWAS-associated TFs targeting DEGs. We identified a set of seven such factors (TCF4, SATB2, MEF2C, FOXG1, SOX5, ZNF536, and ZNF804A) with very strong overlap between their putative target genes and schizophrenia DEGs across all neuronal cell-types/states (Extended Data Fig. 7). Interestingly, these TFs again point to neurodevelopmental functions. TCF4 is a broadly expressed helix-loop-helix transcription factor involved in nervous system development that, when disrupted, causes the severe neurodevelopmental disorder Pitt Hopkins Syndrome; SOX5 is implicated in fate determination and regulation of corticofugal projection neuron development, with loss-of-function perturbations leading to disrupted neuronal proportions and emergence timing⁶⁶; ZNF804A is involved in neurodevelopment, neuronal migration, and protein translation. To test whether additional neurodevelopmentally-associated transcriptional regulators not implicated by schizophrenia genetic studies might be functionally relevant in the context of schizophrenia by sharing converging behavior with these lead factors, we used the set of top-ranked TFs that are also annotated with neurodevelopmental functions to identify modules of coherent TFs. We constructed a TF-TF co-expression network and applied a clustering algorithm to identify five neurodevelopmentally-related, schizophrenia-associated TF modules (Fig. 4b, Supplementary Table 12, 13). In addition to the lead GWAS-associated factors, we found additional factors with known neurodevelopmental roles, including TBR1, which interacts with a SOX5/SATB2 circuit thought to regulate layer specification during fetal neurogenesis⁶⁷.

Among these modules, we found a TF group highly enriched for both schizophrenia- and neurodevelopmental delay (NDD)-associated genetic variants (schizophrenia p-value:1.2×10⁻⁷, NDD p-value: 3×10⁻⁶, Fisher’s exact test). NDD encompasses a broad category of disorders characterized by disrupted development of the nervous system leading to abnormal brain function, including intellectual disability, movement, speech, and tic disorders⁶⁸. We observed a significant overlap between genetically-identified NDD-associated TFs (de novo mutations and CNVs)⁶⁹ and our top-ranked schizophrenia TFs, in particular, TFs mediating upregulated genes in Ex-SZTR and Ex-NRGN, as well as downregulated genes across multiple subtypes of excitatory neurons (Extended Data Fig. 8). We observed a significant over-representation of the schizophrenia-GWAS associated TFs among the top-10 ranked TFs (lead NDD TFs), where NDD-associated factors are sorted based on their overall ChEA score (p-value<7×10⁻⁷, Fisher’s exact test). The majority of the lead NDD TFs (7/10) are also involved in neurodevelopment (Extended Data Fig. 9). Thus, we identified a group of coregulated TFs that are associated with (1) downstream schizophrenia-associated transcriptional dysregulation in neurons of the adult cortex, (2) schizophrenia genetic risk, and (3) genetic risk for disrupted neurodevelopment.

Two of the modules identified herein are consistent with neurodevelopmentally-associated gene modules recently reported by the PsychENCODE Consortium (PEC)⁷⁰. GWAS-enriched modules MEF2C/SATB2 and ZNF804A/ZEB2 strongly map to modules ME37 and ME32 in the PEC dataset, both of which were also reported to be enriched for schizophrenia variants. The latter module contains 19 TFs, 9 of which overlap with the MEF2C/SATB2 module (adjusted p-value <10⁻¹⁰, Fisher’s exact test). The overlap covers the majority of GWAS-associated TFs, including SOX5, SATB2, MEF2C, TCF4, and EMX1 (Fig. 4c).

Experimentally validated TF targets confirm developmental-synaptic axis

To further investigate the link between the prioritized TFs and schizophrenia DEGs, we next tested whether active cis-regulatory elements that are targeted by these TFs show preferential association with schizophrenia DEGs. To this end, we performed Cleavage Under Targets and Tagmentation (CUT&Tag⁷¹) assays to map the binding of MEF2C, SATB2, SOX5, and TCF4 genome-wide in neuronal nuclei sampled from four schizophrenia and four control individuals within our larger cohort. For each transcription factor, we defined a master-set of targeted regions that are identified in at least one sample. These master-sets of regulatory elements were then projected onto the set of regulatory elements defined by the PEC developing human brain dataset⁷⁰, which is a union of enhancers/promoters across different stages of brain development (Extended Data Fig. 10). Finally, we mapped these TF-centric active regulatory elements to their putative target genes, considering both proximal and distal interactions.

Given the set of target genes for each of these factors, we performed functional enrichment analysis to uncover the neuronally-associated categories that are related to each TF’s targets (Fig. 5a). Observed targets of these TFs show a functional enrichment profile highly similar to that of up- and downregulated schizophrenia DEGs (Fig. 2e), with prominent enrichment in categories related to postsynaptic and neurodevelopmental processes, supporting a functional role for these TFs in dysregulation of these events in schizophrenia.

• Download figure
• Open in new tab

Figure 5. Functional characterization of TFs targeting schizophrenia DEGs.

a. The overlap between neuronal cell-type/state-specific schizophrenia DEGs and bindings sites for MEF2C, SATB2, SOX5, and TCF4 directly assessed by CUT&Tag assays in FANS sorted neuronal nuclei from four schizophrenia and four control individuals. b. TF targeted genes enrich neuron-associated functional categories with a pattern highly similar to that seen for schizophrenia DEGs (Fig. 2e). c. Association of TFs targeting up- (red) and downregulated (blue) schizophrenia DEGs with developmental stage-specific enhancers within each neuronal cell-type/state. In all panels significance of enrichments was assessed with Fisher’s exact test.

We then assessed the overlap of target genes with up- and downregulated schizophrenia DEGs independently (Fig. 5b). We found a high degree of overlap between identified target genes and schizophrenia DEGs across a wide-range of excitatory neurons, with the largest overlap between TF targeted genes and Ex-L4/5 upregulated and EX-L5/6-CCb downregulated DEGs. Among inhibitory neuronal populations, PV and SST expression subtypes which originate in the medial ganglionic eminence showed markedly stronger associations than inhibitory neurons originating in the caudal ganglionic eminence (In-Reelin, In-Rosehip, and In-VIP). Finally, we did not observe a strong association between identified target genes and DEGs in non-neuronal cell-types.

Target genes for SOX5 and SATB2, which play prominent roles in early neurodevelopment, exhibit a higher degree of overlap with upregulated genes than MEF2C or TCF4. Conversely, SOX5 has an insignificant association with downregulated genes, whereas MEF2C shows the strongest association with downregulated genes, consistent with MEF2C’s role as a transcriptional repressor⁷². This suggestion of up and downregulated schizophrenia DEGs having distinct associations with neurodevelopmental stages prompted us to investigate how the regulatory landscape of schizophrenia DEGs varies across brain development. We identified cis-regulatory elements associated with different stages of brain development based on H3K27ac-enriched genomic regions profiled in the human brain by the PEC⁷⁰, and linked enhancers with genes using evidence of physical chromatin interaction (HiC profiles) and genomic proximity. Associating fetal and adult regulatory elements independently with schizophrenia DEGs, we found that upregulated genes are predominantly linked to enhancers that are specific to the fetal stage, whereas downregulated genes are linked to enhancers that are either shared between the fetal and adult brain or are unique to the adult brain (Fig. 5c). Among the strongest connections, we found Ex-SZTR upregulated and Ex-L5/6-CCb downregulated genes to be associated with shared enhancers, while Ex-L4/5 upregulated genes are enriched for fetal-specific enhancers.

Assembly of the Tissue Cohort

We obtained postmortem human Brodmann Area 10 tissue from 24 schizophrenia subjects and 24 control individuals matched for age, gender, and postmortem interval from the Harvard Brain Tissue Resource Center at McLean Hospital. Institutional review board approval was obtained by the Harvard Brain Tissue Resource Center for the collection, storage, and distribution of brain tissue and de-identified clinical information for each case. Each case was assigned a consensus diagnosis by two psychiatrists based on a review of medical records and a questionnaire completed by the donor’s family. All cases were examined histologically, and cases with neuropathology diagnosed by tissue examination were excluded. Cases were obtained through family referral, and no cases were referred by a medical examiner’s office. Demographic variables for the assembled cohort are listed in Table S1. Upon arrival at the Harvard Brain Tissue Resource Center, fresh brains were dissected, and Brodmann Area 10 tissue was identified and quickly frozen with liquid nitrogen vapor prior to storage at −80°C.

Assessment of Medication Exposures

To control for the effects of psychotropic medications on disease-associated changes in gene expression, medical records were reviewed and each subject’s history of medication exposure was compiled for inclusion in our differential gene expression analysis. Medications were clustered into categories including selective serotonin reuptake inhibitors, tricyclic antidepressants, mood stabilizers, antiepileptic drugs, benzodiazepines, and neuroleptics. Neuroleptic drugs were further classified by structural category, including benzisoxazoles, butyrophenones, dibenzodiazepines, indoles, phenothiazines, phenylpiperazines, thienobenzodiazepines, and thiothixenes. Within each medication category, the total number of distinct agents prescribed to each subject at any time was input to the linear model (Supplementary Table 1).

Sample processing and single-cell RNA sequencing

A Nissle stained cryosection from each tissue block was examined microscopically to verify that sampled tissue included all six cortical layers and underlying white matter. On dry ice 50 mg of tissue was cut from the original block and stored at −80°C until further use. Tissue samples were then processed in batches of nine (three schizophrenia, three control, and three samples not analyzed in the current study) using a protocol adapted from a previous study²². Tissue was thawed in 0.5 ml ice-cold homogenization buffer (320 mM sucrose, 5 mM CaCl2, 3 mM Mg(CH3COO)2, 10 mM Tris HCl pH 7.8, 0.1 mM EDTA pH 8.0, 0.1% IGEPAL CA-630, 1 mM β-mercaptoethanol, and 0.4 U μl−1 recombinant RNase inhibitor (Clontech)) prior to homogenization with 12 strokes in a 2 ml Wheaton Dounce tissue grinder. Tissue homogenates were passed through a 40 µm cell strainer, mixed with an equal volume of working solution (50% OptiPrep density gradient medium (Sigma-Aldrich), 5 mM CaCl2, 3 mM Mg(CH3COO)2, 10 mM Tris HCl pH 7.8, 0.1 mM EDTA pH 8.0, and 1 mM β-mercaptoethanol), layered on top of an Optiprep density gradient (750 µl 30% Optiprep solution (30% OptiPrep density gradient medium,134 mM sucrose, 5 mM CaCl2, 3 mM Mg(CH3COO)2, 10 mM Tris HCl pH 7.8, 0.1 mM EDTA pH 8.0, 1 mM β-mercaptoethanol, 0.04% IGEPAL CA-630, and 0.17 U μl−1 recombinant RNase inhibitor) on top of 300 µl 40% optiprep solution (35% OptiPrep density gradient medium, 96 mM sucrose, 5 mM CaCl2, 3 mM Mg(CH3COO)2, 10 mM Tris HCl pH 7.8, 0.1 mM EDTA pH 8.0, 1 mM β-mercaptoethanol, 0.03% IGEPAL CA-630, and 0.12 U μl−1 recombinant RNase inhibitor)) and centrifuged at 10,000 g for 10 minutes at 4°C. Nuclei were recovered from the gradient, and resuspended in an equal volume of 1% BSA in phosphate buffered saline (PBS) prior to labeling with sample-specific cholesterol conjugated oligonucleotide hashtags (Integrated DNA Technologies). Labeled nuclei were washed with 1% BSA in PBS and pelleted by centrifugation at 500 g at 4°C for 5 minutes for three washes. After the final wash, nuclei were counted on a hemocytometer, and equal numbers from each sample were combined. The pooled nuclei were then applied to all eight channels of one 10x Genomics Chip B, targeting recovery of 20,000 nuclei per channel. 10x Genomics and hashtag libraries were prepared as per standard 10x Genomics Chromium Single Cell 3’ Reagent Kits v.3 and MULTI-Seq²⁷ protocols.

Libraries were sequenced in batches of two on an Illumina NextSeq500 instrument (two snRNA-Seq libraries and two hashtag libraries per flowcell - average 3.6E7 reads per hashtag library), and an additional round of sequencing was performed for all snRNA-Seq libraries on an Illumina NovaSeq instrument (eight or 16 libraries per NovaSeq S2 flowcell) to achieve an average sequencing depth of 35,142 reads per cell. Gene count matrices were generated by aligning reads (including intronic reads) to the hg38 genome using 10x Genomics Cell Ranger software v3.0.1 (Supplementary Table 2).

Cross-individual doublet-detection using sample hashtags

The deMULTIplex R package²⁷ was used to process hashtag FASTQ files, extracting 10x barcode, sample hashtag, and UMI information for each read. Duplicated UMI and mismatched hashtag reads were excluded and retained reads were converted to a 10x barcode by sample hashtag count matrix. This count matrix was processed with the Seurat R package⁷⁶ using the HTODemux function to cluster cells based on sample hashtag counts and determine a count threshold for each hashtag based on a negative binomial distribution applied to the cluster with the lowest expression for that hashtag. This threshold identified each cell as positive or negative for each sample hashtag, and cells identified as positive for more than one hashtag were assigned as inter-sample doublets and removed from the study (Supplementary Table 2).

Characterization, visualization, and annotation of cell subpopulations using ACTIONet

Filtered data after doublet-detection was used as input to the archetypal analysis for cell-type identification (ACTION) algorithm ⁷⁷ to identify a set of transcriptional archetypes, each representing a potential underlying cell type/state. Using ACTION-decompositions with varying numbers of archetypes, we employed our recently developed ACTION-based network (ACTIONet) framework²⁸ to create a multi-resolution nearest neighbor graph. A modified version of the stochastic gradient descent-based layout method was used in the uniform manifold approximation and projection (UMAP) algorithm (Becht et al., 2018), to visualize the ACTIONet graph. ACTIONet framework identifies dominant transcriptional patterns (or archetypes) and associates cells to these archetypes with different degrees of confidence. To discretize cell associations, we assigned each cell to its most closely associated archetype. We then used a curated set of reproducible cell-type-specific markers in the human prefrontal cortex²⁸ to annotate cells assigned to each archetype.

To filter cells, we performed two independent iterations of the ACTIONet framework, the first one to identify additional low-quality cells and missed doublets, and the second one to identify and annotate cell types/states. In the first stage, we removed two archetypes that were assigned to multiple unrelated marker sets, as well as cells that are not well-connected to other cells in the ACTIONet. To identify the overall connectivity of each cell within the ACTIONet graph that is assigned to a given archetype, we computed its “coreness” ⁷⁸ within the subgraph induced by the cells associated with the same archetype. In the second round, we annotated and grouped archetypes into major cell types in the human prefrontal cortex, including subtypes of excitatory and inhibitory neurons, as well as non-neuronal cell types, as well as two novel schizophrenia-associated cell states (Fig. 1b,d). We used the coreness of cells in the ACTIONet as their transparency to de-emphasize low-quality cells. To verify our annotations, we projected individual marker genes for different subtypes on the ACTIONet (Fig. 1c). We did not observe any batch effect and cells from all batches and both phenotypes were well-mixed in the ACTIONet plot (Extended Data Fig. 1b, c), and none of the archetypes is influenced by a minority of individuals, and the majority of individuals contribute to all archetypes (Fig. 1f).

Differential gene expression analysis using a modified pseudo-bulk approach

Following recent studies showing the success of pseudo-bulk methods for multi-sample multi-group differential analysis of single-cell datasets with complex experimental designs ⁷⁹, we developed a new method based on the analysis of pseudo-bulk profiles to identify perturbed genes in schizophrenia. We aggregated the expression of genes within each archetype/individual combination and used the linear-modeling approach in the Limma package⁸⁰ to include our experimental design in the differential analysis. However, unlike the case of bulk analysis, in which the underlying variance of each gene in each sample is unknown, in single-cell pseudo-bulk analysis we can compute both the mean and the variance of genes directly from single-cells, which provides a more accurate approach than either Limma-trend or Limma-voom⁸¹. We used the inverse of the gene variances as weights in our linear model and used an outlier detection approach to mask out genes that were deemed not reliable, on a per-sample basis. For each pseudo-bulk profile, we required that it has to contain at least 50 cells to be included in our model. Finally, to account for individual-specific differences, we incorporated age, gender, PMI, batch, and medication history as covariates in our model. We filtered results based on a raw p-value cut-off of 0.05 and LFR of 0.1 to declare differential genes.

Computing transcriptional pathology scores

To assess the degree to which different cell types/states in each individual are affected by schizophrenia-associated transcriptional perturbations, we used our pseudo-bulk profiles. For each archetype, we computed a representative expression vector by averaging all pseudo-bulk samples. As a measure of schizophrenia-associated transcriptional effect in each individual/cell type, we used partial Pearson correlation between individual pseudo-bulk profiles with differential perturbation scores, after controlling for the baseline cell type/state-specific representative vector. We then scored and ranked each individual scored according to their average correlation across all archetypes.

Identification of transcriptional regulators using ChEA3

We used ChIP-X Enrichment Analysis 3 (ChEA3)⁶⁵ to prioritize transcription factors (TFs) that are likely to mediate the observed transcriptional dysregulation in schizophrenia. In summary, ChEA3 integrates multiple libraries of putative TF-target lists, gathered from different sources, including TF-gene co-expression from RNA-seq studies, TF-target associations from ChIP-seq experiments, and TF-gene co-occurrence computed from crowd-submitted gene lists, to compute a composite ranked list of TFs. We used the TopRank strategy to combine rankings from individual libraries, in which the best scaled-rank of each TF across all libraries is used to aggregate TF scores. We report the log-transformed scaled-rank as the relevance measure of each TF by ChEA3 analysis.

Gene-centric analysis of common variants using H-MAGMA

Hi-C-coupled MAGMA (H-MAGMA) is a recent extension of the traditional multimarker analysis of genomic annotation (MAGMA), a method developed to prioritize genes by aggregating single nucleotide polymorphism associations to nearest genes. In H-MAGMA, the linking of SNPs to genes is extended to include long-range interactions brought together through the chromatin looping. In our analysis, we used the preprocessed dataset of SNPs-gene links based on the Hi-C data in the adult human brain.

Fluorescence in situ hybridization (FISH)

From our larger cohort of postmortem human BA10 tissue, six cases were selected for FISH experiments including CON10, CON14, SZ5, SZ12, SZ20, and SZ23. Frozen tissue blocks were embedded in optimal cutting temperature medium (OCT), and 10 μm sections were cut at −20 °C with a Microm HM560 cryostat, mounted on Superfrost Plus slides (one control, one schizophrenia, and one transcriptionally-resilient case mounted on each slide to ensure balanced processing) and stored −80°C until FISH was conducted.

Advanced Cell Diagnostics (ACD) designed the in situ hybridization probes (human SHANK2, UNC13A, TCF4, and CLU) as well as the positive and negative control probes (Supplementary Table 14). The RNAscope Multiplex Fluorescent Reagent Kit v2 (ACD) was used for the assay following the manufacturer’s instructions with some modifications. Nuclei were counterstained with DAPI using TSA buffer (ACD) and TSA Plus fluorophores (PerkinElmer). Briefly, frozen sections were fixed for 1 h at 4 °C using freshly prepared ice-cold 10% neutral buffered formalin, and rinsed with phosphate buffered saline (PBS) and dehydrated in 50%, 70%, and two changes of 100% ethanol (5 min each) at room temperature (RT). Sections were air-dried and a hydrophobic barrier drawn around each section with an Immedge pen (Vector Laboratories). When completely dry, sections were treated with hydrogen peroxide for 10 min at RT, washed twice with PBS, incubated with protease IV for 15 min at RT, and washed again twice with PBS.

After diluting the probes for RNA detection 1:50, sections were hybridized with the probes (40°C, 2 h; HybEZ Hybridization System (ACD)), washed twice, and stored overnight at RT in 5x SSC buffer. The following day, slides were washed twice with the wash buffer, followed by three amplification steps (AMP 1 (30 min), AMP 2 (30 min), and AMP 3 (15 min); 40°C). Each amplification step was followed by two washes of 2 min each with the wash buffer. Sections were then incubated sequentially with the HRP reagent corresponding to each channel (e.g. HRP-C1; 40°C, 15 min) followed by the respective TSA Plus fluorophore assigned to the probe channel (Opal dyes, 520, 620 and 690, dilution 1:1500; 40°C, 30 min) and HRP blocker (40°C, 15 min), and each incubation was followed by two wash steps. Lipofuscin autofluorescence was visible in both green and red channels. Since the far-red channel showed less autofluorescence, highly expressed UNC13A-C4 and CLU-C2 probes were assigned to the green fluorescein channel and the CUX2 layer-specific marker was assigned to the red cyanine 5 channel. DAPI (30 s) was used to visualize cell nuclei. Sections were mounted using ProLong Gold mounting medium (Thermo Fisher Scientific) and stored at 4 °C. Two independent experiments were performed, with three biological replicates each, and positive and negative control probes to test for RNA quality and background signal, respectively.

Leica Laser Scanning Confocal Microscopy

Image acquisition, processing, and quantitation were performed blind to diagnosis. All microscopy was performed at the Microscopy Core at McLean Hospital on a Leica TCS-SP8 confocal microscope. Exposure times were set separately for each of the four channels (red for Cy5, blue for DAPI, green for FITC, and orange for Texas Red) and kept similar among the cases to enable comparison. Images at 40x and 63x magnifications were visualized as maximum intensity projections of Z-stacks at 1.5 μm intervals. Two representative fields of view were selected within cortical layers II and III as identified by positive CUX2 staining, and Z-stacks were taken throughout the depth of single cells, with ≥40 single cells per case and six cases per target gene. Adjustment of brightness, contrast, and sharpness was done using Adobe Photoshop. Quantification of transcripts was performed in an automated unbiased manner using dotdotdot⁴⁰. For all six cases, the dotdotdot MATLAB script (https://github.com/LieberInstitute/dotdotdot) was used to process one 40x and one 63x maximum intensity projection image, counting the number of puncta in each channel overlapping with each region of interest (nucleus) identified by DAPI staining. Data presented in figure 2e represent counts of TCF4, CLU, SHANK2, and UNC13A transcripts within individual nuclei containing at least one detected CUX2 transcript, and excluding counts of zero.

CUT&Tag mapping of transcription factor binding in the neuronal genome

Nuclei were isolated from the postmortem human prefrontal cortex as described above and incubated with 1:1000 diluted anti-NeuN antibody (EMD Millipore MAB377X) with 0.5% BSA in PBS at 4°C with end-over-end rotation for 45 minutes. After staining nuclei were counterstained with propidium iodide and sorted on a BD FACSAria III Cell Sorter at Harvard University’s Bauer Core Facility. 100,000 neuronal nuclei were used as input for each Cleavage Under Targets and Tagmentation (CUT&Tag) assay using rabbit primary antibodies targeting MEF2C (Abcam ab211439), SATB2 (Abcam ab92446), SOX5 (Abcam ab94396), TCF4 (ProteinTech 22337-1-AP), H3K27Ac (EMD Millipore MABE647), and mouse-anti-rabbit secondary antibody (Sigma R2655) with the Vazyme pG-Tn5 CUT&Tag kit (Cellagen Technology, San Diego) according to the manufacturer’s protocol. CUT&Tag libraries were sequenced on one NextSeq500 flow cell at the MIT BioMicroCenter. Reads were aligned to the HG38 genome, processed to bedgraph format, and analyzed with the Sparse Enrichment for CUT&Run⁸² tool with “stringent” parameters retaining the top 1% of peaks.