Integrating Genetic and Gene Co-expression Analysis Identifies Gene Networks Involved in Alcohol and Stress Responses

Abstract

Although the link between stress and alcohol is well recognized, the underlying mechanisms of how they interplay at the molecular level remain unclear. The purpose of this study is to identify molecular networks underlying the effects of alcohol and stress responses, as well as their interaction on anxiety behaviors in the hippocampus of mice using a systems genetics approach. Here, we applied a gene co-expression network approach to transcriptomes of 41 BXD mouse strains under four conditions: stress, alcohol, stress-induced alcohol and control. The co-expression analysis identified 14 modules and characterized four expression patterns across the four conditions. The four expression patterns include up-regulation in no restraint stress and given an ethanol injection (NOE) but restoration in restraint stress followed by an ethanol injection (RSE; pattern 1), down-regulation in NOE but rescue in RSE (pattern 2), up-regulation in both restraint stress followed by a saline injection (RSS) and NOE, and further amplification in RSE (pattern 3), and up-regulation in RSS but reduction in both NOE and RSE (pattern 4). We further identified four functional subnetworks by superimposing protein-protein interactions (PPIs) to the 14 co-expression modules, including γ-aminobutyric acid receptor (GABA) signaling, glutamate signaling, neuropeptide signaling, cAMP-dependent signaling. We further performed module specificity analysis to identify modules that are specific to stress, alcohol, or stress-induced alcohol responses. Finally, we conducted causality analysis to link genetic variation to these identified modules, and anxiety behaviors after stress and alcohol treatments. This study underscores the importance of integrative analysis and offers new insights into the molecular networks underlying stress and alcohol responses.

Keywords: BXD mice strains; alcohol; causality analysis; network analysis; stress; transcriptome.

Figure 1

Diagram showing network-based data analysis. (A) The expression data were generated from the BXD recombinant inbred (RI) mice. (B) The transcripts of these four conditions were compared and analyzed. (C) Weighted gene co-expression network analysis by WGCNA. (D) Conservative analysis of gene modules under experimental conditions. (E) Causal analysis.

Figure 2

Co-expression module analyses. (A) The soft thresholding index R² (y-axis) as a function of different powers β (x-axis). (B) The mean connectivity (y-axis) is a strictly decreasing function of the power β (x-axis). (C) Fourteen co-expression modules identified from the combined dataset. WGCNA cluster dendrogram groups genes (n = 6413) measured across BXD Hippocampus into distinct gene modules (M1–14) defined by dendrogram branch cutting. These modules were significantly enriched for gene ontologies linked to discrete cellular functions and/or organelles in the brain. Genes that did not belong to any modules were housed in the gray modules. The gray gene modules were ignored in this study. (D) Four expression patterns. Four expression patterns were found: up-regulation in No restraint stress and given an ethanol injection (NOE) but restoration in Restraint stress followed by an ethanol injection (RSE) (NOE Up→RSE Restore); down-regulation in NOE but rescue in Restraint stress followed by an ethanol injection (NOE Down→RSE Rescue); up-regulation in both Restraint stress followed by a saline injection (RSS) and NOE, and further amplification in RSE (RSS, NOE Up→RSE Amplify); up-regulation in RSS but reduction in NOE and RSE (RSS Up→NOE, RSE Reduce). One-way analysis of variance (ANOVA) was used to determine the conditions which are significantly different from each other for each expression pattern. ANOVA p-values are indicated in each pattern. The error bar represents standard error of the mean (SEM). (E) Heat maps of Pearson correlation and p-value between modules and traits. Each cell represents the correlation coefficient (and p-value) computing from correlating module eigengenes (MEs) (rows) to traits (columns). Only those correlations with |p| < 0.1 are shown.

Figure 3

Subnetworks were constructed by combining co-expression network and protein-protein interaction (PPI) network. A total of 14 co-expression modules identified by co-expression analysis were superposed with PPI networks from the STRING database, leading to four subnetworks. The network was visualized by Cytoscape software.

Figure 4

Module specificity and preservation. (A) Perseveration of co-expression networks between RSE and RSS. Each square in the graph represents the degree of overlap between two modules. The number in the cell represents the probability of module preservation between two conditions. A hypergeometric two-tailed Fisher’s exact test was used to determine the probability. (B) Summary of module preservations for all four conditions.

Figure 5

Causative analysis. (A) Schematic diagram showing SNP→Module→Phenotype causality analysis. Five possible Single Anchor Models were shown. (B–D) Quantitative trait locus (QTL) mapping for three conditions, including No restraint stress (NOS), RSS and RSE. (E–G) The causal network was constructed by SEM methods for NOE, RSS and RSE.