Transcriptome-wide isoform-level dysregulation in ASD, schizophrenia, and bipolar disorder

Abstract

Most genetic risk for psychiatric disease lies in regulatory regions, implicating pathogenic dysregulation of gene expression and splicing. However, comprehensive assessments of transcriptomic organization in diseased brains are limited. In this work, we integrated genotypes and RNA sequencing in brain samples from 1695 individuals with autism spectrum disorder (ASD), schizophrenia, and bipolar disorder, as well as controls. More than 25% of the transcriptome exhibits differential splicing or expression, with isoform-level changes capturing the largest disease effects and genetic enrichments. Coexpression networks isolate disease-specific neuronal alterations, as well as microglial, astrocyte, and interferon-response modules defining previously unidentified neural-immune mechanisms. We integrated genetic and genomic data to perform a transcriptome-wide association study, prioritizing disease loci likely mediated by cis effects on brain expression. This transcriptome-wide characterization of the molecular pathology across three major psychiatric disorders provides a comprehensive resource for mechanistic insight and therapeutic development.

Figure 1.. Gene and isoform expression dysregulation in psychiatric brain

A) Differential expression effect size (|log₂FC|) histograms are shown for protein-coding, lncRNA, and pseudogene biotypes up or downregulated (FDR<0.05) in disease. Isoform-level changes (DTE; blue) show larger effect sizes than at gene level (DGE; red), particularly for protein-coding biotypes in ASD and SCZ. B) A literature-based comparison shows that the number of DE genes detected is dependent on study sample size for each disorder. C) Venn diagrams depict overlap among up or downregulated genes and isoforms across disorders. D) Gene ontology enrichments are shown for differentially expressed genes or isoforms. The top 5 pathways are shown for each disorder. E) Heatmap depicting cell type specificity of enrichment signals. Differentially expressed features show substantial enrichment for known CNS cell type markers, defined at the gene level from single cell RNA-Seq. F) Annotation of 944 unique non-coding RNAs DE in at least one disorder. From left to right: Sequence-based characterization of ncRNAs for measures of human selective constraint; brain developmental expression trajectories are similar across each disorder (colored lines represent mean trajectory across disorders); tissue, and CNS cell type expression patterns.

Figure 2.. Aberrant local splicing and isoform usage in ASD, SCZ and BD

A) Venn diagram showing cross-disorder overlap for 472 genes with significant differentially spliced (DS) intron clusters (FDR< 10%) identified by LeafCutter. P values for hypergeometric tests of pairwise overlaps between each disorder are shown at the bottom. B) Scatter plots comparing percent spliced-in (PSI) changes for all 1,287 introns in 515 significant DS clusters in at least one disorder, for significant disease pairs SCZ vs ASD and SCZ vs BD (Spearman’s ⍴=0.52 and ⍴=0.59, respectively). Principal component regression lines are shown in red, with regressions slopes for ASD and BD delta PSI compared to SCZ in the top-left corner. C) Top 10 gene ontology (GO) enrichments for DS genes in each disorder (see also Fig S8C). D) Significant enrichment for neuronal and astrocyte markers (ASD and SCZ), as well as oligodendrocyte and microglia (SCZ) cell type markers in DS genes. *Odds Ratio (OR) is given only for FDR< 5% and OR> 1. Oligo - oligodendrocytes; OPC - oligodendrocyte progenitor cells. E) A significant DS intron cluster in GRIN1 (clu_35560; chr9:140,040,354–140,043,461) showing increased exon 4 (E4) skipping in both ASD and SCZ. Increased or decreased intron usage in ASD/SCZ cases compared to controls are highlighted in red and blue, respectively. Protein domains are annotated as ANF_receptor - Extracellular receptor family ligand binding domain; Lig_chan - Ionotropic glutamate receptor; Lig_chan-Glu_bd - Ligated ion channel L-glutamate- and glycine-binding site; CaM_bdg_C0 - Calmodulin-binding domain C0 of NMDA receptor NR1 subunit. Visualization of splicing events in cluster clu_35560 with the change in PSI (ΔPSI) for ASD (left) and SCZ (right) group comparisons. FDR-corrected p-values (q) are indicated for each comparison. Covariate-adjusted average PSI levels in ASD or SCZ (red) vs CTL (blue) are indicated at each intron. F) Violin-plots with the distribution of covariate-adjusted PSI per sample for the intron skipping E4 are shown for each disease group comparison. G) DGE for GRIN1 in each disorder (*FDR< 5%). H) Whole-gene view of NRXN1 highlighting (dashed lines) the intron cluster with significant DS in ASD (clu_28264; chr2:50,847,321–50,850,452), as well as transcripts NRXN1–004 and NRXN1–012 that show significant DTU in SCZ and/or BD. Protein domain mappings are shown in purple. DM - Protein domains; Tx - Transcripts. ConA-like_dom_sf - Concanavalin A-like lectin/glucanase domain. EGF-like - Epidermal growth factor-like domain; Laminin_G - Laminin G domain; Neurexin-like - Neurexin/syndecan/glycophorin C domain. I) Left: close-up of exons and protein domains mapped onto the DS cluster, and FDR-corrected p-value (q). Right: visualization of introns in cluster clu_28264 with their change in percent spliced in (ΔPSI). Covariate-adjusted average PSI levels in ASD (red) vs CTL (blue) are indicated for each intron. J) Violin-plots with the distribution of covariate-adjusted PSI per sample for the largest intron skipping exon 8 (E8). K) Bar plots for changes in gene expression and transcript usage for NRXN1–004 and NRXN1–012 (*FDR< 5%).

Figure 3.. Overlaps and genetic enrichment among dysregulated transcriptomic features

A) Scatterplots demonstrate overlap among dysregulated transcriptomic features, summarized by their first principle component across subjects (R² values; *P<0.05). Polygenic risk scores (PRS) show greatest association with differential transcript signal in SCZ. B) SNP-heritability in SCZ is enriched among multiple differentially expressed transcriptomic features, with downregulated isoforms showing most substantial association via stratified LD-score regression. C) Several individual genes and isoforms exhibit genome-wide significant associations with disease PRS. Plots are split by direction of association with increasing PRS. In ASD, most associations localize to the 17q21.31 locus, harboring a common inversion polymorphism. In SCZ, a significant association as observed with C4A in the MHC locus.

Figure 4.. Transcriptome-wide association

Results from TWAS prioritize genes whose cis-regulated expression in brain is associated with disease. Plots show conditionally-independent TWAS prioritized genes, with lighter shade depicting marginal associations. The sign of TWAS Z-scores indicates predicted direction of effect. Genes significantly up or downregulated in disease brain are shown with arrows, indicating directionality. A) In SCZ, 193 genes (164 outside of MHC) are prioritized by TWAS at Bonferroni-corrected P<0.05, including 107 genes with conditionally independent signals. Of these, 23 are also differentially expressed in SCZ brain with 11 in the same direction as predicted. B) Seventeen genes are prioritized by TWAS in BD, of which 15 are conditionally independent. C) In ASD, TWAS prioritizes 12 genes, of which 5 are conditionally independent.

Figure 5.. Gene and isoform co-expression networks capture shared and disease-specific cellular processes and interactions

A) Gene and isoform co-expression networks demonstrate pervasive dysregulation across psychiatric disorders. Hierarchical clustering shows that separate gene- and isoform-based networks are highly overlapping, with greater specificity conferred at the isoform level. Disease associations are shown for each module (linear regression β value, * FDR<0.05, – P<0.05). Module cell type enrichments (*FDR<0.05) are shown for major CNS cell types defined from PsychENCODE UMI single cell clusters. Enrichments are shown for GWAS results from SCZ (59), BD (90), and ASD (38), using stratified LD score regression (* FDR<0.05, – P<0.05). B) Co-expression modules capture specific cellular identities and biological pathways. Colored circles represent module differential expression effect size in disease, with red outline representing GWAS enrichment in that disorder. Modules are organized and labeled based on CNS cell type and top gene ontology enrichments. C) Examples of specific modules dysregulated across disorders, with top 25 hub genes shown. Edges represent co-expression (Pearson correlation > 0.5) and known protein-protein interactions. Nodes are colored to represent disorders in which that gene is differentially expressed (*FDR<0.05).

Figure 6.. Two RBFOX1 isoform modules capture distinct biological and disease associations.

A) Previous studies have identified RBFOX1 as a critical hub of neuronal and synaptic modules downregulated across multiple psychiatric disorders (1, 16, 19, 32). Here, we identify two pairs of modules with distinct RBFOX1 isoforms as hub genes. Plots show the top 25 hub genes of modules isoM2 and isoM17, following the same coloring scheme as Fig 5C. B) Distinct module-eigengene trait associations are observed for isoM2 (downregulated in ASD only) compared with isoM17, which is downregulated in ASD and SCZ. C) Modules show distinct enrichments for nuclear and cytoplasmic RBFOX1 targets, defined experimentally in mouse (32). D) Genes harboring differential splicing events observed in ASD and SCZ show greater overlap with isoM17, consistent with its association with nuclear RBFOX1 targets. E) Modules show distinct patterns of genetic association. isoM2 exhibits broad enrichment for GWAS signal in SCZ, BD, and MDD, as well as for epilepsy risk genes, whereas isoM17 shows no apparent genetic enrichment. GWAS enrichments show FDR-corrected P-values calculated using stratified-LDSC, and rare-variant associations were calculated using logistic regression, controlling for gene length and GC content (21).

Figure 7.. Distinct neural-immune trajectories in disease

A) Co-expression networks provide substantial refinement of the neuro-immune/inflammatory processes upregulated in ASD, SCZ, and BD. Previous work has identified specific contributions to this signal from astrocyte and microglial populations (1, 19). Here, we further identify additional, distinct interferon (IFN)-response and NFkB signaling modules. B) Eigengene-disease associations are shown for each of 4 identified neural-immune module pairs. The astrocyte and IFN-response modules are upregulated in ASD and SCZ. NFkB signaling is elevated across all three disorders. The microglial module is upregulated in ASD and downregulated in SCZ and BD. C) Top hub genes for each module are shown, along with edges supported by co-expression (light grey; Pearson correlation>0.5) and known protein-protein interactions (dark lines). Nodes follow same coloring scheme as in Fig 5C. Hubs in the astrocyte module (geneM3/isoM1) include several canonical, specific astrocyte markers, including SOX9, GJA1, SPON1, and NOTCH2. Microglial module hub genes include canonical, specific microglial markers, including AIF1, CSF1R, TYROBP, TMEM119. The NFkB module includes many known downstream transcription factor targets (JAK3, STAT3, JUNB, FOS) and upstream activators (IL1R1, 9 TNF receptor superfamily members) of this pathway. D) The top 4 GO enrichments are shown for each module. E) Module enrichment for known cell type-specific marker genes, collated from sequencing studies of neural-immune cell types (–95). F) Module eigengene expression across age demonstrates distinct and dynamic neural-immune trajectories for each disorder.

Figure 8.. LncRNA annotation, ANK2 isoform switching & microexon enrichment

A) FISH images demonstrate interneuron expression for two poorly annotated lincRNAs – LINC00643 and LINC01166 – in area 9 of adult human prefrontal cortex. Sections were labeled with GAD1 probe (green) to indicate GABAergic neurons and lncRNA (magenta) probes for LINC00643 (left) or for LINC01166 (right). All sections were counterstained with DAPI (blue) to reveal cell nuclei. Lipofuscin autofluorescence is visible in both the green and red channels and appears yellow/orange. Scale bar, 10 μm. FISH was repeated at least twice on independent samples (Table S9 (21)) with similar results (see also Fig S16). B) ANK2 isoforms ANK2–006 and ANK2–013 show significant DTU in SCZ and ASD, respectively (*FDR<0.05). C) Exon structure of ANK2 highlighting (dashed lines) the ANK2–006 and ANK2–013 isoforms. Inset, these isoforms have different protein domains and carry different microexons. ANK2–006 is hit by multiple ASD DNMs while ANK2–013 could be entirely eliminated by a de novo CNV deletion in ASD. D) Disease-specific coexpressed PPI network. Both ANK2–006 and ANK2–013 interact with NRCAM. The ASD-associated isoform ANK2–013 has two additional interacting partners, SCN4B and TAF9. E) As a class, switch isoforms are significantly enriched in microexon(s). In contrast, exons of average length are not enriched among switch isoforms. Y-axis displays odds ratio on log₂ scale. P-values are calculated using logistic regression and corrected for multiple comparisons. F) Enrichment of 64 genes with switch isoforms in: ASD risk loci (81); CHD8 targets (96); FMRP targets (33); Mutationally constraint genes (97); Syndromic and highly ranked (1 and 2) genes from SFARI Gene database; Vulnerable ASD genes (98); Genes with probability of loss-of-function intolerance (pLI) > 0.99 as reported by the Exome Aggregation Consortium (99); Genes with likely-gene-disruption (LGD) or LGD plus missense de novo mutations (DNMs) found in patients with neurodevelopmental disorders (21).