Single-cell dissection of schizophrenia reveals neurodevelopmental-synaptic axis and transcriptional resilience =============================================================================================================== * W. Brad Ruzicka * Shahin Mohammadi * Jose Davila-Velderrain * Sivan Subburaju * Daniel Reed Tso * Makayla Hourihan * Manolis Kellis ## Abstract Schizophrenia is a devastating mental disorder with a high societal burden, complex pathophysiology, and diverse genetic and environmental risk factors. Its complexity, polygenicity, and small-effect-size and cell-type-specific contributors have hindered mechanistic elucidation and the search for new therapeutics. Here, we present the first single-cell dissection of schizophrenia, across 500,000+ cells from 48 postmortem human prefrontal cortex samples, including 24 schizophrenia cases and 24 controls. We annotate 20 cell types/states, providing a high-resolution atlas of schizophrenia-altered genes and pathways in each. We find neurons are the most affected cell type, with deep-layer cortico-cortical projection neurons and parvalbumin-expressing inhibitory neurons showing significant transcriptional changes converging on genetically-implicated regions. We discover a novel excitatory-neuron cell-state indicative of transcriptional resilience and enriched in schizophrenia subjects with less-perturbed transcriptional signatures. We identify key trans-acting factors as candidate drivers of observed transcriptional perturbations, including MEF2C, TCF4, SOX5, and SATB2, and map their binding patterns in postmortem human neurons. These factors regulate distinct gene sets underlying fetal neurodevelopment and adult synaptic function, bridging two leading models of schizophrenia pathogenesis. Our results provide the most detailed map to date for mechanistic understanding and therapeutic development in neuropsychiatric disorders. ## Introduction Schizophrenia afflicts young adults just as they approach their full potential, manifests as a combination of psychosis, social withdrawal, and cognitive dysfunction, and often leads to a lifetime of profound and chronic disability1. Schizophrenia pathogenesis is thought to begin during neurodevelopment, yet symptoms do not emerge until many years later in early adulthood2. The complex etiology of schizophrenia is driven by both genetic and environmental factors affecting a wide range of brain-related processes, including neurodevelopment3,4, cognition5,6, synaptic function7,8, neuronal excitability9,10, and neuronal connectivity11,12. Decades-long efforts to decipher schizophrenia genetics have yielded 145+ robustly-associated genetic loci13,14, but their target genes and cell-types of action remain largely uncharacterized, hindering the search for mechanistic insights and therapeutics development. Cell-type-specific profiles of reference (non-schizophrenia) brain-sample transcriptomes have shown that genes near schizophrenia risk loci are expressed in pyramidal neurons and specific subsets of inhibitory neurons15, but do not reveal whether schizophrenia-locus-proximal genes are indeed differentially-expressed in schizophrenia, or whether increased or decreased expression is risk-associated or protective. Conversely, schizophrenia case-control gene expression analyses using homogenized cortical tissue have been carried out at the bulk level16–19, revealing several differentially-expressed genes in schizophrenia. However, such bulk-level studies average gene expression across millions of cells, merging together diverse and inconsistent cell types, and thus miss genes that are differentially-expressed only in lower-abundance cell types and genes with opposite changes in different cell populations, and can also result in false positives stemming from cell-type-composition changes between samples. Emerging technologies for single-cell transcriptomics20,21 achieve both cell-type-specificity and reveal disease-associated changes, as demonstrated for Alzheimer’s Disease22, Autism Spectrum Disorder23, Major Depressive Disorder24, and Multiple Sclerosis25, but these have not been applied to schizophrenia to date. Here, we present the first transcriptomic atlas of schizophrenia at single-cell resolution, profiling >500,000 cells from postmortem human prefrontal cortex from 24 schizophrenia, and 24 age-matched control individuals. We annotate 18 distinct cell types, including seven excitatory and six inhibitory neuronal subpopulations and five non-neuronal cell types. We also identify a new excitatory neuron cell-state, Ex-SZTR, which is enriched for differentially-expressed genes and significantly more prevalent in schizophrenia than in control individuals, but surprisingly preferentially found in schizophrenia individuals with non-schizophrenia transcriptional signatures across all other cell types, indicating transcriptional resilience. Schizophrenia-associated transcriptional perturbations enrich in highly-specific cellular populations, including deep-layer cortico-cortical projection neurons and PV-expressing inhibitory neurons, highlighting the importance of cell-type-specific assessments. These changes preferentially perturb postsynaptic organization, synaptic plasticity, and neurodevelopmental pathways, providing mechanistic insights into functional pathways that underlie schizophrenia perturbations. Genes showing differential expression in schizophrenia are significantly enriched in the proximity of GWAS-associated genes14 and linked to schizophrenia-associated genetic variants via enhancer-promoter loops in adult dorsolateral prefrontal cortex26, providing candidate target genes, the directionality of effect, and cell-type-of-action for 68 of 145 schizophrenia loci. Searching for common upstream regulators for differentially-expressed genes, we find that schizophrenia-associated transcriptional perturbations converge on a small number of transcription factors, which are themselves encoded in schizophrenia-associated genetic loci, including MEF2C, SATB2, SOX5, and TCF4. These factors act in both early fetal brain development and in adult brain synaptic processes, thus linking two leading models of schizophrenia pathogenesis. Our results provide both a unique resource for understanding the molecular biology of schizophrenia at single-cell resolution, and numerous new insights for understanding candidate cell-type-specific driver genes, biological pathways, master regulators, and resilience mechanisms, and for prioritizing target genes for therapeutic development against this devastating disorder. ### 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)27. 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. ![Figure 1.](http://medrxiv.org/http://medrxiv.stage.highwire.org/content/medrxiv/early/2020/11/09/2020.11.06.20225342/F1.medium.gif) [Figure 1.](http://medrxiv.org/content/early/2020/11/09/2020.11.06.20225342/F1) 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 associations82. **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, ACTIONet28, 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 cortex28 (**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 population29. 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 signatures30. 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 brain23, 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 individuals31. Cells associated with the Ex-NRGN cell state highly express mitochondrial genes (**Extended Data Fig. 2c**), an observation reproducible in an independent study23 (**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 work32, 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**). ![Figure 2.](http://medrxiv.org/http://medrxiv.stage.highwire.org/content/medrxiv/early/2020/11/09/2020.11.06.20225342/F2.medium.gif) [Figure 2.](http://medrxiv.org/content/early/2020/11/09/2020.11.06.20225342/F2) 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 tissue28.**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 DotDotDot76. 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 samples19. We observed a significant overlap between the 2,450 upregulated (*p-*value<10−8, Fisher’s exact test) and 2,371 downregulated (*p*-value<10−29) 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**14*, 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-Seq18. *TCF4* plays critical roles in both neurodevelopment and in regulating the excitability of prefrontal cortex neurons33, and was previously predicted as a “master regulator” of schizophrenia gene network perturbations34. Additionally, the molecular chaperone *CLU*, which also lies in a schizophrenia-associated genetic locus14 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 samples35, and overexpressed in laser-capture-microdissected excitatory neurons and in parvalbumin-expressing inhibitory neurons36,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 function38. In contrast, *SHANK* genes encode postsynaptic scaffolding proteins essential for the organization of glutamatergic synapses39. ### 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 dotdotdot40) 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 organization41, synaptic plasticity2, and neurodevelopment42 (**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 density43,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**45*, *VAV3, WASF3**46*), and axon guidance (*ABLIM1, NCK1**47*, *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 III48–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 schizophrenia14, 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. ![Figure 3.](http://medrxiv.org/http://medrxiv.stage.highwire.org/content/medrxiv/early/2020/11/09/2020.11.06.20225342/F3.medium.gif) [Figure 3.](http://medrxiv.org/content/early/2020/11/09/2020.11.06.20225342/F3) 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 loci14 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 excitability51, 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 schizophrenia14, we further considered four psychiatric disorders that are known to share genetic risk factors52: major depressive disorder (MDD)53, bipolar disorder (BD)54, autism spectrum disorder (ASD)55, and attention-deficit/hyperactivity disorder (ADHD)56. This approach allows us to distinguish schizophrenia from general psychiatric illness-related associations. As a point of contrast, we included Alzheimer’s disease (AD)57, 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-MAGMA52, 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 cortex26. 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 overlaps52. We found weaker relationships between schizophrenia transcriptional perturbations and ASD risk52,58, consistent with the known lower correlation of genetic risk between these disorders58. 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 V59. 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 schizophrenia60,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 pathways62; *CLIP1*, a linker protein involved in microtubule trafficking of cargo within axons and dendrites63; *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 regulators64. 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 DEGs65 (**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. ![Figure 4.](http://medrxiv.org/http://medrxiv.stage.highwire.org/content/medrxiv/early/2020/11/09/2020.11.06.20225342/F4.medium.gif) [Figure 4.](http://medrxiv.org/content/early/2020/11/09/2020.11.06.20225342/F4) 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 timing66; 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 neurogenesis67. 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−7, **NDD** p-value: 3×10−6, 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 disorders68. We observed a significant overlap between genetically-identified NDD-associated TFs (de novo mutations and CNVs)69 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−7, 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)70. 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−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&Tag71) 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 dataset70, 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. ![Figure 5.](http://medrxiv.org/http://medrxiv.stage.highwire.org/content/medrxiv/early/2020/11/09/2020.11.06.20225342/F5.medium.gif) [Figure 5.](http://medrxiv.org/content/early/2020/11/09/2020.11.06.20225342/F5) 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 repressor72. 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 PEC70, 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. ## Discussion In this work, we presented the first single-cell transcriptomic case-control dissection of schizophrenia, producing an invaluable high-resolution and high-quality dataset of more than 500,000 single-cell transcriptomes. We annotated 18 neuronal and glial cell types and two excitatory-neuron cross-cutting cellular states, which we used to investigate cell-type-specific schizophrenia-dysregulated genes, pathways, and regulators. The vast majority of differentially-expressed genes were cell-type-specific and undetectable in bulk datasets, including 387 genes that were both upregulated and downregulated in distinct cell populations, highlighting the importance of our single-cell datasets. Within neuronal subpopulations, we uncovered both wide-spread and subtype-specific transcriptional changes in both glutamatergic and GABAergic neurons, in particular PV-expressing interneurons, that are not readily detectable in bulk tissue. Our results reveal the central role in schizophrenia for genes and regulators involved in synapse formation/structure/function, with synaptic dysregulation of different neuronal subtypes implicating distinct but functionally-related genes. Unexpectedly, the transcription factors targeting these synapse-related genes enrich primarily neurodevelopmental processes, linking disruption of early neurodevelopment and adult synaptic dysfunction and providing a novel bridge between two prominent theories of schizophrenia pathogenesis. Furthermore, across multiple neuronal populations, upregulated genes implicated fetal enhancers while downregulated genes implicated adult enhancers, suggesting impairment of epigenomic differentiation of brain-specific enhancers across development as an underlying mechanistic paradigm. Cutting across multiple subtypes of excitatory neurons within multiple layers, we uncovered a new cellular state, Ex-SZTR, overrepresented in schizophrenic individuals but surprisingly associated with transcriptional signatures of non-schizophrenia individuals across all other neuronal cell types, indicative of transcriptional resilience. This relationship between gene-regulatory changes in only one cellular sub-state and global transcriptional dysregulation in other neuronal cell types is reminiscent of the role of PV-expressing inhibitory neurons in regulating synchronous firing of assemblies of excitatory neurons, a process implicated in multiple neurodegenerative and psychiatric disorders73,74. Indeed, we found that PV interneurons showed the strongest correlation with non-schizophrenia (“resilience”) transcriptional signatures for the most transcriptionally-resilient individuals, suggesting a potential interplay between Ex-SZTR and PV-expressing interneurons. While the genetic and environmental factors influencing the Ex-SZTR cell state require further investigation, these tantalizing findings suggest a new potential mechanism of resilience to schizophrenia molecular pathologies that may be exploitable for new whole-brain therapeutic interventions against schizophrenia. Our cell-type-specific DEGs showed a highly-significant enrichment in genetic loci associated with schizophrenia, with schizophrenia-differentially-expressed genes in approximately half of all schizophrenia GWAS loci14, including many loci previously lacking any mechanistic hypothesis. This substantial overlap suggests potential mediating roles of these genes in linking genetic causation and molecular phenotypic manifestation and helps reveal putative target genes that may be mediating the genetic effects of these loci, the cell type in which these genes may be acting, and the directionality of their effect to distinguish protective vs. risk-increasing gene expression changes. These insights can be invaluable in guiding the development of therapeutic interventions, deciding on genes to target, the cell types in which to observe their impact, and whether pharmaceutical interventions should be inhibitory or inducing, especially given the greatly-increased success of therapeutic efforts with genetic support75. Across GWAS loci, deep-layer excitatory neurons, cortico-cortical projection neurons, and PV-expressing basket interneurons showed the strongest enrichments, indicating central roles in mediating causal genetic effects. These enrichments between our differentially-expressed genes and GWAS loci held across multiple psychiatric disorders, consistent with the shared heritability between them58, but were absent from Alzheimer’s GWAS loci, providing confidence about the specificity of our analyses. These high-resolution cell-type-specific dysregulated gene sets also enabled insights into the upstream transcription factors (TFs) most likely to influence these changes. Unlike dysregulated genes, these factors are predominantly enriched for neurodevelopmental pathways. In particular, multiple lines of evidence reveal a central role for four master regulators, namely TCF4, MEF2C, SOX5, and SATB2, which: (1) regulate schizophrenia dysregulated genes primarily in neurons; (2) are genetically associated with both schizophrenia and NDD; and (3) are key neurodevelopmental regulatory factors70. We experimentally validated that the target genes of these factors are schizophrenia DEGs, confirming their relevance to schizophrenia pathology. Because these TFs are major regulators of neurodevelopment with strong evidence for their action as upstream regulators of neuronal schizophrenia DEGs in adult synaptic function, we propose that their pleiotropic roles may represent a link between early neurodevelopmental disruptions and adult brain function. How genetic or early environmental perturbations to these TFs impact function to increase the risk for schizophrenia, and why it is only later in life that the phenotypic effects manifest, cannot be answered at present. However, our observations open up opportunities for future mechanistic studies tracking the developmental consequences of direct genetic perturbations to the candidate schizophrenia regulators and target genes put forward here. While providing numerous biological insights, we acknowledge that the scope of this work is limited by current technologies. As is necessary for investigating postmortem brain tissue, we measured transcriptomes of isolated single-nuclei, and the differences between nuclear and whole-tissue RNA content must be considered in the interpretation of this work. Additionally, there are many mechanisms critical to brain function not visible to this methodology, such as mRNA splicing or trafficking of resident dendritic mRNAs. Another technical limitation of this work is the inability to spatially resolve differential expression events, and we have addressed this problem computationally and histologically where possible. These technology-specific limitations provide an opportunity for future research when newer, more advanced technologies are routinely available. Additionally, while our case-control design is a strength of the study, this focus on patient tissue does not allow experimental manipulation to investigate the causality of our validated observations, and experiments in model systems are needed. Individuals with schizophrenia have very diverse experiences of treatment and clinical outcomes. Our study design is well powered for comparison of schizophrenia vs. control, but identification of trends within subgroups of schizophrenia subjects affected by specific environmental, genetic, and treatment factors will require larger numbers of subjects. This work utilized a rigorous study design and analysis plan to address potential confounds that are common to investigations of postmortem human brain tissue. Through balanced inclusion of schizophrenia and control subjects in each multiplexed sample preparation and sequencing library, we observed a remarkably low batch effect in our data, and almost no doublet contamination, two common problems with microfluidics-based snRNA-seq. We balanced diagnostic groups for demographic variables, including age, gender, and postmortem interval, and also controlled for these variables by including them as covariates in our analyses. We also controlled for psychiatric medications by including them as quantitative covariates in our analysis, as nearly all chronic schizophrenia patients receive longstanding pharmacologic treatment while individuals without psychiatric illness do not, and thus psychiatric diagnostic groups are not easily balanced for psychiatric medication exposures. The data presented here offer a cell-type-specific reframing of schizophrenia transcriptional pathology by revealing specific cell populations impacted by schizophrenia genes, variants, and regulators. Identification of pleiotropic transcriptional regulators linking developmental and adult schizophrenia-associated pathologies, and the newly discovered transcriptionally-resilient cell state, provide avenues for future research to unravel the genetic and environmental underpinnings of this complex and heterogeneous disease, with promise for improving outcomes and quality of life for patients and their families. ## Materials and methods ### 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 study22. 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-Seq27 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 package27 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 package76 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 77 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) framework28 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 cortex28 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” 78 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 79, 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 package80 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-voom81. 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)65 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 dotdotdot40. For all six cases, the dotdotdot MATLAB script ([https://github.com/LieberInstitute/dotdotdot](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&Run82 tool with “stringent” parameters retaining the top 1% of peaks. ## Supporting information TableS2_SequencingMetrics [[supplements/225342_file02.xlsx]](pending:yes) TableS3\_Combined\_marker_tables [[supplements/225342_file03.xlsx]](pending:yes) TableS4\_CCab\_scores [[supplements/225342_file04.xlsx]](pending:yes) TableS5\_Combined\_marker\_tables\_enrichments [[supplements/225342_file05.xlsx]](pending:yes) TableS6a\_DE\_up_genes [[supplements/225342_file06.xlsx]](pending:yes) TableS6b\_DE\_down_genes [[supplements/225342_file07.xlsx]](pending:yes) TableS7c\_TPS\_scores [[supplements/225342_file08.xlsx]](pending:yes) TableS8a_NFP [[supplements/225342_file09.xlsx]](pending:yes) TableS8b\_DE\_up_gprofiler [[supplements/225342_file10.xlsx]](pending:yes) TableS8c\_DE\_down_gprofiler [[supplements/225342_file11.xlsx]](pending:yes) TableS9\_GWAS\_SNP_table [[supplements/225342_file12.xlsx]](pending:yes) TableS10\_ChEA\_scores [[supplements/225342_file13.xlsx]](pending:yes) TableS11a\_ChEA\_TFs\_enrichment\_NFP [[supplements/225342_file14.xlsx]](pending:yes) TableS11b\_ChEA\_TFs\_up\_gprofiler [[supplements/225342_file15.xlsx]](pending:yes) TableS11c\_ChEA\_TFs\_down\_gprofiler [[supplements/225342_file16.xlsx]](pending:yes) TableS13_Regulons [[supplements/225342_file17.xlsx]](pending:yes) TableS14_RNAScopeProbes [[supplements/225342_file18.xlsx]](pending:yes) TableS7ac\_Pseudobulk\_mean_Ex-L2-3 [[supplements/225342_file19.xlsx]](pending:yes) TableS7aa\_Pseudobulk\_mean_Ex-NRGN [[supplements/225342_file20.xlsx]](pending:yes) TableS7ab\_Pseudobulk\_mean_Ex-SZTR [[supplements/225342_file21.xlsx]](pending:yes) TableS7ad\_Pseudobulk\_mean_Ex-L4 [[supplements/225342_file22.xlsx]](pending:yes) TableS7ae\_Pseudobulk\_mean_Ex-L4-5 [[supplements/225342_file23.xlsx]](pending:yes) TableS7af\_Pseudobulk\_mean_Ex-L5 [[supplements/225342_file24.xlsx]](pending:yes) TableS7ag\_Pseudobulk\_mean_Ex-L5-6 [[supplements/225342_file25.xlsx]](pending:yes) TableS7ah\_Pseudobulk\_mean_Ex-L5-6CCa [[supplements/225342_file26.xlsx]](pending:yes) TableS7ai\_Pseudobulk\_mean_Ex-L5-6CCb [[supplements/225342_file27.xlsx]](pending:yes) TableS7aj\_Pseudobulk\_mean_In-Rosehip [[supplements/225342_file28.xlsx]](pending:yes) TableS7ak\_Pseudobulk\_mean_In-VIP [[supplements/225342_file29.xlsx]](pending:yes) TableS7al\_Pseudobulk\_mean_In-Reelin [[supplements/225342_file30.xlsx]](pending:yes) TableS7am\_Pseudobulk\_mean_In-PV-Basket [[supplements/225342_file31.xlsx]](pending:yes) TableS7an\_Pseudobulk\_mean_In-PV-Chandelier [[supplements/225342_file32.xlsx]](pending:yes) TableS7ao\_Pseudobulk\_mean_In-SST [[supplements/225342_file33.xlsx]](pending:yes) TableS7ap\_Pseudobulk\_mean_AST [[supplements/225342_file34.xlsx]](pending:yes) TableS7aq\_Pseudobulk\_mean_OPC [[supplements/225342_file35.xlsx]](pending:yes) TableS7as\_Pseudobulk\_mean_Mic [[supplements/225342_file36.xlsx]](pending:yes) TableS7at\_Pseudobulk\_mean_Endo [[supplements/225342_file37.xlsx]](pending:yes) TableS7ar\_Pseudobulk\_mean_Oli [[supplements/225342_file38.xlsx]](pending:yes) TableS7ba\_Pseudobulk\_variance_Ex-NRGN [[supplements/225342_file39.xlsx]](pending:yes) TableS7bb\_Pseudobulk\_variance_Ex-SZTR [[supplements/225342_file40.xlsx]](pending:yes) TableS7bc\_Pseudobulk\_variance_Ex-L2-3 [[supplements/225342_file41.xlsx]](pending:yes) TableS7bd\_Pseudobulk\_variance_Ex-L4 [[supplements/225342_file42.xlsx]](pending:yes) TableS7be\_Pseudobulk\_variance_Ex-L4-5 [[supplements/225342_file43.xlsx]](pending:yes) TableS7bf\_Pseudobulk\_variance_Ex-L5 [[supplements/225342_file44.xlsx]](pending:yes) TableS7bg\_Pseudobulk\_variance_Ex-L5-6 [[supplements/225342_file45.xlsx]](pending:yes) TableS7bh\_Pseudobulk\_variance_Ex-L5-6CCa [[supplements/225342_file46.xlsx]](pending:yes) TableS7bi\_Pseudobulk\_variance_Ex-L5-6CCb [[supplements/225342_file47.xlsx]](pending:yes) TableS7bj\_Pseudobulk\_variance_In-Rosehip [[supplements/225342_file48.xlsx]](pending:yes) TableS7bk\_Pseudobulk\_variance_In-VIP [[supplements/225342_file49.xlsx]](pending:yes) TableS7bl\_Pseudobulk\_variance_In-Reelin [[supplements/225342_file50.xlsx]](pending:yes) TableS7bm\_Pseudobulk\_variance_In-PV-Basket [[supplements/225342_file51.xlsx]](pending:yes) TableS7bn\_Pseudobulk\_variance_In-PV-Chandelier [[supplements/225342_file52.xlsx]](pending:yes) TableS7bo\_Pseudobulk\_variance_In-SST [[supplements/225342_file53.xlsx]](pending:yes) TableS7bp\_Pseudobulk\_variance_AST [[supplements/225342_file54.xlsx]](pending:yes) TableS7bq\_Pseudobulk\_variance_OPC [[supplements/225342_file55.xlsx]](pending:yes) TableS7br\_Pseudobulk\_variance_Oli [[supplements/225342_file56.xlsx]](pending:yes) TableS7bs\_Pseudobulk\_variance_Mic [[supplements/225342_file57.xlsx]](pending:yes) TableS7bt\_Pseudobulk\_variance_Endo [[supplements/225342_file58.xlsx]](pending:yes) TableS12\_FilteredNorm\_JointPseudobulk_part01 [[supplements/225342_file59.zip]](pending:yes) TableS12\_FilteredNorm\_JointPseudobulk_part02 [[supplements/225342_file60.zip]](pending:yes) TableS12\_FilteredNorm\_JointPseudobulk_part03 [[supplements/225342_file61.zip]](pending:yes) TableS12\_FilteredNorm\_JointPseudobulk_part04 [[supplements/225342_file62.zip]](pending:yes) TableS12\_FilteredNorm\_JointPseudobulk_part05 [[supplements/225342_file63.zip]](pending:yes) TableS12\_FilteredNorm\_JointPseudobulk_part06 [[supplements/225342_file64.zip]](pending:yes) TableS12\_FilteredNorm\_JointPseudobulk_part07 [[supplements/225342_file65.zip]](pending:yes) TableS12\_FilteredNorm\_JointPseudobulk_part08 [[supplements/225342_file66.zip]](pending:yes) TableS12\_FilteredNorm\_JointPseudobulk_part09 [[supplements/225342_file67.zip]](pending:yes) TableS12\_FilteredNorm\_JointPseudobulk_part10 [[supplements/225342_file68.zip]](pending:yes) TableS12\_FilteredNorm\_JointPseudobulk_part11 [[supplements/225342_file69.zip]](pending:yes) TableS12\_FilteredNorm\_JointPseudobulk_part12 [[supplements/225342_file70.zip]](pending:yes) TableS12\_FilteredNorm\_JointPseudobulk_part13 [[supplements/225342_file71.zip]](pending:yes) TableS12\_FilteredNorm\_JointPseudobulk_part14 [[supplements/225342_file72.zip]](pending:yes) TableS1_SubjectDemographics&Medication [[supplements/225342_file73.xlsx]](pending:yes) ## Data Availability MULTI-seq data are available at Synapse ([https://www.synapse.org/#!Synapse:syn22963646](https://www.synapse.org/#!Synapse:syn22963646)). The data are available under controlled use conditions set by human data privacy regulations, and access requires a data use agreement to ensure the anonymity of the donors of postmortem human brain tissue. A data use agreement can be pursued with SAGE, who maintains Synapse and can be downloaded from their website. ## Author Contributions This study was designed by W.B.R., and directed and coordinated by W.B.R. and M.K. W.B.R., S.S., and D.R.T. performed the snRNA-seq experiment, S.S. performed the RNAScope experiment and W.B.R. and S.S. performed the CUT&Tag experiment. W.B.R. performed data processing and S.M. and J.D.-V. performed the computational analysis. D.R.T. and M.H. reviewed medical records under supervision of W.B.R. W.B.R., S.M., J.D.-V., and M.K. wrote the manuscript. ## Data Availability The MULTI-seq data are available at Synapse ([https://www.synapse.org/#!Synapse:syn22963646](https://www.synapse.org/#!Synapse:syn22963646)). The data are available under controlled use conditions set by human data privacy regulations, and access requires a data use agreement to ensure the anonymity of the donors of postmortem human brain tissue. A data use agreement can be pursued with SAGE, who maintains Synapse and can be downloaded from their website. ## Supplementary Tables 1. Subject demographics and medications (TableS1_SubjectDemographics&Medication.xlsx) 2. Sequencing statistics (TableS2_SequencingMetrics.xlsx) 3. Cell-type/state-specific marker genes (TableS3_Combined_marker_tables.xlsx) 4. Marker gene comparison between Ex-L5/6-CCa and Ex-L5/6-CCb (TableS4_CCab_scores.xlsx) 5. Pathways enriched in top-ranked marker genes from each cell-type/state (TableS5\_Combined_marker_tables_enrichments.xlsx) 6. Up-(TableS6a\_up\_regulated\_genes.xlsx) and down- (TableS6b_down_regulated_genes.xlsx) regulated genes in schizophrenia 7. Pseudo-bulk mean (TableS7a\_Pseudobulk\_mean.xlsx) and variance (TableS7b\_Pseudobulk_variance.xlsx) profile per cell type/state, and computed Transcriptional Pathology Scores for each subject (TableS7c_TPS_scores.xlsx) 8. Curated set of neuro-associated functional pathways (NFP) and member genes (TableS8a\_NFP.xlsx), and functional pathway enrichment of up- (TableS8b\_DE\_up_gprofiler.xlsx) and down- (TableS8c_DE_down_gprofiler.xlsx) regulated genes 9. Schizophrenia GWAS loci linked to DEGs (TableS9_GWAS_SNP_table.xlsx) 10. ChEA scores for association of 1632 TFs with up- and downregulated genes (TableS10_ChEA_scores.xlsx) 11. NFP enrichment of TFs associated with dysregulated genes (TableS11a\_ChEA\_TFs\_enrichment\_NFP.xlsx), and functional pathway enrichment of TFs associated with the set of up- (TableS11b\_ChEA\_TFs\_up\_gprofiler.xlsx) and down- (TableS11c\_ChEA_TFs_down_gprofiler.xlsx) regulated genes 12. Filtered and normalized profile of the merged pseudo-bulk mean profiles (TableS12_FilteredNorm_JointPseudobulk.xlsx), used for TF expression correlation analysis 13. TF membership of individual regulons (TableS13_Regulons.xlsx) 14. RNAscope targets, probes, channels, and fluorophores (TableS14_RNAScopeProbes.xlsx) ![Extended Data Figure 1a.](http://medrxiv.org/http://medrxiv.stage.highwire.org/content/medrxiv/early/2020/11/09/2020.11.06.20225342/F6.medium.gif) [Extended Data Figure 1a.](http://medrxiv.org/content/early/2020/11/09/2020.11.06.20225342/F6) Extended Data Figure 1a. ACTIONet cell-cell similarity network depicting the footprint of all 20 identified transcriptional archetypes **ACTIONet cell-cell similarity network**. Network-based two-dimensional visualization of all cells considered in the analysis (n=560,020) indicating association with identified transcriptional archetypes. ![Extended Data Figure 1b,c.](http://medrxiv.org/http://medrxiv.stage.highwire.org/content/medrxiv/early/2020/11/09/2020.11.06.20225342/F7.medium.gif) [Extended Data Figure 1b,c.](http://medrxiv.org/content/early/2020/11/09/2020.11.06.20225342/F7) Extended Data Figure 1b,c. ACTIONet cell-cell similarity network colored by 10x batch (left) or phenotype (right) ACTIONet cell-cell similarity network. Network-based two-dimensional visualization of all cells considered in the analysis (n=560,020) indicating association with 10x batch (b), or phenotype (c) (SZ:schizophrenia, n=266,431; CON: control, n=293,589). ![Extended Data Figure 2a.](http://medrxiv.org/http://medrxiv.stage.highwire.org/content/medrxiv/early/2020/11/09/2020.11.06.20225342/F8.medium.gif) [Extended Data Figure 2a.](http://medrxiv.org/content/early/2020/11/09/2020.11.06.20225342/F8) Extended Data Figure 2a. Gene capture across cell-types **Cell type gene and cell statistics**. Gene count distribution across cells of each type. Each point represents in log scale the number of genes detected to have a read count x > 0 in a given cell. Box plots are centred around the median, with the interquartile range (IQR) defining the box. The upper whisker extends to the largest value no further than 1.5 × IQR from the end of the box. The lower whisker extends to the smallest value at most 1.5 × IQR from the end of the box. ![Extended Data Figure 2b.](http://medrxiv.org/http://medrxiv.stage.highwire.org/content/medrxiv/early/2020/11/09/2020.11.06.20225342/F9.medium.gif) [Extended Data Figure 2b.](http://medrxiv.org/content/early/2020/11/09/2020.11.06.20225342/F9) Extended Data Figure 2b. UMI capture across cell-types **Cell type gene and cell statistics**. Total UMI count distribution across cells of each type. Each point represents the total UMI count across all genes in a given cell. Box plots are centred around the median, with the interquartile range (IQR) defining the box. The upper whisker extends to the largest value no further than 1.5 × IQR from the end of the box. The lower whisker extends to the smallest value at most 1.5 × IQR from the end of the box. ![Extended Data Figure 2c.](http://medrxiv.org/http://medrxiv.stage.highwire.org/content/medrxiv/early/2020/11/09/2020.11.06.20225342/F10.medium.gif) [Extended Data Figure 2c.](http://medrxiv.org/content/early/2020/11/09/2020.11.06.20225342/F10) Extended Data Figure 2c. Mitochondrial gene capture across cell-types **Cell type gene and cell statistics**. Distribution of UMI percentages that map to mitochondrially encoded genes across cells of each type. Distributions are shown for the SCZ dataset reported herein. Box plots are centred around the median, with the interquartile range (IQR) defining the box. The upper whisker extends to the largest value no further than 1.5 × IQR from the end of the box. The lower whisker extends to the smallest value at most 1.5 × IQR from the end of the box. ![Extended Data Figure 2d.](http://medrxiv.org/http://medrxiv.stage.highwire.org/content/medrxiv/early/2020/11/09/2020.11.06.20225342/F11.medium.gif) [Extended Data Figure 2d.](http://medrxiv.org/content/early/2020/11/09/2020.11.06.20225342/F11) Extended Data Figure 2d. Mitochondrial gene capture across cell-types in the Velmeshev 2019 dataset **Cell type gene and cell statistics**. Distribution of UMI percentages that map to mitochondrially encoded genes across cells of each type. Distributions are shown for the dataset reported in Velmeshev et al. 2019. Box plots are centred around the median, with the interquartile range (IQR) defining the box. The upper whisker extends to the largest value no further than 1.5 × IQR from the end of the box. The lower whisker extends to the smallest value at most 1.5 × IQR from the end of the box. ![Extended Data Figure 3.](http://medrxiv.org/http://medrxiv.stage.highwire.org/content/medrxiv/early/2020/11/09/2020.11.06.20225342/F12.medium.gif) [Extended Data Figure 3.](http://medrxiv.org/content/early/2020/11/09/2020.11.06.20225342/F12) Extended Data Figure 3. Rank-Rank Hypergeometric Overlap Plot of cell-type-specific trancriptomic changes Rank-Rank Hypergeometric Overlap plot depicting the similarity of transcriptional perturbations between all pairs of cell-types/states. ![Extended Data Figure 4.](http://medrxiv.org/http://medrxiv.stage.highwire.org/content/medrxiv/early/2020/11/09/2020.11.06.20225342/F13.medium.gif) [Extended Data Figure 4.](http://medrxiv.org/content/early/2020/11/09/2020.11.06.20225342/F13) Extended Data Figure 4. Transcriptional pathology scores across subject classification. Transcriptional pathology scores across all individuals demonstrate clear gradation across classifications, with SZTR individuals (Ex-SZTR cell fraction > 1%) ranking below the majority of CON individuals, away from the SZ group. ![Extended Data Figure 5.](http://medrxiv.org/http://medrxiv.stage.highwire.org/content/medrxiv/early/2020/11/09/2020.11.06.20225342/F14.medium.gif) [Extended Data Figure 5.](http://medrxiv.org/content/early/2020/11/09/2020.11.06.20225342/F14) Extended Data Figure 5. Fluorescence *in situ* hybridization of selected SZ DEGs. Representative photomicrographs from one case in each phenotype (SZTR, CON, SZ) are shown of in situ hybridization for the *CUX2* (Red, layer II and III excitatory neuron marker), *TCF4* (top, blue), *CLU* (middle, green), *SHANK2* (middle, blue), and *UNC13A* genes (bottom, green). Each image was captured at 40x magnification and for each image the full field of view is shown to the left, and the area boxed in blue is shown to the right. Each gene was assessed in two cases from each phenotype for six total cases. ![Extended Data Figure 6a.](http://medrxiv.org/http://medrxiv.stage.highwire.org/content/medrxiv/early/2020/11/09/2020.11.06.20225342/F15.medium.gif) [Extended Data Figure 6a.](http://medrxiv.org/content/early/2020/11/09/2020.11.06.20225342/F15) Extended Data Figure 6a. Volcano plots of SZ GWAS significance versus cell-type-specific transcriptional perturbation in excitatory neuronal populations **Genetic versus transcriptional SZ-associated perturbations**. Volcano plots showing the relationship between SZ GWAS scores18 (y-axis, −log10 association p-value) and SZ transcriptional perturbations measured and reported herein (x-axis, Log2 fold change of expression values in SZ relative to control samples). Plots are shown independently for subpopulations of excitatory neurons. ![Extended Data Figure 6b.](http://medrxiv.org/http://medrxiv.stage.highwire.org/content/medrxiv/early/2020/11/09/2020.11.06.20225342/F16.medium.gif) [Extended Data Figure 6b.](http://medrxiv.org/content/early/2020/11/09/2020.11.06.20225342/F16) Extended Data Figure 6b. Volcano plots of SZ GWAS significance versus cell-type-specific transcriptional perturbation in inhibitory neuronal populations **Genetic versus transcriptional SZ-associated perturbations**. Volcano plots showing the relationship between SZ GWAS scores18 (y-axis, −log10 association p-value) and SZ transcriptional perturbations measured and reported herein (x-axis, Log2 fold change of expression values in SZ relative to control samples). Plots are shown independently for subpopulations of inhibitory neuron subpopulations. ![Extended Data Figure 6c.](http://medrxiv.org/http://medrxiv.stage.highwire.org/content/medrxiv/early/2020/11/09/2020.11.06.20225342/F17.medium.gif) [Extended Data Figure 6c.](http://medrxiv.org/content/early/2020/11/09/2020.11.06.20225342/F17) Extended Data Figure 6c. Volcano plots of SZ GWAS significance versus cell-type-specific transcriptional perturbation in non-neuronal populations **Genetic versus transcriptional SZ-associated perturbations**. Volcano plots showing the relationship between SZ GWAS scores18 (y-axis, −log10 association p-value) and SZ transcriptional perturbations measured and reported herein (x-axis, Log2 fold change of expression values in SZ relative to control samples). Plots are shown independently for subpopulations of non-neuronal subpopulations. ![Extended Data Figure 7.](http://medrxiv.org/http://medrxiv.stage.highwire.org/content/medrxiv/early/2020/11/09/2020.11.06.20225342/F18.medium.gif) [Extended Data Figure 7.](http://medrxiv.org/content/early/2020/11/09/2020.11.06.20225342/F18) Extended Data Figure 7. Enrichment of SZ GWAS TF targets within SZ DEGs. Overrepresentation analysis (hypergeometric test) within targets of TFs genetically associated with SZ in GWAS (columns) of genes detected as differentially expressed in SZ relative to controls (rows). Overrepresentation analysis was performed independently for SZ upregulated (left) and downregulated genes (right). ![Extended Data Figure 8.](http://medrxiv.org/http://medrxiv.stage.highwire.org/content/medrxiv/early/2020/11/09/2020.11.06.20225342/F19.medium.gif) [Extended Data Figure 8.](http://medrxiv.org/content/early/2020/11/09/2020.11.06.20225342/F19) Extended Data Figure 8. Enrichment of neurodevelopmental delay-associated TF targets within SZ DEGs. Overrepresentation analysis (hypergeometric test) within targets of TFs genetically associated with neurodevelopmental delay65 (de novo mutations and CNVs) (columns) of genes detected as differentially expressed in SZ relative to controls (rows). Overrepresentation analysis was performed independently for SZ upregulated (left) and downregulated genes (right). ![Extended Data Figure 9.](http://medrxiv.org/http://medrxiv.stage.highwire.org/content/medrxiv/early/2020/11/09/2020.11.06.20225342/F20.medium.gif) [Extended Data Figure 9.](http://medrxiv.org/content/early/2020/11/09/2020.11.06.20225342/F20) Extended Data Figure 9. Neurodevelopment and neurodevelopmental delay associated genes overlap. Overlap of genes genetically associated with neurodevelopmental delay (NDD) and genes functionally annotated as related with neurodevelopment. SZ TFs are associated with both genesets. The majority of the lead NDD TFs (7/10) are also involved in neurodevelopment. ![Extended Data Figure 10.](http://medrxiv.org/http://medrxiv.stage.highwire.org/content/medrxiv/early/2020/11/09/2020.11.06.20225342/F21.medium.gif) [Extended Data Figure 10.](http://medrxiv.org/content/early/2020/11/09/2020.11.06.20225342/F21) Extended Data Figure 10. Overlap of PEC developing human brain regulatory elements and observed TF binding sites in neuronal nuclei. The overlap of the total number of peaks for each transcription factor, defined as the union of top 1% peaks observed in any of the samples, with the masterset of PEC developing human brain dataset, containing H3K27Ac peaks for human brain samples at different stages of development. ## Acknowledgements We thank the individuals and families whose donation of human brain tissue made this work possible. This work was supported by the Wilf Family Foundations and by NIH grant K08MH109759 (W.B.R.) and by NIH grants 1U01MH119509, R01MH109978, and R01AG062335 (M.K.). ## Footnotes * * e-mail: wruzicka{at}mclean.harvard.edu; manoli{at}mit.edu * Received November 6, 2020. * Revision received November 6, 2020. * Accepted November 9, 2020. * © 2020, Posted by Cold Spring Harbor Laboratory The copyright holder for this pre-print is the author. All rights reserved. The material may not be redistributed, re-used or adapted without the author's permission. ## References 1. 1.Kahn, R. S. et al. Schizophrenia. Nat Rev Dis Primer 1, 15067 (2015). 2. 2.Forsyth, J. K. & Lewis, D. A. Mapping the Consequences of Impaired Synaptic Plasticity in Schizophrenia through Development: An Integrative Model for Diverse Clinical Features. Trends Cogn. Sci. 21, 760–778 (2017). 3. 3.Weinberger, D. R. On the plausibility of ‘the neurodevelopmental hypothesis’ of schizophrenia. Neuropsychopharmacology 14, 1S–11S (1996). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/0893-133X(95)00199-N&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=8866738&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 4. 4.Weinberger, D. R. Future of Days Past: Neurodevelopment and Schizophrenia. Schizophr. Bull. 43, 1164–1168 (2017). 5. 5.Kahn, R. S. & Keefe, R.S.E. Schizophrenia is a cognitive illness: time for a change in focus. JAMA Psychiatry 70, 1107–1112 (2013). 6. 6.Toulopoulou, T. et al. Polygenic risk score increases schizophrenia liability through cognition-relevant pathways. Brain 142, 471–485 (2019). 7. 7.Glessner, J. T. et al. Strong synaptic transmission impact by copy number variations in schizophrenia. Proc. Natl. Acad. Sci. U. S. A. 107, 10584–10589 (2010). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NDoicG5hcyI7czo1OiJyZXNpZCI7czoxMjoiMTA3LzIzLzEwNTg0IjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjAvMTEvMDkvMjAyMC4xMS4wNi4yMDIyNTM0Mi5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 8. 8.Schwarz, E., Izmailov, R., Liò, P. & Meyer-Lindenberg, A. Protein Interaction Networks Link Schizophrenia Risk Loci to Synaptic Function. Schizophr. Bull. 42, 1334–1342 (2016). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/schbul/sbw035&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=27056717&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 9. 9.Nanou, E. & Catterall, W. A. Calcium Channels, Synaptic Plasticity, and Neuropsychiatric Disease. Neuron 98, 466–481 (2018). 10. 10.Mäki-Marttunen, T. et al. Functional Effects of Schizophrenia-Linked Genetic Variants on Intrinsic Single-Neuron Excitability: A Modeling Study. Biol Psychiatry Cogn Neurosci Neuroimaging 1, 49–59 (2016). 11. 11.Cao, H., Zhou, H. & Cannon, T. D. Functional connectome-wide associations of schizophrenia polygenic risk. Mol. Psychiatry (2020) doi:10.1038/s41380-020-0699-3. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41380-020-0699-3&link_type=DOI) 12. 12.Sekar, A. et al. Schizophrenia risk from complex variation of complement component 4. Nature 530, 177–183 (2016). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nature16549&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=26814963&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 13. 13.Schizophrenia Working Group of the Psychiatric Genomics Consortium. Biological insights from 108 schizophrenia-associated genetic loci. Nature 511, 421–427 (2014). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nature13595&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25056061&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000339335700037&link_type=ISI) 14. 14.Pardiñas, A. F. et al. Common schizophrenia alleles are enriched in mutation-intolerant genes and in regions under strong background selection. Nat. Genet. 50, 381–389 (2018). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41588-018-0059-2&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=29483656&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 15. 15.Skene, N. G. et al. Genetic identification of brain cell types underlying schizophrenia. Nat. Genet. 50, 825–833 (2018). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41588-018-0129-5&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=29785013&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 16. 16.Guidotti, A. et al. Decrease in reelin and glutamic acid decarboxylase67 (GAD67) expression in schizophrenia and bipolar disorder: a postmortem brain study. Arch. Gen. Psychiatry 57, 1061–1069 (2000). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1001/archpsyc.57.11.1061&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=11074872&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000165055700006&link_type=ISI) 17. 17.Akbarian, S. et al. Gene expression for glutamic acid decarboxylase is reduced without loss of neurons in prefrontal cortex of schizophrenics. Arch. Gen. Psychiatry 52, 258–266 (1995). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1001/archpsyc.1995.03950160008002&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=7702443&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1995QQ87500002&link_type=ISI) 18. 18.Fromer, M. et al. Gene expression elucidates functional impact of polygenic risk for schizophrenia. Nat. Neurosci. 19, 1442–1453 (2016). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nn.4399&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=27668389&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 19. 19.Gandal, M. J. et al. Transcriptome-wide isoform-level dysregulation in ASD, schizophrenia, and bipolar disorder. Science 362, (2018). 20. 20.Zheng, G. X. Y. et al. Massively parallel digital transcriptional profiling of single cells. Nat. Commun. 8, 14049 (2017). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/ncomms14049&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=28091601&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 21. 21.Aldridge, S. & Teichmann, S. A. Single cell transcriptomics comes of age. Nat. Commun. 11, 4307 (2020). 22. 22.Mathys, H. et al. Single-cell transcriptomic analysis of Alzheimer’s disease. Nature 570, 332–337 (2019). 23. 23.Velmeshev, D. et al. Single-cell genomics identifies cell type-specific molecular changes in autism. Science 364, 685–689 (2019). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6Mzoic2NpIjtzOjU6InJlc2lkIjtzOjEyOiIzNjQvNjQ0MS82ODUiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyMC8xMS8wOS8yMDIwLjExLjA2LjIwMjI1MzQyLmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 24. 24.Nagy, C. et al. Single-nucleus transcriptomics of the prefrontal cortex in major depressive disorder implicates oligodendrocyte precursor cells and excitatory neurons. Nat. Neurosci. (2020) doi:10.1038/s41593-020-0621-y. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41593-020-0621-y&link_type=DOI) 25. 25.Schafflick, D. et al. Integrated single cell analysis of blood and cerebrospinal fluid leukocytes in multiple sclerosis. Nat. Commun. 11, 247 (2020). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41467-019-14118-w&link_type=DOI) 26. 26.Wang, D. et al. Comprehensive functional genomic resource and integrative model for the human brain. Science 362, (2018). 27. 27.McGinnis, C. S. et al. MULTI-seq: sample multiplexing for single-cell RNA sequencing using lipid-tagged indices. Nat. Methods 16, 619–626 (2019). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41592-019-0433-8&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=31209384&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 28. 28.Mohammadi, S., Davila-Velderrain, J. & Kellis, M. A multiresolution framework to characterize single-cell state landscapes. doi:10.1101/746339. [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoiYmlvcnhpdiI7czo1OiJyZXNpZCI7czo4OiI3NDYzMzl2MiI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIwLzExLzA5LzIwMjAuMTEuMDYuMjAyMjUzNDIuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 29. 29.Boldog, E. et al. Transcriptomic and morphophysiological evidence for a specialized human cortical GABAergic cell type. Nat. Neurosci. 21, 1185–1195 (2018). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41593-018-0205-2&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=30150662&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 30. 30.Lim, L., Mi, D., Llorca, A. & Marín, O. Development and Functional Diversification of Cortical Interneurons. Neuron 100, 294–313 (2018). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/J.NEURON.2018.10.009&link_type=DOI) 31. 31.Demars, F. et al. Dysregulation of peripheral expression of the YWHA genes during conversion to psychosis. Sci. Rep. 10, 9863 (2020). 32. 32.Enwright, J.F., Iii et al. Transcriptome alterations of prefrontal cortical parvalbumin neurons in schizophrenia. Mol. Psychiatry 23, 1606–1613 (2018). 33. 33.Rannals, M. D. et al. Psychiatric Risk Gene Transcription Factor 4 Regulates Intrinsic Excitability of Prefrontal Neurons via Repression of SCN10a and KCNQ1. Neuron 90, 43–55 (2016). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.neuron.2016.02.021&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=26971948&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 34. 34.Doostparast Torshizi, A. et al. Deconvolution of transcriptional networks identifies TCF4 as a master regulator in schizophrenia. Sci Adv 5, eaau4139 (2019). [FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6MzoiUERGIjtzOjExOiJqb3VybmFsQ29kZSI7czo4OiJhZHZhbmNlcyI7czo1OiJyZXNpZCI7czoxMjoiNS85L2VhYXU0MTM5IjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjAvMTEvMDkvMjAyMC4xMS4wNi4yMDIyNTM0Mi5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 35. 35.Jaffe, A. E. et al. Mapping DNA methylation across development, genotype and schizophrenia in the human frontal cortex. Nat. Neurosci. 19, 40–47 (2016). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nn.4181&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=26619358&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 36. 36.Pietersen, C. Y. et al. Molecular profiles of pyramidal neurons in the superior temporal cortex in schizophrenia. J. Neurogenet. 28, 53–69 (2014). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.3109/01677063.2014.882918&link_type=DOI) 37. 37.Pietersen, C. Y. et al. Molecular profiles of parvalbumin-immunoreactive neurons in the superior temporal cortex in schizophrenia. J. Neurogenet. 28, 70–85 (2014). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.3109/01677063.2013.878339&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=24628518&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 38. 38.Südhof, T. C. Neuroligins and neurexins link synaptic function to cognitive disease. Nature 455, 903–911 (2008). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nature07456&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=18923512&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000260038300038&link_type=ISI) 39. 39.Guilmatre, A., Huguet, G., Delorme, R. & Bourgeron, T. The emerging role of SHANK genes in neuropsychiatric disorders. Dev. Neurobiol. 74, 113–122 (2014). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/dneu.22128&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=24124131&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 40. 40.Maynard, K. R. et al. dotdotdot: an automated approach to quantify multiplex single molecule fluorescent in situ hybridization (smFISH) images in complex tissues. Nucleic Acids Res. 48, e66 (2020). 41. 41.Berdenis van Berlekom, A. et al. Synapse Pathology in Schizophrenia: A Meta-analysis of Postsynaptic Elements in Postmortem Brain Studies. Schizophr. Bull. 46, 374–386 (2020). 42. 42.Birnbaum, R. & Weinberger, D. R. Genetic insights into the neurodevelopmental origins of schizophrenia. Nat. Rev. Neurosci. 18, 727–740 (2017). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nrn.2017.125&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=29070826&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 43. 43.Föcking, M. et al. Proteomic and genomic evidence implicates the postsynaptic density in schizophrenia. Mol. Psychiatry 20, 424–432 (2015). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/mp.2014.63&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25048004&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 44. 44.Kristiansen, L. V., Beneyto, M., Haroutunian, V. & Meador-Woodruff, J. H. Changes in NMDA receptor subunits and interacting PSD proteins in dorsolateral prefrontal and anterior cingulate cortex indicate abnormal regional expression in schizophrenia. Mol. Psychiatry 11, 737–47, 705 (2006). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/sj.mp.4001844&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=16702973&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000239309800006&link_type=ISI) 45. 45.Moon, M.-S. & Gomez, T. M. Balanced Vav2 GEF activity regulates neurite outgrowth and branching in vitro and in vivo. Molecular and Cellular Neuroscience vol. 44 118–128 (2010). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.mcn.2010.03.001&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=20298788&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 46. 46.Pilpel, Y. & Segal, M. Rapid WAVE dynamics in dendritic spines of cultured hippocampal neurons is mediated by actin polymerization. J. Neurochem. 95, 1401–1410 (2005). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1111/j.1471-4159.2005.03467.x&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=16190876&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000233170200017&link_type=ISI) 47. 47.Fawcett, J. P. et al. Nck adaptor proteins control the organization of neuronal circuits important for walking. Proc. Natl. Acad. Sci. U. S. A. 104, 20973–20978 (2007). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NDoicG5hcyI7czo1OiJyZXNpZCI7czoxMjoiMTA0LzUyLzIwOTczIjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjAvMTEvMDkvMjAyMC4xMS4wNi4yMDIyNTM0Mi5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 48. 48.Cullen, T. J. et al. Anomalies of asymmetry of pyramidal cell density and structure in dorsolateral prefrontal cortex in schizophrenia. Br. J. Psychiatry 188, 26–31 (2006). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MTA6ImJqcHJjcHN5Y2giO3M6NToicmVzaWQiO3M6ODoiMTg4LzEvMjYiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyMC8xMS8wOS8yMDIwLjExLjA2LjIwMjI1MzQyLmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 49. 49.Garey, L. J. et al. Reduced dendritic spine density on cerebral cortical pyramidal neurons in schizophrenia. Journal of Neurology, Neurosurgery & Psychiatry vol. 65 446–453 (1998). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NDoiam5ucCI7czo1OiJyZXNpZCI7czo4OiI2NS80LzQ0NiI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIwLzExLzA5LzIwMjAuMTEuMDYuMjAyMjUzNDIuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 50. 50.Wagstyl, K. et al. Multiple markers of cortical morphology reveal evidence of supragranular thinning in schizophrenia. Transl. Psychiatry 6, e780 (2016). 51. 51.Guet-McCreight, A., Skinner, F. K. & Topolnik, L. Common Principles in Functional Organization of VIP/Calretinin Cell-Driven Disinhibitory Circuits Across Cortical Areas. Front. Neural Circuits 14, 32 (2020). 52. 52.Sey, N. Y. A. et al. A computational tool (H-MAGMA) for improved prediction of brain-disorder risk genes by incorporating brain chromatin interaction profiles. Nat. Neurosci. 23, 583–593 (2020). 53. 53.Howard, D. M. et al. Genome-wide meta-analysis of depression identifies 102 independent variants and highlights the importance of the prefrontal brain regions. Nat. Neurosci. 22, 343–352 (2019). [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 54. 54.Stahl, E. A. et al. Genome-wide association study identifies 30 loci associated with bipolar disorder. Nat. Genet. 51, 793–803 (2019). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41588-019-0397-8&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=31043756&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 55. 55.Grove, J. et al. Identification of common genetic risk variants for autism spectrum disorder. Nat. Genet. 51, 431–444 (2019). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41588-019-0344-8&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=30804558&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 56. 56.Demontis, D. et al. Discovery of the first genome-wide significant risk loci for attention deficit/hyperactivity disorder. Nat. Genet. 51, 63–75 (2019). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41588-018-0269-7&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=30478444&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 57. 57.Jansen, I. E. et al. Genome-wide meta-analysis identifies new loci and functional pathways influencing Alzheimer’s disease risk. Nat. Genet. 51, 404–413 (2019). [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 58. 58.Cross-Disorder Group of the Psychiatric Genomics Consortium et al. Genetic relationship between five psychiatric disorders estimated from genome-wide SNPs. Nat. Genet. 45, 984–994 (2013). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/ng.2711&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=23933821&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 59. 59.Maynard, K. E., Collado-Torres, L., Weber, L. M. & Uytingco, C. Transcriptome-scale spatial gene expression in the human dorsolateral prefrontal cortex. bioRxiv (2020). 60. 60.Metzner, C., Zurowski, B. & Steuber, V. The Role of Parvalbumin-positive Interneurons in Auditory Steady-State Response Deficits in Schizophrenia. Sci. Rep. 9, 18525 (2019). 61. 61.Bartos, M., Vida, I. & Jonas, P. Synaptic mechanisms of synchronized gamma oscillations in inhibitory interneuron networks. Nat. Rev. Neurosci. 8, 45–56 (2007). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nrn2044&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=17180162&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000242994200015&link_type=ISI) 62. 62.Forrest, M. P. et al. The Psychiatric Risk Gene Transcription Factor 4 (TCF4) Regulates Neurodevelopmental Pathways Associated With Schizophrenia, Autism, and Intellectual Disability. Schizophr. Bull. 44, 1100–1110 (2018). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/schbul/sbx164&link_type=DOI) 63. 63.Larti, F. et al. A defect in the CLIP1 gene (CLIP-170) can cause autosomal recessive intellectual disability. Eur. J. Hum. Genet. 23, 416 (2015). [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25672248&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 64. 64.Liu, X., Li, Y. I. & Pritchard, J. K. Trans Effects on Gene Expression Can Drive Omnigenic Inheritance. Cell 177, 1022–1034.e6 (2019). 65. 65.Keenan, A. B. et al. ChEA3: transcription factor enrichment analysis by orthogonal omics integration. Nucleic Acids Res. 47, W212–W224 (2019). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/nar/gkz446&link_type=DOI) 66. 66.Lai, T. et al. SOX5 controls the sequential generation of distinct corticofugal neuron subtypes. Neuron 57, 232–247 (2008). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.neuron.2007.12.023&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=18215621&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000252758400008&link_type=ISI) 67. 67.Zahr, S. K., Kaplan, D. R. & Miller, F. D. Translating neural stem cells to neurons in the mammalian brain. Cell Death Differ. 26, 2495–2512 (2019). 68. 68.1. R. Reynolds Phillips, J. M. Neurodevelopmental and Genetic Disorders: A Book Review - Handbook of Neurodevelopmental and Genetic Disorders in Children. Second edition, Sam Goldstein, Cecil R. Reynolds (Eds.). (2011). New York: The Guilford Press, 588 pp., $75.00 (HB). Journal of the International Neuropsychological Society vol. 17 1167–1168 (2011). 69. 69.Coe, B. P. et al. Neurodevelopmental disease genes implicated by de novo mutation and copy number variation morbidity. Nat. Genet. 51, 106–116 (2019). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41588-018-0288-4&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=30559488&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 70. 70.Li, M. et al. Integrative functional genomic analysis of human brain development and neuropsychiatric risks. Science 362, (2018). 71. 71.Kaya-Okur, H. S. et al. CUT&Tag for efficient epigenomic profiling of small samples and single cells. Nat. Commun. 10, 1930 (2019). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41467-019-09982-5&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=31036827&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 72. 72.Harrington, A. J. et al. MEF2C regulates cortical inhibitory and excitatory synapses and behaviors relevant to neurodevelopmental disorders. Elife 5, (2016). 73. 73.Verret, L. et al. Inhibitory interneuron deficit links altered network activity and cognitive dysfunction in Alzheimer model. Cell 149, 708–721 (2012). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.cell.2012.02.046&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=22541439&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000303443100020&link_type=ISI) 74. 74.Gonzalez-Burgos, G., Cho, R. Y. & Lewis, D. A. Alterations in cortical network oscillations and parvalbumin neurons in schizophrenia. Biol. Psychiatry 77, 1031–1040 (2015). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.biopsych.2015.03.010&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25863358&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 75. 75.Nelson, M. R. et al. The support of human genetic evidence for approved drug indications. Nat. Genet. 47, 856–860 (2015). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/ng.3314&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=26121088&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 76. 76.Butler, A., Hoffman, P., Smibert, P., Papalexi, E. & Satija, R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat. Biotechnol. 36, 411–420 (2018). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nbt.4096&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=29608179&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 77. 77.Mohammadi, S., Ravindra, V., Gleich, D. F. & Grama, A. A geometric approach to characterize the functional identity of single cells. Nat. Commun. 9, 1516 (2018). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41467-018-03933-2&link_type=DOI) 78. 78.Batagelj, V. An O(m) Algorithm for Cores Decomposition of Networks. (2001). 79. 79.Crowell, H. L. et al. On the discovery of population-specific state transitions from multi-sample multicondition single-cell RNA sequencing data. bioRxiv 713412 (2019) doi:10.1101/713412. [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NzoiYmlvcnhpdiI7czo1OiJyZXNpZCI7czo4OiI3MTM0MTJ2MyI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIwLzExLzA5LzIwMjAuMTEuMDYuMjAyMjUzNDIuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 80. 80.Ritchie, M. E. et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43, e47 (2015). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/nar/gkv007&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25605792&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 81. 81.Law, C. W., Chen, Y., Shi, W. & Smyth, G. K. voom: Precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biol. 15, R29 (2014). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1186/gb-2014-15-2-r29&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=24485249&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 82. 82.Meers, M. P., Tenenbaum, D. & Henikoff, S. Peak calling by Sparse Enrichment Analysis for CUT&RUN chromatin profiling. Epigenetics Chromatin 12, 42 (2019). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1186/s13072-019-0287-4&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=31300027&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom) 83. 83.He, Z. et al. Comprehensive transcriptome analysis of neocortical layers in humans, chimpanzees and macaques. Nat. Neurosci. 20, 886–895 (2017). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nn.4548&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=28414332&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F11%2F09%2F2020.11.06.20225342.atom)