
==== Front
Biol Psychiatry Glob Open Sci
Biol Psychiatry Glob Open Sci
Biological Psychiatry Global Open Science
2667-1743
Elsevier

S2667-1743(24)00078-8
10.1016/j.bpsgos.2024.100365
100365
Archival Report
Integrated Long Noncoding RNA and Messenger RNA Expression Analysis Identifies Molecules Specifically Associated With Resiliency and Susceptibility to Depression and Antidepressant Response
Wang Qingzhong wangqingzhong3@gmail.com
a∗
Wang Huizhen a
Dwivedi Yogesh ydwivedi@uab.edu
b∗
a Institute of Chinese Materia Medica, Shanghai University of Traditional Chinese Medicine, Shanghai, China
b Department of Psychiatry and Behavioral Neurobiology, University of Alabama at Birmingham, Birmingham, Alabama
∗ Address correspondence to Yogesh Dwivedi, Ph.D. ydwivedi@uab.edu
∗ Qingzhong Wang, Ph.D. wangqingzhong3@gmail.com
20 7 2024
11 2024
20 7 2024
4 6 10036523 4 2024
18 6 2024
2 7 2024
© 2024 The Authors
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Background

Depression involves maladaptive processes impairing an individual’s ability to interface with the environment appropriately. Long noncoding RNAs (lncRNAs) are gaining traction for their role in higher-order brain functioning. Recently, we reported that lncRNA coexpression modules may underlie abnormal responses to stress in rats showing depression-like behavior. The current study explored the global expression regulation of lncRNAs and messenger RNAs (mRNAs) in the hippocampus of rats showing susceptibility (learned helplessness [LH]) or resiliency (non-LH) to depression and fluoxetine response to LH (LH+FLX).

Methods

Multiple comparison analysis was performed with an analysis of variance via the aov and summary function in the R platform to identify the differential expression of mRNAs and lncRNAs among LH, non-LH, tested control, and LH+FLX groups. Weighted gene coexpression network analysis was used to identify distinctive modules and pathways associated with each phenotype. A machine learning analysis was conducted to screen the critical target genes. Based on the combined analysis, the regulatory effects of lncRNAs on mRNA expression were explored.

Results

Multiple comparison analyses revealed differentially expressed mRNAs and lncRNAs with each phenotype. Integrated bioinformatics analysis identified novel transcripts, specific modules, and regulatory pairs of mRNA-lncRNA in each phenotype. In addition, the machine learning approach predicted lncRNA-regulated Spp2 and Olr25 genes in developing LH behavior, whereas joint analysis of mRNA-lncRNA pairs identified Mboat7, Lmod1, Il18, and Rfx5 genes in depression-like behavior and Adam6 and Tpra1 in antidepressant response.

Conclusions

The study shows a novel role for lncRNAs in the development of specific depression phenotypes and in identifying newer targets for therapeutic development.

Plain Language Summary

We explored transcriptional signatures and regulatory patterns of mRNA and lncRNA in the rat hippocampus of a learned helplessness animal model, including stress-induced depression susceptibility and resilience and changes after antidepressant fluoxetine treatment to learned helpless rats. With the help of integrated bioinformatics analysis, we identified novel transcripts, specific modules, and mRNA-lncRNA regulatory pairs in each phenotype. This study built the foundation for the identification of specific drug targets for depression susceptibility and resilience.

Plain Language Summary

We explored transcriptional signatures and regulatory patterns of mRNA and lncRNA in the rat hippocampus of a learned helplessness animal model, including stress-induced depression susceptibility and resilience and changes after antidepressant fluoxetine treatment to learned helpless rats. With the help of integrated bioinformatics analysis, we identified novel transcripts, specific modules, and mRNA-lncRNA regulatory pairs in each phenotype. This study built the foundation for the identification of specific drug targets for depression susceptibility and resilience.

Keywords

Antidepressants
Depression
lncRNAs
Machine learning
Resilience
==== Body
pmcMajor depressive disorder (MDD) is a common mental health condition that affects approximately 2% to 5% of the global population and is associated with significantly higher levels of morbidity, disability, and mortality (1, 2, 3, 4, 5, 6). Although treatments such as antidepressants and psychotherapy are available, less than half of patients with MDD benefit from them, and a significant number of patients do not respond at all (7, 8, 9). One of the reasons for this could be the unique response to stress in each individual (10). These differences can also be seen in how depression behavior develops (11). Some individuals are more susceptible to depressive symptoms when exposed to early-life trauma or adverse environmental factors while others exhibit resilience (12). Similar patterns have been recapitulated in rodents, where greater vulnerability or resistance to develop depression-like behavior has been shown under chronic stress (13). The learned helplessness (LH) rodent model of depression-like behavior can successfully categorize animals into susceptible and resilient groups when subjected to random, inescapable shocks, which could aid in exploring the molecular and neurobiological mechanisms responsible for individual differences in stress responsiveness during the development of depressive behavior (14,15).

Previous studies have shown that the neurobiology of susceptibility and resilience to depression could involve alterations in neuronal circuits, neuroinflammation, and exacerbation of the endocrine system (16). For example, stressful conditions can cause immune activation and elevation of cytokines that can lead to metabolic disorders associated with the tryptophan kynurenine pathway. Dysregulation of the kynurenine pathway may result in an imbalance of neuroactive metabolites, which has been suggested to serve as a biomarker for depression and suicidal behavior in patients with MDD (17). More recently, increasing attention has been paid to identifying the transcriptomic and epigenome regulators that can affect depression susceptibility and resiliency (18). In this respect, noncoding RNAs are gaining momentum for their role in disease pathogenesis, including major depression (19). Interestingly, approximately 98% of all transcriptional output in humans is noncoding RNA, whereas only 2% belong to protein-coding genes. Among various noncoding RNAs, long noncoding RNAs (lncRNAs) are defined as those that are >200 nt long (19). Although lncRNAs are expressed throughout the body, approximately 40% are expressed specifically in the brain, suggesting brain-specific roles for lncRNAs (20). Based on genomic organization, lncRNAs can be categorized as intragenic, intergenic, and enhancer (21). Depending on their subcellular localization in the nucleus or cytoplasm, lncRNAs are involved in key transcriptomic and epigenetic processes, including alternative splicing, RNA subcellular localization, RNA stabilization, competing endogenous RNA function, chromosome silencing, chromatin modification, and interference (22). In addition, lncRNAs mediate transcription by regulating the interaction of transcription factors and enhancers by binding enhancers to regulate their activity (23). Some lncRNAs function as microRNA sponges where they bind microRNA response elements to alleviate target messenger RNA (mRNA) suppression mediated by microRNAs (24).

While there have been some reports of associations between lncRNAs and psychiatric phenotypes, little is known about their role in depression (25). A few studies have found lncRNAs as useful biomarkers and therapeutic targets for depression (26). Increased expression of lncRNA NONHSAG045500, whose function is to reduce the expression of the 5-HT transporter, thereby inhibiting the transmission of 5-HT, has been reported in patients with MDD (27). Another study showed that lncRNA TCONS 00019174 may serve as a potential diagnostic and therapeutic biomarker for MDD (26). Viral-mediated lncRNA TCONS 00019174 overexpression in hippocampal neurons improved depression-like behaviors in mice exposed to chronic ultramild stress. We recently showed that an overrepresented class of lncRNAs in the hippocampus was associated with resiliency (28). In a subsequent study, we reported that lncRNA coexpression modules may underlie normal and aberrant responses to stress, providing evidence that lncRNA-associated complex trait-specific networks may play a crucial role in developing depression-like behavior (29). In the present study, we hypothesized that there is a specific regulatory relationship between mRNA and lncRNA associated with stress resiliency, which would mediate certain biological functions, giving rise to specific phenotypes.

Using mRNA and lncRNA transcription profiling datasets from our previous studies (28,29), we examined 1) global perspective of depression phenotype based on multiple comparisons among rats showing depression susceptibility (termed as LH), resiliency (termed as non-LH [NLH]), and antidepressant response (LH rats treated with fluoxetine [LH+FLX]); 2) depression resiliency–, susceptibility-, and FLX treatment–specific lncRNA coexpression patterns; and 3) key genes and regulatory pairs of lncRNAs and mRNAs associated with different phenotypes using a machine learning approach.

Methods and Materials

Animals

Male Sprague Dawley rats (Holtzman strain; ages 6–8 weeks) were provided ad libitum food and water and were acclimatized for 1 week before the experiment. The protocol to induce LH behavior was approved by the Institutional Animal Care and Use Committee of the University of Alabama at Birmingham.

Induction of LH Behavior

The detailed protocol for the induction of LH behavior has been described in our earlier publication (28). A total of 100 inescapable tail shocks were given to randomly selected rats with an intensity of 1 mA for 5 seconds, at an average interval of 60 seconds. Twenty-four hours later, the escape latency test was performed using 2 independent trials. Based on escape latency in the fixed ratio trial, rats were divided into 2 groups: LH (showing escape latency ≥20 seconds) and NLH (showing escape latency <20 seconds). Twenty-four hours after the final escape latency test, rats were decapitated, and brains were dissected out. Hippocampi were isolated from 6 tested controls (TCs), 7 NLH rats, and 7 LH rats and flash frozen in liquid nitrogen. Tissues were stored at −80 °C until they were analyzed.

FLX Treatment to LH Rats

As shown in Figure S1 in Supplement 1 and reported earlier (28), intraperitoneal injections of FLX were given to 7 randomly selected LH rats (LH+FLX) with a dose of 5 mg/kg once daily for 13 days. All the rats were decapitated 24 hours after the final escape latency test, and their hippocampi were dissected. Detailed escape latencies in each group are provided in Supplement 1 (28).

Datasets and Phenotypic Groups in the LH Rat Model

In the present study, hippocampal mRNA and lncRNA transcription profiling datasets were sourced from the previous literature generated in our laboratory (28,29). The phenotypic information includes the following 4 groups: depression-resilient (NLH), depression-susceptible (LH), TC, and LH rats treated with FLX (LH+FLX). Unlike our previously published analysis, in this study we conducted multiple comparisons rather than a groupwise analysis to identify the differential expression of mRNAs and lncRNAs among the 4 groups (29). In addition, we focused on the targets and modules involved in developing specific phenotypes. A schematic diagram shows the experimental procedures and analysis plan (Figure S1 in Supplement 1).

Multiple Comparison Analysis

To identify the most significant changes among 4 phenotypes from lncRNA and mRNA datasets, we conducted a multiple comparison statistical analysis with an analysis of covariance via the aov and summary function in the R platform for the mRNA and lncRNA expression profiling. The function of aov was used to determine whether or not there is a statistically significant difference between the means of 3 or more independent groups. The summary function was used to summarize the analysis of the variance model, including the p value associated with the F statistic, which provided significant differences in the mean values between the 4 groups. To perform analysis of variance for the genome-wide transcripts, we also combined for-loop programming language in R. The for-loop is one of the primary control flow structures that can be used to iterate over a collection of objects, such as a vector, list, matrix, or data frame, which apply the same set of operations to each item of a given data structure (https://www.statology.org/interpret-anova-results-in-r/). Similarly, multiple comparisons statistical analysis was performed on the lncRNA expression profiling among the groups. A p value < .05 was defined as significant.

Weighted Gene Coexpression Network Analysis

To identify distinctive gene coexpression patterns (modules) and pathways associated with the LH, NLH, and LH+FLX phenotypes, we performed weighted gene coexpression network analysis (WGCNA) to obtain the significant module and hub genes. Unlike the previous study (29), we focused on the coexpression pattern of the differential expression mRNAs and lncRNAs calculated from multiple comparisons instead of groupwise analysis and conducted an association between the module and multiple phenotype analysis to explore the susceptibility-, resiliency-, and FLX-specific associated modules. We extracted the differentially expressed gene (DEG) expression values as the expression analysis matrix and conducted WGCNA. The function pickSoftThreshold was used to determine the right number as a soft threshold which can power the correlation of the genes. Then, after calculating the adjacencies, the adjacency matrix was transformed into the topological overlap matrix using the function TOMsimilarity. Cluster analysis was produced based on the hclust function and dynamicTreeCut, which helped classify a group of genes with a high topological overlap index into 1 module labeled with 1 color name. Considering the number of DEGs and the default value, the MinModuleSize of each module was set as 30 genes. Meanwhile, the module labeled gray as a color or as 0 in numeric corresponds to the set of genes that have not been clustered in any module. Thus, the genes in the gray modules were not involved in subsequent research. Subsequently, the Spearman correlation coefficient from each module’s Moduleeigengene and phenotype trait association analysis was used to analyze the correlation between the module and phenotypic characteristics. Then, we selected the key modules with the most significant correlation coefficients and the smallest p values for further study. Finally, the chooseTopHubInEachModule function was used to identify the hub genes in the associated modules. The clusterProfile package was used to enrich the genes in the modules into signal pathways.

Machine Learning Analysis

To screen critical target genes of each phenotype, mRNA and lncRNA datasets were analyzed from different perspectives. First, the limma package was used to perform differential expression analyses of 3 comparisons, including TC vs. LH, TC vs. NLH, and TC vs. LH+FLX. The overlapping DEGs between each comparison were plotted as a Venn diagram. Differential gene expression changes specific to each type were sourced from the nonoverlapping section. The specifically enriched genes via machine learning analysis were considered as candidate novel targets. The machine learning model was built with a random forest algorithm and cross-validation analysis. We first optimized the parameters of mtry (number of randomly drawn candidate variables at each split) and ntree (number of trees) with the loop algorithm to obtain the lowest estimate error rate to measure the discrimination capability between classifiers. The Gini index (GI) was applied to evaluate the importance of different transcripts and lncRNAs, which means that the higher the values of the GI, the more critical transcripts and lncRNAs were. The selected candidate genes with higher GI were validated with cross-validation methods to ensure the accurate prediction for the altered genes. To avoid an overfitting analysis, we used the cross-validation approach in this study. In the cross-validation method, we split the dataset into multiple subsets built and trained on different combinations. We evaluated the model performance in both trained and untrained subsets. This cross-validation can identify and avoid potential overfitting.

Combined Analysis of Key Genes and lncRNAs

We next investigated the regulatory effects of lncRNAs on mRNA expression. The chromosomal positions of all the mRNAs and lncRNAs were retrieved from the UCSC Known Genes, Ensembl, and RefSeq databases. Considering that the lncRNAs have potential cis-regulatory expression, we extended the chromosome region of mRNA into the upstream 10 kb and downstream 10 kb. By means of BedTools software, we retrieved the lncRNA-mRNA pairs that overlapped the chromosome regions between the lncRNAs and extended mRNAs. After the chromosomal regions of lncRNAs were mapped into the corresponding transcripts, correlation coefficients between lncRNA and expression levels of transcripts mRNA in each sample were calculated using the COR function. By comparing the mean correlation coefficients between the case and control groups, we focused on lncRNA-mRNA pairs with absolute correlation coefficients >0.8, which were defined as highly correlated lncRNAs with gene expression levels. We also investigated the lncRNA-mRNA pairs with opposite correlation coefficients between the 2 groups.

Results

Identification of DEGs and lncRNAs in Developing Phenotypes Using Multiple Comparison Analysis

A total of 2842 differentially expressed lncRNAs and 3782 DEGs were screened out from multiple comparison analyses. The biological information of lncRNAs and mRNAs is presented in Tables S1 and S2 in Supplement 2, respectively. Heatmaps showing expression values of the top 100 significantly altered lncRNAs and mRNAs are shown in Figure 1A and B, respectively. Uc.354+ (p = .00016, 100% conserved in human OBI1-AS), MRAK035806 (p = .00016, 89.4% conserved in human RPRD1B), and X64411(p = .00032, 90% conserved in human UBR5) were ranked as the top significantly altered lncRNAs in the global 4-group analysis. Among the top significantly altered genes, olfactory receptor (OR) genes such as Olr594 (p = 5.48 × 10−5), Olr1397 (p = .000073), Olr601 (p = .00011), Olr1117 (p = .0002201), and Olr901 (p = .00031) were present, suggesting their prominent role in the development of different phenotypes.Figure 1 Heatmaps of the top 100 dysregulated lncRNAs (A) and mRNAs (B) that are differentially expressed across TC, LH, NLH, and LH+FLX phenotypes. Data are from 6 TCs, 7 NLH, and 7 LH rats. FLX, fluoxetine; LH, learned helplessness; lncRNA, long noncoding RNA; mRNA, messenger RNA; NLH, nonlearned helplessness; TC, tested control.

Screening of Key lncRNAs and mRNAs as Candidate Targets for Each Phenotype

To identify the key mRNAs and lncRNAs that can influence the development of depression resiliency and susceptibility and FLX treatment response, we performed the comparison and contrast analysis among TC, LH, NLH, and LH+FLX groups. Venn diagrams showing the distribution of lncRNAs and mRNAs are shown in Figure 2A and B, respectively. We found that 658 lncRNAs were associated with depression susceptibility, 285 with resiliency, and 2678 with FLX response in the susceptible group. Similarly, 654 mRNAs were associated with depression susceptibility, 368 with resiliency, and 3663 with FLX response in the susceptible group.Figure 2 Venn diagram showing overlapping and distinct differentially expressed long noncoding RNAs (A) and messenger RNAs (B) from various comparisons. The LH group was compared with the tested control group, the NLH group was compared with the tested control group, and the LH+FLX group was compared with the tested control group. FLX, fluoxetine; LH, learned helplessness; NLH, nonlearned helplessness.

Machine Learning Approach to Identify lncRNAs and mRNAs in Predicting Specific Phenotypic Development

The expression levels of the abovementioned lncRNAs (Tables S3–S5 in Supplement 2) and mRNAs (Tables S6–S8 in Supplement 2) were used to establish the machine learning model. We primarily focused on the GI of each gene and the area under the curve (AUC) for machine learning models. Mean decrease in GI of MRAK080245 (GI = 0.440, 90% conserved human GPR161) and MRAK052585 (GI = 0.273, 87% conserved human FBXL6) were the top-ranked lncRNAs in predicting depressive-like behavior (AUC = 0.6875) (Figure 3A), whereas Spp2 (GI = 0.109) and Olr25 (GI = 0.085) were the top genes in predicting depressive behavior (AUC = 0.96) (Figure 3B).Figure 3 Mean decrease in Gini index using random forest model for lncRNAs and mRNAs in various comparison groups. (A, B) Learned helplessness vs. tested control group, (C, D) nonlearned helplessness vs. tested control group, and (E, F) learned helplessness + fluoxetine-treated group vs. tested control group. lncRNA, long noncoding RNA; mRNA, messenger RNA.

For the resiliency phenotype, AY301282 (GI = 0.228, 88.1% conserved human CX3CL1), MRAK084646 (GI = 0.199, 91.4% conserved human NAGK), and BC064030 (GI = 0.168, 89.2% conserved human RPS16) were the top significantly predictive risk lncRNAs for the occurrence of resiliency (AUC = 0.75) (Figure 3C). In contrast, Pttg1 (GI = 0.164), Cesl1 (GI = 0.161), and Ccl9 (GI = 0.162) had the highest GI (importance) for mRNAs, which suggested their possible risk prediction for resiliency (AUC = 0.95) (Figure 3D).

For the FLX-responsive group, the lncRNAs MRAK050995 (GI = 0.170, 92.5% conserved human LINC01738) and U81826 (GI = 0.137, 94.4% conserved human DPH1), and the genes Gfi1 (GI = 0.116), and Srsf12 (GI = 0.113) had the highest GI scores (Figure 3E, F).

Identification of DEGs, Coexpression Patterns, and Hub Genes for 4 Phenotypes

We conducted the WGCNA to further investigate the coexpressed lncRNAs and mRNAs as modules significantly associated with specific phenotypes. DEGs and lncRNAs from multiple comparison analysis were assigned in 8 and 10 modules, respectively (Figure 4A, B). We focused on the most positive and most negative modules associated with depression resiliency, susceptibility, and antidepressant treatment phenotypes. The darkturquoise module was the top negative module in the resilient (NLH) group. The multiple significant gene ontology terms were mainly related to the nucleus, heterocycle metabolic process, and receptor signaling activities, for example, the vasopressin receptor (Figure 5A), indicating the possible mechanism by which rats can resist LH. In contrast, the lightyellow module was the top positively correlated module with susceptibility to depression phenotype, which mainly involved membrane protein complex, synaptic transmission, and acetylcholine activities (Figure 5B). For the FLX treatment phenotype (LH+FLX) group, the blue module had a significant positive association focusing on organelle- and lumen-related terms (Figure 5C). In contrast, the green module served as the most negative module for the treatment group, mainly mediating extracellular space and signaling pathways for sexual organic development and differentiation (Figure 5D).Figure 4 Module-trait relationship analysis using weighted gene coexpression network analysis. (A) Module-phenotype relationship for the assigned module and mRNA expression for different phenotypes (TC, LH, NLH, and LH+FLX). (B) Module-phenotype relationship for the assigned module and lncRNA expression for 4 different phenotypes (TC, LH, NLH, and LH+FLX). FLX, fluoxetine; LH, learned helplessness; lncRNA, long noncoding RNA; mRNA, messenger RNA; NLH, nonlearned helplessness; TC, tested control.

Figure 5 GO terms analysis of significant modules. (A) The darkturquoise module was negatively correlated with the resiliency (nonlearned helplessness) group (r = −0.43, p = .02). (B) The lightyellow module was positively correlated with the susceptibility (learned helplessness) group (r = 0.46, p = .02). (C) The blue module was positively correlated with the fluoxetine group (r = 0.79, p = 9 × 10−7). (D) The green module was negatively correlated with the fluoxetine group (r = −0.84, p = 6 × 10−8). GO, gene ontology; mRNA, messenger RNA; rRNA, ribosomal RNA.

Similarly, the association between coexpression patterns of lncRNAs and different phenotypes was conducted (Figure S2 in Supplement 1). Hub genes play a central role in the coexpression networks. We found that Aamp in the darkturquoise module in the resiliency (NLH) and Cspg4b in the lightyellow module of depression-susceptible (LH) group were hub genes having high connectivity with nodes. In contrast, Gzf1 in the blue module and Rad51 in the green module were hub genes significantly associated with the FLX-treated (LH+FLX) group. Through the connectivity level of the coexpression network, we obtained the hub lncRNAs in each module (Figure S3 in Supplement 1). AB075604 in the brown module and AF487544 in the black module were found to be hub lncRNAs in the module significantly associated with the treatment group. In contrast, MAK1580591 in the red module acted as the role of hub lncRNA, which was associated with susceptibility to depression-like behavior.

Regulatory Relationship Analysis Between mRNAs and lncRNAs

According to the chromosomal location of mRNAs and lncRNAs, we identified a total of 861 pairs of mRNA_lncRNA shared with the flanking regions (Table S9 in Supplement 2). Multiple comparison analyses using mRNA-lncRNA pairs and phenotypes identified 861, 769, and 787 pairs from the LH vs. TC, LH vs. TC, and FLX+LH vs. TC groups, respectively (Tables S10–S12 in Supplement 2). We also considered the correlation between mRNA and lncRNA and calculated correlation coefficients in the LH, NLH, and LH+FLX groups, respectively. RNA_lncRNA pairs with large significant differences in correlation coefficient in the LH vs. TC, NLH vs. TC, and FLX+LH vs. TC groups are presented in Table 1.Table 1 mRNA_lncRNA Pairs With Large Significant Differences in Correlation Coefficient in the LH vs. TC, NLH vs. TC, and FLX+LH vs. TC Groups

Comparisons	Chromosome	mRNA–lncRNA Pairs	Gene Symbol	R∗	TC R	R∗–TC R	p Value, mRNA–lncRNA	p Value, Phenotype × mRNA–lncRNA	
LH vs. TC	chr1	NM_001134978–MRAK038998	Mboat7	−0.949	0.752	−1.700	2.26 × 10−18	4.14 × 10−16	
chr13	NM_001107179–MRAK143269	Lmod1	−0.897	0.748	−1.645	2.38 × 10−26	3.33 × 10−30	
chr8	NM_019165–DQ230327	Il18	−0.744	0.727	−1.471	1.18 × 10−19	1.43 × 10−8	
chr2	NM_001107694–MRAK083603	Rfx5	0.850	−0.871	1.720	1.14 × 10−29	1.16 × 10−11	
NLH vs. TC	chr2	NM_001134548–MRAK007090	Rimoc1	0.859	−0.734	1.593	6.78 × 10−19	9.48 × 10−15	
chr10	NM_001038595–MRAK136707	Dnaja3	−0.726	0.754	−1.479	.00000114	8.5 × 10−16	
chr13	NM_001107179–MRAK143269	Lmod1	−0.708	0.748	−1.457	4.57 × 10−17	5.28 × 10−24	
chr19	NM_001109128–MRAK039312	Fbxl8	−0.324	0.926	−1.250	3.20 × 10−23	.000124	
FLX vs. TC	chr6	NM_138906–BC088254	Adam6	−0.572	0.664	−1.236	1.03 × 10−18	3.72 × 10−6	
chr4	NM_053534–MRBC012634	Tpra1	−0.468	0.740	−1.209	1.24 × 10−9	1.07 × 10−7	
R∗ represents the correlation coefficients of the LH, NLH, and FLX groups compared with LH vs. TC, NLH vs. TC, and FLX vs. TC, respectively. TC R represents the correlation coefficient of the TC group in each comparison. R∗–TC R represents the difference between the 2 numbers of R∗ and TC R. p Value (mRNA–lncRNA) shows the statistical significance between mRNA and lncRNA correlation. p value (phenotype × mRNA–lncRNA) indicates the significant level of the interaction between phenotype, mRNA, and lncRNA.

Chr, chromosome; FLX, fluoxetine; LH, learned helplessness; lncRNA, long noncoding RNA; mRNA, messenger RNA; NLH, nonlearned helplessness; TC, tested control.

From Table 1, 4 groups of mRNA-lncRNA pairs were obtained based on LH vs. TC screening. The corresponding genes were Mboat7, Lmod1, Il18, and Rfx5. Among them, Mboat7, Lmod1, and Il18 were in the LH combination lncRNA, and mRNA had a negative regulatory relationship (R < −0.7); in the TC group, there was a positive correlation (R > 0.7). Rfx5 had a positive correlation in the LH (r = 0.85) and a negative correlation in the TC (r = −0.87) group. Three groups of mRNA-lncRNA pairs were also identified in the NLH vs. TC groups, and the corresponding genes were Rimoc1 (NLH: r = −0.85, TC: r = 0.73), Dnaja3 (NLH: r = −0.73, TC: r = 0.75), Lmod1 (NLH: r = −0.71, TC: r = 0.75), and Fbxl8 (NLH: r = 0.32, TC: r = 0.92). Screening identified 2 lncRNA-mRNA pairs in FLX vs. TC and their expression. The genes were Adam6 (FLX r = −0.51; TC r = 0.66) and Tpra1 (FLX r = −0.47; TC r = 0.74). Interestingly, Lmod1-related lncRNA-mRNA pairs were found not only in the LH group but also in the NLH group.

Discussion

In this study, we focused on susceptibility and resiliency to develop depression and FLX response to differentially expressed mRNAs and lncRNAs from multiple comparison analyses using genome-wide transcriptional profiling data. Potential targets were screened with machine learning to predict different phenotypes. From the view of gene coexpression patterns, we identified and annotated the biological functions of significant modules associated with depression susceptibility, resiliency, and FLX response. Considering that lncRNAs are the essential epigenetic modifiers, we conducted a joint analysis of lncRNAs and mRNAs and found several novel genes responsible for susceptibility (Mboat7, Lmod1, Il18, and Rfx5), resiliency (Rimoc1, Dnaja3, Lmod1, and Fbxl8), and FLX-responsive (Adam6 and Tpra1) phenotypes.

At present, studies are focusing on the biological mechanisms associated with depression susceptibility to identify effective targets to develop potential drugs (30). However, there are relatively few studies on depression resiliency, and there is a lack of an effective animal model that can recapitulate both susceptibility and resiliency to develop depression-like behavior in rodents (13,31, 32, 33). This is critical given that certain individuals are prone to develop depressive behaviors under stressful conditions, whereas some individuals show resistance to depression-like behavior under similar conditions (34). The mechanisms associated with such phenomena are not clearly understood (35). We have modified the classic LH rodent model to construct a dichotomous depression animal model, which can effectively simulate the behavior of individuals showing susceptibility or resiliency to develop depression-like behavior.

The study included 4 different phenotypes, namely depression susceptibility (LH), depression resilience (NLH), antidepressant response to depression (LH+FLX), and nonstressed controls that underwent similar behavioral testing but were not given any stress (TC). In our previous studies, we reported important differentially expressed mRNAs and lncRNAs through groupwise analysis (28,29); however, this analysis was intended to examine the overall expression differences of any transcripts among the 4 groups. In the present study, we implemented the overall expression difference screening of transcriptomic profiling data among 4 groups through the analysis of variance difference for thousands of transcripts simultaneously. We also applied these methods to perform a multiomics integrated analysis of the lncRNAs and mRNAs in-between groups. Finally, we screened the highly correlated mRNA-lncRNAs pairs, which showed significant association with depression, resiliency, and antidepressant response. We have used a similar approach previously to mine the transcriptomic data from participants with depression who did or did not die by suicide and participants without depression who died by suicide to distinguish these 3 phenotypes (36).

A few studies have reported the importance of lncRNAs in depression, particularly those associated with suicide (37). Zhou et al. (25) identified 23 significantly dysregulated lncRNAs through differential expression and weighted coexpression network analyses. They found that 6 lncRNAs were significantly associated with antisense or overlapping protein genes, which were related to interferon signaling as a component of the innate immune response (25). In addition, Punzi et al. (38) revealed a correlation between violent suicide and the expression of lncRNA LOC28758 in the dorsolateral prefrontal cortex (PFC) of participants with schizophrenia, which was located in the MARCKS gene and found to be significantly higher in violent suicide than in nonsuicide or nonviolent suicide groups. Recently, Issler et al. (39) reported that lncRNA FEDORA was enriched in oligodendrocytes and neurons and was upregulated in the PFC of women with depression. Viral delivery of FEDORA selectively promoted depressive-like behavior in female mice and mediated cell type–specific modulation of synaptic properties, myelin thickness, and gene expression (39). Another study also found the primate-specific neuronal enrichment lncRNA gene LINC00473 in the PFC of females with depression, but not males. Using virus-mediated gene transfer of LINC00473 in PFC neurons simulated the human sex-specific phenotype in female mice (40). The abovementioned studies confirmed that lncRNAs, especially FEDORA and LINC00473, can shape the sex-specific landscape of the brain and promote differences in depression in a sex-dependent manner.

From the perspective of lncRNA and mRNA regulation, we found multiple combinations of lncRNAs-mRNAs to be significantly associated with resiliency or depression phenotypes. For example, the correlation coefficients of Rfx5-MRAK083603 in the NLH and TC groups were 0.85 and −0.87, respectively, indicating that there is an opposite regulatory relationship in the NLH and TC groups. Earlier, Zimmermann et al. (41) identified a coexpression network of glucocorticoid-responsive genes that were closely related to stress-induced depression in rodents, Rfx5 being one of the prominent genes in the network. Our results suggest that lncRNA regulates the expression of the Rfx5 gene and may be a key factor in developing resiliency. In terms of susceptibility to develop depressive behavior, we found that the functions of modules significantly related to depressive behavior were primarily related to synaptic transmission. Preclinical and clinical evidence from rodent models of depression and patients with depression repeatedly show that stress exposure can lead to a reduction in the number of synapses in the PFC and hippocampus, atrophy, and loss of neurons and nerve cells (42, 43, 44).

Using a machine learning approach, we found that several OR-related genes (Olr25, Olr529, Olr1561, Olr1064) had higher AUCs that can distinguish specific phenotypes. Molecular correlation analysis based on lncRNA-mediated regulation of olfactory genes supports the role of complex epigenetic regulatory networks that may help in understanding the involvement of lncRNAs in sensory-cognitive dysfunction at the systems level (29). It is pertinent to mention that although OR genes were generally expressed in olfactory-related tissues, more and more studies have demonstrated that OR genes are expressed in multiple tissues and have important biological functions in those nonolfactory tissues (45, 46, 47). For example, it was reported that OR genes (Olfr1505 and Olfr287) are expressed in selected tissues and brain areas outside of the olfactory epithelium and downregulated in Parkinson’s disease postmortem brain tissues (48).

From the regulatory relationship between mRNA and lncRNA, we found that 3 regulatory relationship combinations Mboat7–MRAK038998, Lmod1–MRAK143269, and Il18–DQ230327 mediated the occurrence of depression behavior. These combinations were negatively correlated with the depression-susceptible group, whereas it was positively associated with the TC group. Previously, Yamanishi et al. (49) reported that IL18−/− mice exhibited depression-like behavior and identified Fgfr1, Ptpn1, and Ucn3 genes to be possibly involved in developing depression phenotype. Our present study indicates that lncRNA DQ230327 may inhibit the expression of the IL18 gene, thereby causing susceptibility to developing depression. When we conducted machine learning analysis on mRNA-lncRNA combinations in the FLX-treated LH rats, we identified signaling pathways related to the lumen and sexual organic development and differentiation significantly related to this group. We also found that CREB1 was a promising biomarker of FLX response in LH rats. This supports an earlier study that suggests that a significant gene-gene interaction involving CREB1 can affect response to paroxetine in patients with depression where CREB acts as a transcription factor in regulating plasticity-associated genes (50).

The current study has a few limitations. First, the sample size is relatively small, which may lead to possible discrepancies in machine learning analysis. Second, we lacked biological function and rescue phenotype experiments to identify lncRNAs and mRNAs as targets of antidepressants. Third, all the experiments were done in males, so we do not know whether the changes are sex specific. Further studies are needed to clarify whether there are sex differences in lncRNA species and related coexpression networks in identifying specific phenotypes. Finally, we did not validate the expression level of screened lncRNAs and mRNAs in other animal models. The social defeat stress model is another animal model that has been used to study susceptibility and resiliency to develop depression phenotype. The susceptible animals show anhedonia, social avoidance, and metabolic impairment, which can be reversed by chronic antidepressant treatment (51). It will be interesting to show whether similar changes can be replicated in this animal model. Similarly, the study should also be done in the human postmortem brain to establish relevance to human depression.

Our study has clinical relevance. For example, some people who experience stress are susceptible to developing depression while others display resiliency. Therefore, it is important to know the factors associated with these phenotypes, which can help identify individuals who are prone to develop depression phenotype once they are exposed to stress. Our current study identified several transcripts, modules, biological pathways, and regulatory pairs of mRNA-lncRNA transcripts that may be associated with depression resiliency, susceptibility, and treatment response, which can not only provide their role in developing specific phenotypes but also identify novel targets that can be used to develop therapeutic strategies.

Supplementary Material

Key Resources Table

Supplement 1

Supplement 2

Acknowledgments and Disclosures

This work was supported by the 10.13039/100000025 National Institute of Mental Health (R01MH130539 , R01MH124248 , R01MH118884 , R01MH128994 , R01MH107183 [to YD]) and the 10.13039/501100001809 National Natural Science Foundation of China (31871281 [to QW]), 10.13039/100017950 Shanghai Municipal Health Commission (2020XGKY12 [to QW]), and Scientific Research Foundation for Advanced Talents of Shanghai University of Traditional Chinese Medicine (to QW).

QW and HW analyzed the data. YD and QW designed this study and drafted and revised the article.

The authors declare that all data supporting the findings of this study are available within the article and its supplemental information files. The authors declare that all codes supporting this study’s findings are available from the corresponding author upon reasonable request.

The authors report no biomedical financial interests or potential conflicts of interest.

Supplementary material cited in this article is available online at https://doi.org/10.1016/j.bpsgos.2024.100365.
==== Refs
References

1 Gu L. Xie J. Long J. Chen Q. Chen Q. Pan R. Epidemiology of major depressive disorder in mainland China: A systematic review PLoS One 8 2013 e65356
2 Insel T.R. Charney D.S. Research on major depression: Strategies and priorities JAMA 289 2003 3167 3168 12813123
3 Judd L.L. Mood disorders in the general population represent an important and worldwide public health problem [published correction appears in Int Clin Psychopharmacol 1996;11:153] Int Clin Psychopharmacol 10 1995 5 10
4 Gureje O. Kola L. Afolabi E. Epidemiology of major depressive disorder in elderly Nigerians in the Ibadan Study of Ageing: A community-based survey Lancet 370 2007 957 964 17869636
5 Otte C. Gold S.M. Penninx B.W. Pariante C.M. Etkin A. Fava M. Major depressive disorder Nat Rev Dis Primers 2 2016 16065
6 Baxter A.J. Patton G. Scott K.M. Degenhardt L. Whiteford H.A. Global epidemiology of mental disorders: What are we missing? PLoS One 8 2013 e65514
7 Pandarakalam J.P. Challenges of treatment-resistant depression Psychiatr Danub 30 2018 273 284 30267518
8 Connolly K.R. Thase M.E. If at first you don’t succeed: A review of the evidence for antidepressant augmentation, combination and switching strategies Drugs 71 2011 43 64 21175239
9 Al-Harbi K.S. Treatment-resistant depression: Therapeutic trends, challenges, and future directions Patient Prefer Adherence 6 2012 369 388 22654508
10 Koehn R.K. Bayne B.L. Towards a physiological and genetical understanding of the energetics of the stress response Biol J Linn Soc Lond 37 1989 157 171
11 Zhao L. Han G. Zhao Y. Jin Y. Ge T. Yang W. Gender differences in depression: Evidence from genetics Front Genet 11 2020 562316
12 Remes O. Mendes J.F. Templeton P. Biological, psychological, and social determinants of depression: A review of recent literature Brain Sci 11 2021 1633 34942936
13 Krishnan V. Nestler E.J. Animal models of depression: Molecular perspectives Curr Top Behav Neurosci 7 2011 121 147 21225412
14 Stepanichev M. Dygalo N.N. Grigoryan G. Shishkina G.T. Gulyaeva N. Rodent models of depression: Neurotrophic and neuroinflammatory biomarkers BioMed Res Int 2014 2014 932757
15 Landgraf D. Long J. Der-Avakian A. Streets M. Welsh D.K. Dissociation of learned helplessness and fear conditioning in mice: A mouse model of depression PLoS One 10 2015 e0125892
16 Cathomas F. Murrough J.W. Nestler E.J. Han M.H. Russo S.J. Neurobiology of resilience: Interface between mind and body Biol Psychiatry 86 2019 410 420 31178098
17 Serafini G. Adavastro G. Canepa G. Capobianco L. Conigliaro C. Pittaluga F. Abnormalities in kynurenine pathway metabolism in treatment-resistant depression and suicidality: A systematic review CNS Neurol Disord Drug Targets 16 2017 440 453 28412922
18 Jiang S. Postovit L. Cattaneo A. Binder E.B. Aitchison K.J. Epigenetic modifications in stress response genes associated with childhood trauma Front Psychiatry 10 2019 808 31780969
19 Dwivedi Y. Emerging role of microRNAs in major depressive disorder: Diagnosis and therapeutic implications Dialogues Clin Neurosci 16 2014 43 61 24733970
20 Derrien T. Johnson R. Bussotti G. Tanzer A. Djebali S. Tilgner H. The GENCODE v7 catalog of human long noncoding RNAs: Analysis of their gene structure, evolution, and expression Genome Res 22 2012 1775 1789 22955988
21 Hou Y. Zhang R. Sun X. Enhancer LncRNAs influence chromatin interactions in different ways Front Genet 10 2019 936 31681405
22 Mattick J.S. Amaral P.P. Carninci P. Carpenter S. Chang H.Y. Chen L.L. Long non-coding RNAs: Definitions, functions, challenges and recommendations Nat Rev Mol Cell Biol 24 2023 430 447 36596869
23 Statello L. Guo C.J. Chen L.L. Huarte M. Gene regulation by long non-coding RNAs and its biological functions [published correction appears in Nat Rev Mol Cell Biol 2021;22:159] Nat Rev Mol Cell Biol 22 2021 96 118 33353982
24 Ratti M. Lampis A. Ghidini M. Salati M. Mirchev M.B. Valeri N. Hahne J.C. MicroRNAs (miRNAs) and long non-coding RNAs (lncRNAs) as new tools for cancer therapy: First steps from bench to bedside Target Oncol 15 2020 261 278 32451752
25 Zhou Y. Lutz P.E. Wang Y.C. Ragoussis J. Turecki G. Global long non-coding RNA expression in the rostral anterior cingulate cortex of depressed suicides Transl Psychiatry 8 2018 224 30337518
26 Cui X. Sun X. Niu W. Kong L. He M. Zhong A. Long non-coding RNA: Potential diagnostic and therapeutic biomarker for major depressive disorder Med Sci Monit 22 2016 5240 5248 28039689
27 Cui X. Xu Y. Zhu H. Wang L. Zhou J. Long noncoding RNA NONHSAG045500 regulates serotonin transporter to ameliorate depressive-like behavior via the cAMP-PKA-CREB signaling pathway in a model of perinatal depression J Matern Fetal Neonatal Med 36 2023 2183468
28 Roy B. Wang Q. Dwivedi Y. Long noncoding RNA-associated transcriptomic changes in resiliency or susceptibility to depression and response to antidepressant treatment Int J Neuropsychopharmacol 21 2018 461 472 29390069
29 Wang Q. Roy B. Dwivedi Y. Co-expression network modeling identifies key long non-coding RNA and mRNA modules in altering molecular phenotype to develop stress-induced depression in rats Transl Psychiatry 9 2019 125 30944317
30 Brigitta B. Pathophysiology of depression and mechanisms of treatment Dialogues Clin Neurosci 4 2002 7 20 22033824
31 Wang Q. Timberlake M.A. 2nd Prall K. Dwivedi Y. The recent progress in animal models of depression Prog Neuropsychopharmacol Biol Psychiatry 77 2017 99 109 28396255
32 Planchez B. Surget A. Belzung C. Animal models of major depression: Drawbacks and challenges J Neural Transm (Vienna) 126 2019 1383 1408 31584111
33 Mineur Y.S. Obayemi A. Wigestrand M.B. Fote G.M. Calarco C.A. Li A.M. Picciotto M.R. Cholinergic signaling in the hippocampus regulates social stress resilience and anxiety- and depression-like behavior Proc Natl Acad Sci U S A 110 2013 3573 3578 23401542
34 Dai Q. Smith G.D. Resilience to depression: Implication for psychological vaccination Front Psychiatry 14 2023 1071859
35 Botha M.E. Theory development in perspective: The role of conceptual frameworks and models in theory development J Adv Nurs 14 1989 49 55 2926015
36 Peng S. Zhou Y. Xiong L. Wang Q. Identification of novel targets and pathways to distinguish suicide dependent or independent on depression diagnosis Sci Rep 13 2023 2488 36781900
37 Yoshino Y. Dwivedi Y. Non-coding RNAs in psychiatric disorders and suicidal behavior Front Psychiatry 11 2020 543893
38 Punzi G. Ursini G. Shin J.H. Kleinman J.E. Hyde T.M. Weinberger D.R. Increased expression of MARCKS in post-mortem brain of violent suicide completers is related to transcription of a long, noncoding, antisense RNA Mol Psychiatry 19 2014 1057 1059 24821221
39 Issler O. van der Zee Y.Y. Ramakrishnan A. Xia S. Zinsmaier A.K. Tan C. The long noncoding RNA FEDORA is a cell type- and sex-specific regulator of depression Sci Adv 8 2022 eabn9494
40 Issler O. van der Zee Y.Y. Ramakrishnan A. Wang J. Tan C. Loh Y.E. Sex-specific role for the long non-coding RNA LINC00473 in depression Neuron 106 2020 912 926.e5 32304628
41 Zimmermann C.A. Arloth J. Santarelli S. Löschner A. Weber P. Schmidt M.V. Stress dynamically regulates co-expression networks of glucocorticoid receptor-dependent MDD and SCZ risk genes Transl Psychiatry 9 2019 41 30696808
42 Duman R.S. Aghajanian G.K. Sanacora G. Krystal J.H. Synaptic plasticity and depression: New insights from stress and rapid-acting antidepressants Nat Med 22 2016 238 249 26937618
43 Duman R.S. Pathophysiology of depression and innovative treatments: Remodeling glutamatergic synaptic connections Dialogues Clin Neurosci 16 2014 11 27 24733968
44 Pittenger C. Duman R.S. Stress, depression, and neuroplasticity: A convergence of mechanisms Neuropsychopharmacology 33 2008 88 109 17851537
45 Bradley J. Zhang Y. Bakin R. Lester H.A. Ronnett G.V. Zinn K. Functional expression of the heteromeric “olfactory” cyclic nucleotide-gated channel in the hippocampus: A potential effector of synaptic plasticity in brain neurons J Neurosci 17 1997 1993 2005 9045728
46 Jundi D. Coutanceau J.P. Bullier E. Imarraine S. Fajloun Z. Hong E. Expression of olfactory receptor genes in non-olfactory tissues in the developing and adult zebrafish Sci Rep 13 2023 4651 36944644
47 Wu C. Xu M. Dong J. Cui W. Yuan S. The Structure and Function of Olfactory Receptors Trends Pharmacol Sci 45 2024 268 280 2024 38296675
48 Grison A. Zucchelli S. Urzì A. Zamparo I. Lazarevic D. Pascarella G. Mesencephalic dopaminergic neurons express a repertoire of olfactory receptors and respond to odorant-like molecules BMC Genomics 15 2014 729 25164183
49 Yamanishi K. Hashimoto T. Miyauchi M. Mukai K. Ikubo K. Uwa N. Analysis of genes linked to depressive-like behaviors in interleukin-18-deficient mice: Gene expression profiles in the brain Biomed Rep 12 2020 3 10 31839943
50 Murphy G.M. Jr. Sarginson J.E. Ryan H.S. O’Hara R. Schatzberg A.F. Lazzeroni L.C. BDNF and CREB1 genetic variants interact to affect antidepressant treatment outcomes in geriatric depression Pharmacogenet Genomics 23 2013 301 313 23619509
51 Golden S.A. Covington H.E. 3rd Berton O. Russo S.J.J. A standardized protocol for repeated social defeat stress in mice [published correction appears in Nat Protoc 2015;10:643] Nat Protoc 6 2011 1183 1191 21799487
