Integrative gene and isoform co-expression networks reveal regulatory rewiring in stress-related psychiatric disorders

Abstract

Isoform-specific expression patterns have been linked to stress-related psychiatric disorders such as major depressive disorder (MDD). To further explore their involvement, we constructed co-expression networks using total gene expression (TE) and isoform ratio (IR) data from affected (n = 210, 81% with depressive symptoms) and unaffected (n = 95) individuals. Networks were validated using advanced graph generation methods. Our analysis revealed distinct differences in network topology and structure. Shared hubs exhibited unique co-regulatory patterns in each network, with key master hubs in the affected network showing association with psychiatric disorders. Gene Ontology enrichment highlighted condition-specific biological processes linked to each network's master hubs. Notably, isoform-level data uncovered unique co-regulatory interactions and enrichments not observed at the gene level. This is the first study to show network-level differences of gene and isoform co-expression between affected and unaffected individuals of stress-related psychiatric disorders, emphasizing the importance of isoforms in understanding the molecular mechanisms of these conditions.

Keywords: Genetics; Molecular biology; Neuroscience; Omics.

Graphical abstract

Figure 1

Results of differential gene expression (DGE, orange) and differential transcript expression (DTE, gray) analysis (A and B) Volcano plots visualize differentially expressed genes (A) and transcripts (B). See Tables S3, S4, and S7. The dashed line indicates the significance threshold of 5% FDR (n = 450 genes and 269 transcripts). The top 5 up- and down-regulated entities (based on FDR) are labeled. (C and D) Venn diagrams illustrate the overlap of upregulated (C) and downregulated (D) entities, see Figure S3. Orange represents genes, and gray represents genes of transcripts, where transcripts were mapped to their corresponding genes for overlap analysis. (C) Up-regulated: 53 genes found both at gene and transcript level, and 22 genes corresponding to 23 up-regulated transcripts found only at the transcript level. (D) Down-regulated: 101 genes found both at the gene and transcript level, and 81 genes corresponding to 81 down-regulated transcripts found only at the transcript level. (E) Significantly enriched (maximum of 10) Gene Ontology (GO) biological processes (based on BH adjusted p-values) are shown for upregulated genes and transcripts. Rich factor quantifying the degree of enrichment is indicated by the dot size, and the -log10 BH adjusted p-values are represented by the color code, ranging from blue (not significant) to red (significant). Non-significant enrichment values are also indicated by filled circles where available, otherwise left empty (no dots). (F) Enrichment of differentially expressed (DE) genes and transcripts in genes carrying SNPs associated with the GWAS traits schizophrenia, post-traumatic stress disorder, major depressive disorder, bipolar disorder, autistic spectrum disorder, attention deficit hyperactivity disorder, and a psychiatric cross-disorder GWAS. See Table S6. The GWAS trait of height was used as a baseline for comparison. Absolute beta values from the MAGMA analysis are indicated by the dot size, with negative beta values represented by a triangle and positive beta values by a circle. -log 10 of BH adjusted p-values are represented by the color code, ranging from blue (not significant) to red (significant).

Figure 2

Results of network inference and investigation of common hubs (A) Barplot shows the number of total expression (TE, orange) and isoform ratio (IR, blue) nodes, separated by network: affected individuals (AIN, yellow) and unaffected individuals (UIN, green) after filtering and thresholding (AIN: 4115 TE nodes, 3881 IR nodes, 7996 total nodes; UIN: 3514 TE nodes, 4758 IR nodes, 8272 total nodes). (B) PCA plot of the Graph2vec embeddings of the simulated networks generated using the graph generative approach ARROWDiff. The PCA result shows a separate clustering of the 100 simulated UINs (green) and the 100 simulated AINs (yellow), validating the group-specific network topology and the robustness of network inference results. Dots indicate individual simulated networks. (C) Venn diagram shows the overlap of hub genes (degree ≥10) between AIN (Table S8) and UIN (Table S9), with 127 common hubs (Table S10). (D) Venn diagram illustrates the number of distinct connections of the 127 common hubs in the AIN and UIN, and their overlap. (E and F) Visualization of the full affected individuals’ network (AIN, e-bottom left) and unaffected individuals’ network (UIN, f-bottom right). Total gene expression of genes is represented by orange nodes. Isoform ratio (IR) nodes are represented by blue nodes. Zoomed-in views (top panels in E and F) focus on the first-order neighbors of two common hub genes associated with psychiatric disorders: KCD12 and TNFRSF1B. Edges representing connections unique to either the AIN or UIN are represented by solid lines. Transcription and splicing factors are highlighted by red borders around the nodes.

Figure 3

Master hub nodes in the AIN Heatmap of the 61 master hub nodes in the AIN (yellow) exhibiting a greater than 2-fold increase in degree compared to the unthresholded UIN (green). Node degree is indicated by a color gradient from blue (low) to red (high). The 61 master hubs are also annotated for different key information including DE status, roles as TFs or SFs, whether the master hub is also a common hub (degree ≥10 in both networks), and the existence of evidence associating each master hub to psychiatric disorders based either on DisGeNet or on our manual PubMed search (Table S11).

Figure 4

Combined analysis of master hubs in the AIN and the UIN (A) Significantly enriched (maximum of 10) GO biological processes (BH adj p-values) for the 61 master hubs and their first-order neighbors in the AIN and UIN. See Table S12. (B and C) Visualization of the Top 2 master hubs in each network. Total gene expression of genes is represented by orange nodes. Isoform ratio (IR) nodes are represented by blue nodes. Transcription and splicing factors are highlighted by red borders around the nodes. Edges representing connections unique to either the AIN or UIN are represented by solid lines. Network nodes that are differentially expressed in our DE analysis are represented by a triangle shape. (B) Top 2 master hubs (highest degree) in the AIN RC3H1-202 and SRSF6-201 (left) showing a distinct connection pattern compared to their corresponding nodes in the UIN (right). (C) Top 2 master hubs (highest degree) in the UIN MEGF9 and SYNCRIP-202 show distinct connection patterns compared to their corresponding nodes in the AIN (left). (D) Enrichment of the top 2 master hubs and their two-hop neighbors in each network in genes that carry SNPs associated with the GWAS traits schizophrenia, post-traumatic stress disorder, major depressive disorder, bipolar disorder, autistic spectrum disorder, attention deficit hyperactivity disorder, and a psychiatric cross-disorder GWAS. See Table S13. The GWAS trait of height was used as a baseline for comparison. Absolute beta values from the MAGMA analysis are indicated by the dot size, with negative beta values represented by a triangle and positive beta values by a circle. -log10 of BH adjusted p-values are represented by the color code, ranging from blue (not significant) to red (significant).

Figure 5

Isoform-specific network differences (A) Bar plots show the distribution of Isoform Ratio (IR) and Total Expression (TE) nodes in hubs with varying degree thresholds (≥10, ≥15, ≥20) in both the affected individuals’ network (AIN) and the unaffected individuals’ network (UIN). (B and C) HNRNPH1 and HNRNPH1-227 network visualizations: Total gene expression of genes is represented by orange nodes. Isoform ratio (IR) nodes are represented by blue nodes. Transcription and splicing factors are highlighted by red borders around the nodes. Edges represent connections unique to either the AIN or UIN are represented by solid lines. Network nodes that are differentially expressed in our DE analysis are represented by a triangle shape. (B) Network visualization of the heterogeneous nuclear ribonucleoprotein H1 (HNRNPH1) shows a distinct connection pattern in each network. HNRNPH1 is an RNA-binding protein involved in pre-mRNA splicing and RNA processing and is known to be associated with neurodevelopmental disorders. (C) Network visualization of the hub IR node HNRNPH1-227, a transcript of the gene HNRNPH1, and its first-order neighbors. In both networks, first-order neighbors of HNRNPH1-227 are substantially different from the first-order neighbors of its TE node HNRNPH1. (D and E) Significantly enriched (maximum of 10) GO biological processes (BH adj p-values) for master IR nodes and master TE nodes separately and for each network. Non-significant enrichment values are also indicated by filled circles where available, otherwise left empty (no dots). (D) In the AIN, IR master hubs and their first-order neighbors (36 IR master hubs and 627 neighbors) showed significant enrichment in GO biological processes, including in mRNA processing, mRNA metabolic process, positive regulation of catabolic process, and regulatory ncRNA-mediated gene silencing. No significant enrichments were found at the gene-level (TE master hubs and their neighbors). See Table S14. (E) In the UIN, considering IR master hubs and their first-order neighbors (30 IR master hubs and 413 neighbors) showed additional enrichment not found at the gene level, namely in mRNA metabolic processes. See Table S15.