Transcriptomic pathology of neocortical microcircuit cell types across psychiatric disorders

Human subjects

We obtained brain samples during autopsies conducted at the Allegheny County Medical Examiner’s Office in Pittsburgh, PA, after we procured consent from the next of kin. An independent committee of clinical research scientists at the University of Pittsburgh made consensus DSM-IV diagnoses using clinical records, toxicology results, and standardized psychological autopsies. We included 76 subjects in our study, divided into 19 groups of 4 subjects (“tetrads”): one with major depressive disorder (MDD), one with bipolar disorder (BD), one with schizophrenia (SCZ), and one without any psychiatric disorder (control). The four subjects in each tetrad were matched for sex and as closely as possible for age. Groups did not differ in mean postmortem interval (PMI), pH, or RIN. We sectioned fresh-frozen tissue blocks containing the subgenual anterior cingulate cortex (sgACC) at 12 μm thickness, mounted these on RNase-free laser capture microdissection (LCM)-compatible slides, and stored them at −80°C until we harvested the cells. All procedures were approved by the University of Pittsburgh Committee for Oversight of Research and Clinical Training Involving Decedents and Institutional Review Board for Biomedical Research. We describe the demographic and clinical characteristics of our cohort in Supplemental Table 1.

Cell density analysis

We performed double-label fluorescent in situ hybridization (FISH) on sgACC sections, targeting Slc17a7 (a marker of pyramidal cells) + Sst, or Pvalb + Vip, counterstained with DAPI to determine changes in cell type density. We performed RNAscope staining in accordance with the manufacturer’s instructions. We adopted a stereological approach to randomly sample 20 sites of the sgACC for imaging with an Olympus IX3 confocal microscope (Olympus, Tokyo, Japan). We distinguished L2/3 and L5/6 pyramidal (PYR) cells through the selection of sites superficial and deep to a reference line midway between the pial surface and white matter. We obtained Z-stacks of images that represented transverse depth into tissue sections at a step size of 0.25 μm. We imaged each site to a depth of 12 Z-stacks, reflecting the two boundaries of the tissue beyond which our microscope focus could not be achieved. We deconvolved images using the AutoQuant Blind algorithm in Slidebook version 6.0 (Intelligent Imaging Innovations, Denver, CO) to reduce background; we collapsed Z-stacks into a composite image, normalizing intensities across all samples. We manually quantified cells containing 10 or more cell-marker grains and overlapping with DAPI. We adjusted cell type densities for age, sex, and PMI (Supplemental Figure 1).

Fluorescent labeling and capture of microcircuit neuronal subtypes

We used a rapid FISH protocol, a shortened version of the RNAscope Multiplex Fluorescent V2 kit protocol (Advanced Cell Diagnostics, Newark, CA), to stain sgACC sections with cDNA probes specific to Slc17a7, Sst, Pvalb, or Vip, thereby visualizing microcircuit cell types, as previously described (39,40). Just as in our cell density analyses, we identified L2/3 and L5/6 PYR-cells using Slc17a7 and distinguished them with a reference line placed halfway between the pial surface and white matter. We used one probe per tissue section and made all probes visible with the Opal 520 dye (FITC range) to minimize detection bias across cell types. We counter-stained sections with Neurotrace Red (ThermoFisher, Waltham, MA), a fluorescent Nissl stain, to visualize cell bodies. We then collected 130 of each cell type from each subject using an LMD7 laser microdissection microscope (Leica Microsystems). Altogether, we dissected and captured 49,400 neurons across all samples.

Library preparation and RNA-sequencing

We extracted RNA using a PicoPure RNA isolation kit (Thermo Fisher Scientific, Waltham, MA), and prepared libraries with the SMARTer Stranded Total RNA-Seq Pico Input Kit V2 (Clontech, Mountain View, CA). We determined library fragment size distribution with a Bioanalyzer (Agilent Technologies, Santa Clara, CA), and carried out initial quality control sequencing on an Illumina MiSeq sequencer (Illumina, San Diego, CA). Given the limited RNA quantity from each pool of our microdissected neurons, we did not assess RNA integrity number (RIN) values in the neuron samples used for RNA-seq profiling. However, we have previously shown that our Nissl staining and laser microdissection method results in RIN >7 in all samples, and values are nearly identical to those obtained in tissue homogenates from the same subjects (35). If libraries showed < 80% reads aligning to the genome, or < 15% of reads aligning to exons, we re-collected them (i.e., repeated the RNAscope-LCM procedure using new tissue sections) and sequenced again. We loaded 3 or 4 tetrads on individual flow cells, using 6 flow-cells in total: 5 for initial sequencing and 1 for re-collected cells. Libraries that passed quality control (N = 379/380, 1 subject had insufficient gray matter for PV-cell collection) were sequenced on an Illumina NovaSeq sequencer (Illumina, San Diego, CA) to generate 2 x 125 bp paired-end reads. We pooled 3 or 4 tetrads and loaded onto each flow cell, using 5 NovaSeq S2 flow-cells.

We aligned paired-end reads to GRCh38 using the Ensembl annotation (Release 98) with alignment performed using STAR v2.7.2 with default parameters (41). We used the summarizeOverlaps function from the GenomicAlignments R package to quantify genes, with separate quantifications for exon- and intron-aligning reads (42). Before gene quantification, we utilized the ScanBamParam function from Rsamtools to extract read information, exclude PCR/optical duplicates, and filter for uniquely aligning reads that mapped to their proper read pair. NovaSeq sequencing yielded an average of 48,457,414 reads per cell type per subject (SE ± 377,012), with an average genomic mapping rate of 78.39% (SE ± 0.41%; 38,124,400 reads ± 388,752), exonic alignment rate of 31.48% (SE ± 0.39%), and intronic alignment rate of 49.26% (SE ± 0.15%). We combined exonic and intronic reads and removed mitochondrially encoded genes (genes with the MT- prefix) for downstream analysis steps.

Differential expression analysis

We performed gene filtering and differential expression (DE) analysis separately for each neuronal cell type. We included genes with > 1 count in ≥ 80% of biological replicates. After normalizing counts and calculating precision weights using voom, we computed the variance of gene expression using linear models implemented in the limma package (43). The main independent variable for our models was the diagnostic group (MDD, BD, SCZ versus control), and we also adjusted for sex and intergenic rate.

In addition to diagnostic group, we considered whether to include each of 18 biological and technical covariates in our models. Biological covariates comprised sex, age, race, method of death (MOD), tobacco use, medication use at time of death (antipsychotics, antidepressants, lithium, benzodiazepines and/or anticonvulsants), while technical covariates comprised pH, RIN, PMI, pool (batch), number of detected genes, library size, intronic rate, intergenic rate, and intragenic rate. Many of these covariates were highly correlated with each other (Supplemental Figure 2A) and were correlated with principal components (PCs) of the normalized gene expression data (Supplemental Figure 2B). We selected the final covariates empirically within each cell type using the Bayesian Information Criterion (BIC) (44) with the mvForwardStepwise function from the mvIC R package. One phase of covariate selection was performed on biological covariates, and a second phase, including significant covariates from the first phase, was performed on technical covariates. Of the 18 covariates, the ones selected were sex (all five cell types), intergenic rate (all cell types except L2/3 PYR), race (PVALB and L5/6 PYR only), and MOD (L2/3 PYR and L5/6 PYR only). We included only sex and intergenic rate in the final linear models to ensure a consistent interpretation of DE across cell types. Importantly, medication use covariates were not significantly correlated with gene expression PCs or selected by mvForwardStepwise.

We considered genes with a Benjamini-Hochberg adjusted p-value (FDR) < 0.20 to be DE. This exploratory threshold enabled us to assess the overlap between cell types and disorders. We also compared our DE genes to those from a high-quality meta-analysis of bulk cortex microarray and RNA-seq data from cerebral cortex (45). Additionally, we evaluated the specificity of cell type-specific DE genes using an external single nucleus RNA-seq dataset of human medial temporal gyrus from the Allen Institute (46). First, we mapped each of the five LCM-seq neuronal subtypes to its closest snRNA-seq ‘subclass’ based on maximal expression of subclass-specific marker genes. We considered DE genes cell type-specific if at least 10% of cells of the respective snRNA-seq subclass expressed that gene (Supplemental Figure 3B). We also calculated the gene-wise correlation of expression levels between LCM-seq controls and the snRNA-seq reference to assess concordance between these sequencing methods (Supplemental Figure 3C).

We further compared the concordance between LCM-seq gene signatures and recent snRNA-seq datasets of MDD and SCZ (30,34). We obtained snRNA-seq count matrices, filtered out low-quality cells and genes, and performed pseudobulk DE testing using limma-voom (43). LCM-seq cell types were manually matched to the authors’ original cell type annotations based on semantic similarities (e.g., InN9_PV mapped to PVALB) and correlation of DE t-statistics across genes. To assess the overlap of gene signatures between LCM-seq and snRNA-seq in specific cell types, we employed the rank-rank hypergeometric overlap test (RRHO) (47). Genes were ranked according to their p-value and effect size direction, allowing us to identify significantly overlapping genes across a continuous range of significance levels (Supplemental Figure 4).

Differential pathway activity analysis

We curated disease-associated pathways (gene sets) by filtering Gene Ontology (GO) biological processes (downloaded from the Mayaan laboratory; https://maayanlab.cloud/Enrichr) to pathways whose names contained manually curated keywords related to biological processes known to be affected in brain disorders, including bioenergetics, metabolism, transcription and translation, growth factors, cytoskeleton, cellular trafficking, neuron structure and function, and immune system and inflammation (see Supplemental Table 2 for the list of keywords). We manually curated lists to ensure the alignment of the pathways with the chosen biological themes. We performed differential pathway activity analysis similarly to methods published previously (48). We used Gene Set Variation Analysis (GSVA) (49) to summarize normalized gene expression into activity scores for each pathway, set to the following parameters: mx.diff = TRUE, kcdf = “Gaussian”, min.sz = 5, max.sz = 150. Briefly, GSVA calculates a normalized relative expression level per gene across samples. This expression level is then rank-ordered for each sample and aggregated into a pathway-level activity score by calculating sample-wise enrichment scores using a Kolmogorov–Smirnov-like rank statistic. To minimize false positives, we applied certain filters: genes that were not expressed in a given cell type (<1 CPM in ≥70% of samples) were excluded, as well as pathways for which >30% of the genes in the set were not expressed. For each cell type and pathway, we modeled pathway activity as a function of diagnosis, sex, and intergenic rate using lmFit in limma (43). Pathways demonstrating evidence of association with diagnosis at p < 0.01 in at least one cell type and disorder were included in the final analyses.

Next, we identified coordinated alterations in biological pathways within the cortical microcircuit associated with each disorder. To achieve this, we examined the DE t-statistic of all genes defining the pathways within each biological theme and correlated them across cell types pairs. A significant (p < 0.05) positive Pearson correlation between two cell types was interpreted as a coordinated biological change for that disease.

Overlap with genetic risk loci

We evaluated the intersection of cell type-specific DE genes and psychiatric disorder-associated risk genes using PoPS and PsyOPS, advanced gene-prioritization methods for GWAS (61,62). PoPS integrates multiple gene features, such as cell-type-specific expression and biological pathways, using aggregated gene-level p-values from GWAS summary statistics. PsyOPS is trained on three features that hypothesize causal genes are mutationally constrained, brain-specific, and overlap with known neurodevelopmental genes. For this analysis, we provided both tools with GWAS summary statistics from recent large-scale studies on MDD, BD, and SCZ (50–52). After running the tools as prescribed, we merged their outputs with our DE gene list (FDR < 0.20). The final dataset included only GWAS loci that were highlighted by at least one tool and overlapped with a disorder-specific DE gene (Figure 3B).

LCM-seq generates robust cell type-specific transcriptomes from cortical microcircuit neuronal subtypes in postmortem human brain tissues

We performed laser capture microdissection (LCM) followed by RNA sequencing (LCM-seq) on postmortem samples from the subgenual anterior cingulate cortex (sgACC) of 76 subjects. These subjects comprised 19 groups of 4 subjects (“tetrads”) matched closely on age, sex, and tissue quality: one with major depressive disorder (MDD), one with bipolar disorder (BD), one with schizophrenia (SCZ), and one with no psychiatric disorder (control; Figure 1A, Table 1). We achieved cell type specificity using RNAscope-based FISH probes to identify each of five neuronal subtypes: Layer 2/3 and layer 5/6 glutamatergic pyramidal neurons (L2/3, L5/6 PYR), as well as vasoactive intestinal peptide- (VIP), somatostatin- (SST), and parvalbumin- (PVALB) expressing inhibitory interneurons. We then used LCM to collect 130 cells for each neuron type per subject prior to performing RNA sequencing. Thus, we generated 76 (subjects) × 5 (cell types) = 380 bulk transcriptomes, each containing the pooled RNA from 130 cells.

In UMAP space, we observed a clustering of transcriptomes primarily by cell type (93.8% accuracy, k-means clustering, k = 5), with some minor overlap between L2/3 and L5/6 PYR cells (Figure 1B). As expected, each cell type had enriched expression of cell type-specific marker genes (PVALB, SST, VIP, CUX2 = L2/3 PYR, FEZF2 = L5/6 PYR, all p < 1.0 × 10⁻⁴) (Figure 1C). Further demonstrating the cell type-specificity of the data, gene expression levels from LCM-seq profiled cell types exhibited a strong correlation with those from a single-nucleus human neocortical RNA-seq reference atlas (46) (r = 0.85 ± 0.0021, all p < 2.6 × 10⁻¹⁶) (Supplemental Figure 3B). In parallel tissue sections from those used for cell type collection, we determined cell densities using FISH. No cell types showed significant differences in microscopy-based cellular density between cases and controls after adjusting for sex, age and PMI, although trending decreases in interneuron densities were observed for SCZ (Supplemental Figure 1).

Differential expression analysis highlights distinct transcriptional alterations in inhibitory interneuron subtypes, partially shared across psychiatric disorders

To identify disease-associated transcriptional alterations, we performed case-control differential expression (DE) testing separately for each cell type. A total of 460 DE genes (N = 371 unique) were identified across cell types and disorders at FDR < 0.20, with the majority (N = 398, 87.3%) observed in inhibitory interneurons (Figure 1D). PVALB cells displayed the greatest number of DE genes (N = 239), with 91 for SCZ (25 upregulated, 66 downregulated), 33 for BD (5 up, 28 down), and 113 for MDD (28 up, 85 down). For SST cells, there were 38 DE genes for SCZ (31 up, 7 down), 40 for BD (27 up, 13 down), and 12 for MDD (6 up, 6 down). For VIP cells, there were 24 DE genes for SCZ (7 up, 17 down), 11 for BD (6 up, 5 down), and 36 for MDD (5 up, 31 down). The majority of DE in excitatory neurons was observed in L5/6 PYR for SCZ (N = 41, 13 up, 28 down). See Supplemental Table 1 for the complete DE analysis results.

Figure 1E presents the top altered genes across neuronal subtypes in each disorder, ordered by significance (FDR). Although, the majority of DE genes (N = 314, 84.6%) were perturbed in a single neuronal subtype, 9.2% (N = 34) of DE genes were altered in the same direction and cell type across at least two disorders (p = 2.6 × 10⁻³³, hypergeometric test). In other words, disease-associated expression changes were highly cell type-specific; however, these changes were shared in part across MDD, BD, and SCZ.

Nearly all cell type-specific DE genes shared across disorders are involved in aspects of the formation and maintenance of neuronal circuits. For example, NRXN3 (coding for neurexin 3) was downregulated in PVALB cells in BD (log₂ fold change [logFC] = −0.38, FDR = 0.03) and MDD (logFC = −0.56, FDR = 9.0 × 10⁻¹²), and downregulated in VIP cells in SCZ (logFC = 0.33, FDR = 1.4 × 10⁻⁴) and MDD (logFC = −0.23, FDR = 3.6 × 10⁻¹⁰). NRXN3 is part of a family of presynaptic single-pass transmembrane proteins that act as synaptic organizers, and genomic alterations in NRXN3 have been identified in a number of neuropsychiatric disorders, including autism spectrum disorders (ASD), schizophrenia, intellectual disability (ID), and addiction (53). CNTNAP2 (coding for Contactin associated Protein 2) also belongs to the neurexin family and is required for the organization of myelinated axons, and was downregulated in VIP cells in SCZ (logFC = −0.42, FDR = 5.1 × 10⁻⁴) and MDD (logFC = −0.27, FDR = 0.18). Heterozygous missense and loss-of-function (LOF) CNTNAP2 mutations cause several neurodevelopmental disorders, including ASD and ID (54), likely through the disruption of synaptic transmission, neuronal connectivity, and network-level activity as identified in LOF mouse models (55,56). As another example, ERBB4 (coding for Erb-B2 Receptor Tyrosine Kinase 4) was downregulated in VIP cells in SCZ (logFC = −0.56, FDR = 5.1 × 10⁻¹²) and MDD (logFC = −0.36, FDR = 0.014), and NRG3 (coding for Neuregulin 3) was downregulated in SST cells in BD (logFC = −0.20, FDR = 5.0 × 10⁻²⁰) and MDD (logFC = −0.23, FDR = 1.2 × 10⁻¹⁵). ERBB4 codes for different isoforms of the receptor tyrosine-protein kinase, which binds its ligand, NRG3, a trophic factor. Together, ERBB4 and NRG3 play a critical role in the establishment and function of neuronal circuits (57), and their disruption increases susceptibility to SCZ and BD (58). Other examples of affected genes regulating neuronal circuits are NPAS3, FGF14 and AUTS2. Genes regulating neuronal signaling, including DLG2, GRIK2 and KCNIP4, were also DE across disorders. See Supplemental Figure 5 for the full list of shared genes.

We compared our cell type-specific patterns of differential expression to those from a prior meta-analysis of microarray and RNA-seq gene expression in bulk neocortical samples (45) (Figure 1E). A total of 43 cell type-specific DE genes from our analyses were convergent (same logFC sign, FDR < 0.20) with DE changes in bulk tissue (p = 7.08 × 10⁻⁵⁷, hypergeometric test), with the most overlaps in PVALB cells in SCZ (N = 13). As one example, bulk RNA-seq identified the upregulation of NPAS3 in SCZ (logFC = 0.14, p = 5.5 × 10⁻⁴), a change that was uniquely observed in SST cells in our analyses (logFC = 0.36, p = 7.4 × 10⁻⁶). NPAS3 plays a broad role in neurogenesis, and mutations in this gene have been associated with SCZ and neurodevelopmental disorders (59,60). LCM-seq cell-specific disease signatures were also moderately consistent with those surveyed in recent human snRNA-seq datasets from MDD and SCZ (30,34), as assessed by rank-rank hypergeometric overlap analysis (Supplemental Figure 4, see Methods).

Neuronal subtypes exhibit distinct and coordinated dysregulation of biological processes across psychiatric disorders

To characterize the biological effects of transcriptomic dysregulation in neuronal subtypes, we performed differential pathway activity analysis. We focused on 1881 pathways representing eight broad biological themes, including bioenergetics, neuron structure and function, and inflammation, known a priori to be relevant to psychiatric disorders. For each cell type and pathway, we defined an “activity score” summarizing the collective expression levels of genes in the pathway for that cell type, and diagnosis-specific differences were calculated using multivariate linear models (see Methods).

Figure 2A presents all significantly altered pathways (N = 276, p < 0.01), organized by biological theme. The greatest number of perturbed pathways were in VIP (N = 115), SST (N = 93), and L5/6 PYR (N = 100) cells in SCZ. Although alterations varied between cell types, on average, 11.3% of significant pathway changes within an individual cell type were also significantly altered in the same direction in another disorder. The highest overlap of significant pathways between disorders was observed between MDD and SCZ in VIP cells (31.6%, p = 0.026, hypergeometric test).

• Download figure
• Open in new tab

Figure 2. Cortical neuron subtypes exhibit disruption of biological pathways across psychiatric disorders.

(A, B) Differential pathway activity analysis results. A) Heatmap of disease-associated GO biological pathways (rows) altered in various cell type-disorder contrasts (columns), grouped by biological theme. Red and blue colors represent the significance (signed −log10 p-value) of up- and down-regulated pathway activities, respectively. See Supplemental Table 2 for the full list of pathways. B) Differential pathway activity (x-axis, signed −log10 p-value) of select biological processes in each cell type (rows) and disorder (shape) (dashed line = p < 0.05). C) Similarity of biological pathway alterations between neuronal subtypes across disorders, indicated by the significance (line width, −log10 p-value) of the Pearson correlation of genes in various biological themes (same as A,B) between cell types and across disorders (edge color). All correlations are positive. See Supplemental Figure 6 for all eight biological themes.

Figure 2B highlights cell type-specific differential activities for pathways that reflect known and novel aspects of cellular pathology in psychiatric disorders. For SST cells, this included the downregulation of genes involved in ATPase activity in SCZ (p = 5.5 × 10⁻³), glial-cell derived neurotrophic factor receptor signaling in SCZ and BD (p = 6.0 × 10⁻³, p = 2.3 × 10⁻³), lysosomal protein catabolic process (p = 2.3 × 10⁻³), and retinoic acid receptor signaling in SCZ and MDD (p = 0.025, p = 1.8 × 10⁻³). Notably, there was an upregulation of genes involved in response to endoplasmic reticulum (ER) stress in SST cells for all three disorders (SCZ: p = 0.034, BD: p = 0.02, MDD: p = 0.023), corresponding to cell type-specific findings in a chronic psychosocial stress model of depression in rodents (25). For VIP cells, there was increased negative regulation of genes involved in cellular amide metabolism in BD and MDD (p = 8.2 × 10⁻³, p = 8.5 × 10⁻³), translation in SCZ and BD (p = 0.011, p = 5.0 × 10⁻⁴), and cell-matrix adhesion mediator activity in SCZ and BD (p = 1.9 10⁻³, p = 6.1 × 10⁻³). For PVALB cells, there was an upregulation of genes involved in long-term synaptic depression in BD and MDD (p = 7.1 × 10⁻³, p = 0.029) and a downregulation of genes involved in oxidoreductase activity in BD and MDD (p = 4.8 × 10⁻⁴, p = 5.6 × 10⁻³). A full list organized per cell type, themes, and biological pathways is available in Supplemental Table 2.

We next assessed the degree of shared alterations in biological pathways between cortical microcircuit neuronal cell types within each disorder. For each of the eight biological themes, we tabulated the set of genes found in at least one pathway associated with that theme. Then, for each theme and disorder, we calculated the Pearson correlation of DE t-statistics for these genes between each pair of neuronal subtypes. There were no significant negative correlations; only coordinated changes were observed (Figure 2C). Several findings emerge from this analysis. First, L2/3 PYR–L5/6 PYR and SST–PV cell pairs showed frequent coordinated changes for 6/8 and 5/8 biological themes, respectively (Supplemental Figure 6). Second, genes related to transcription and translation, cellular trafficking, and neuron structure and function were correlated between L2/3 PYR and L5/6 PYR cells for all three disorders, as were growth factor-related genes between SST and PVALB cells. Third, genes in neuron structure and function pathways had the most pairwise correlations between cell types, indicating coordinated changes at the microcircuit level. In contrast, other pairs of cells showed much fewer correlated gene changes across disorders, such as VIP–PV, VIP–PYRs, SST–PYRs or SST–PVALB.

Cell type-specific differentially expressed genes overlap with known genetic psychiatric risk genes

We asked whether cell type-specific DE genes overlapped with risk genes identified by genome-wide association studies (GWAS) of psychiatric disorders. To determine candidate causal genes that mediate the effects of disease-associated genetic variants from recent large-scale GWAS of MDD, BD, and SCZ (51,52,61) we utilized PsyOPS and PoPS, two methods that prioritize casual genes at GWAS lo ci (62,63) (see Methods). Our analysis revealed 22 candidate causal genes exhibiting DE within one or more neuronal subtypes – almost all of which were expressed in interneurons rather than in PYR cells (Figure 3B). Candidate causal genes were especially enriched for DE effects in PVALB cells in SCZ and MDD (p = 0.015, p = 0.026; permutation test).

• Download figure
• Open in new tab

Figure 3. Cell type-specific differentially expressed genes overlap with known genetic risk variants from psychiatric genome-wide association studies.

A) Overview of approach. Two recent methods, PsyOPS and PoPS (62,63), are used to predict genes that are causally associated with psychiatric GWAS loci, which are then assessed for overlap with LCM-seq differential expression (DE). B) 29 unique single-nucleotide polymorphisms (SNPs) associated with MDD, BD or SCZ (diagnosis; Dx) were within 500 kb of a cell type-specific DE gene (FDR < 0.20) from that respective disorder and prioritized as causal by either PsyOPS and/or PoPS. Overlap of SNPs and DE genes was significant for PVALB interneurons in SCZ (p = 0.015) and MDD (p = 0.026, permutation test). Multiple SNPs can be prioritized to the same gene (red text), and three genes (RBFOX1, OPCML, CACNA1C) were shared between two disorders (bold text).

For example, we found that the lead variant (i.e., the most significant GWAS variant within the link disequilibrium block) rs12447542, associated with SCZ, and rs7193263 and rs11077206, two independent lead variants associated with MDD, were predicted by PoPS and PsyOPS to casually effect the RBFOX1 gene (Figure 3B). RBFOX1 codes for a neuron-specific RNA-binding protein and splicing factor, and the RBFOX1 locus is highly pleiotropic, being associated with many psychiatric diseases (64,65). We found that RBFOX1 is upregulated uniquely in VIP cells – and no other cell type – across all three disorders (SCZ logFC = 0.41, FDR = 4.4 × 10⁻⁷; BD logFC = 0.53, FDR = 2.97 × 10⁻⁶; MDD logFC = 0.22, FDR = 0.025), supporting a role for VIP cells as a key mediator of RBFOX1’s effects on diverse psychiatric disorders. However, it is worth noting that RBFOX1 has been shown through experimental studies to also regulate splicing and neuronal function in other types of interneurons, including PV and SST (66,67). Similarly, rs2053079, associated with SCZ, is predicted to have ZNF563 as the causal gene. ZNF563 is a zinc finger nuclease that regulates neuronal differentiation in the developing brain (68), and we observe that this risk gene is specifically upregulated in PVALB cells for SCZ (logFC = 0.61, FDR = 4.45 × 10⁻³). As a third example, rs7758630, associated with MDD, is predicted to have GRIK2 as the causal gene. GRIK2, which encodes a subunit of the kainate glutamate receptor, was upregulated in SST cells in SCZ (logFC = 0.61, FDR = 6.4 × 10^{− 23}), BD (logFC = 0.61, FDR = 2.8 × 10⁻²³), and MDD (logFC = 0.32, FDR = 1.4 × 10⁻¹²) (Fig 1F).