Abstract Despite studies providing insight into the neurobiology of chronic stress, depression and anxiety, long noncoding RNA (lncRNA)-mediated mechanisms underlying the common and distinct pathophysiology of these stress-induced disorders remain nonconclusive. In a previous study, we used the chronic mild stress paradigm to separate depression-susceptible, anxiety-susceptible and insusceptible rat subpopulations. In the current study, lncRNA and messenger RNA (mRNA) expression was comparatively profiled in the hippocampus of the three stress groups using microarray technology. Groupwise comparisons identified distinct sets of lncRNAs and mRNAs associated with the three different behavioral phenotypes of the stressed rats. To investigate the regulatory roles of the dysregulated lncRNAs upon mRNA expression, correlations between the differential lncRNAs and mRNAs were first analyzed by combined use of weighted gene coexpression network analysis and ceRNA theory-based methods. Subsequent functional analysis of strongly correlated mRNAs indicated that the dysregulated lncRNAs were involved in various biological pathways and processes to specifically induce rat susceptibility or resiliency to depression or anxiety. Further intersectional analysis of phenotype-associated and drug-associated lncRNA-mRNA networks and subnetworks assisted in identifying 16 hub lncRNAs as potential targets of anti-depression/anxiety drugs. Collectively, our study established the molecular basis for understanding the similarities and differences in pathophysiological mechanisms underlying stress-induced depression or anxiety and stress resiliency, revealing several important lncRNAs that represent potentially new therapeutic drug targets for depression and anxiety disorders. Keywords: Depression, Anxiety, Chronic mild stress, lncRNA, Drug target Highlights * • LncRNAs/mRNAs and pathways related to three CMS-induced phenotypes were identified. * • Specific and common molecular basis of three CMS-induced phenotypes was provided. * • Phenotype- and drug-associated lncRNA-mRNA networks were intersectionally applied. * • 16 hub lncRNAs as potential targets of anti-depression/anxiety drugs were revealed. 1. Introduction Depression and anxiety are two common and chronic psychiatric diseases, having a significant impact on the socio-occupational well-being of patients, family, and society ([41]Almeida et al., 2012; [42]Hamilton et al., 2015; [43]Larson et al., 2007). The core symptoms of depression and anxiety disorders are differentially exhibited but frequently coexist clinically ([44]Brodbeck et al., 2011; [45]Melton et al., 2016). Due to their considerable comorbidity and pathophysiological overlap, results in most basic and clinical studies are frequently mixed ([46]Chiba et al., 2012; [47]Yun et al., 2016), potentially leading to our obscure understanding of the factors moderating these two disorders. Common and specific features in the neural system are still largely unknown, even though some recent studies have slowly begun to separately investigate noncomorbid subjects with these two disorders ([48]Chen et al., 2018; [49]Frick, 2017; [50]Lotan et al., 2014; [51]Zhao et al., 2017). Mounting evidence indicates that depression and anxiety disorders share several putative risk factors, such as life event stress and chronic stress ([52]Leuner and Shors, 2013; [53]Mathew et al., 2011; [54]Pittenger and Duman, 2008). In a given society, chronic stressful life events have been proposed as adverse environmental factors underlying the etiologies of depressive and anxiety disorders ([55]Chang and Grace, 2014; [56]Yun et al., 2016). However, many individuals exposed to chronic stressful events do not display depressive or anxious symptoms ([57]Henningsen et al., 2012; [58]Krishnan et al., 2007; [59]Russo et al., 2012; [60]Uchida et al., 2011). To model some of the environmental factors affecting humans, chronic mild stress (CMS) is commonly employed in rodents ([61]Chang and Grace, 2014; [62]Henningsen et al., 2012). To explore gene-environment interactions in these disorders, genetic and epigenetic dysregulations in adult rats exposed to CMS is generally investigated to gain insight into the biological basis of stress-induced behavioral variations ([63]Argentieri et al., 2017; [64]Golden et al., 2013; [65]Kang et al., 2012; [66]Krishnan et al., 2007; [67]Pena and Nestler, 2018). At the epigenetic level, mRNA transcription, mRNA post-transcription, and long noncoding RNAs (lncRNAs) are critical players in regulating coding gene expression ([68]Yang et al., 2016). Although lncRNAs do not code for proteins, they are annotated due to the presence of a cryptic open reading frame. With the advancement of biomedical research on lncRNAs, increasing evidence indicates that lncRNAs participate in various critical events, such as allosterically mediating enzymatic activity, impacting chromosome conformation, and genomic imprinting ([69]Quinn and Chang, 2016). Moreover, mutations and dysregulation of lncRNA expression have been linked with various developmental processes and disease pathogenesis ([70]Huang et al., 2017; [71]Sunwoo et al., 2016; [72]Yang et al., 2016). Interestingly, lncRNAs are highly expressed in the brain and contribute to key neuronal functions, including neurogenesis, brain patterning, synaptic efficiency, and neural plasticity ([73]Yang et al., 2016). In this vein, lncRNAs may play powerful regulatory roles in CMS-associated depression and anxiety pathologies ([74]Huang et al., 2017; [75]Spadaro et al., 2015), potentially representing a new class of therapeutic targets that can be exploited for disease treatment ([76]Qureshi and Mehler, 2013). Therefore, investigating lncRNA-directed regulatory networks should contribute to our understanding of the mechanisms underlying susceptibility and resilience to depression or anxiety and to drug target discovery ([77]Matsui and Corey, 2017). The hippocampus is well-known as a malleable brain region with respect to stress stimulation and is closely related with learning and memory, as well as emotions. Chronic stress has been shown to detrimentally affect hippocampal neurogenesis and neuroplasticity, consequently influencing learning and memory abilities in depression and anxiety ([78]Bannerman et al., 2014; [79]Malykhin and Coupland, 2015). Previously, we performed quantitative proteomics analysis to explore common and distinct hippocampal proteins associated with depressive and anxiety phenotypes in the CMS rat model ([80]Tang et al., 2019). This model distinguishes between depression-susceptible, anxiety-susceptible, and insusceptible subpopulations, representing the three different responses to CMS ([81]Tang et al., 2019). In the present study, hippocampal tissues from the same batch of CMS rats from the previous publication were used to explore dysregulated lncRNAs that may contribute to chronic stress-induced depression or anxiety. Hippocampal lncRNA and mRNA profiles were comparatively assessed in depression-susceptible, anxiety-susceptible and insusceptible phenotypes, providing a unique chance to decipher the molecular profiles associated with susceptibility or resilience to CMS. To predict the potential functions of these dysregulated lncRNAs, strongly correlated differentially expressed mRNAs were subjected to gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis using the OmicsBean tool. By coexpression analysis, phenotype-associated and drug-associated lncRNA-mRNA networks and subnetworks were constructed to provide new insight into disease pathophysiologies and drug targets. Further integrated analysis of the subnetworks revealed that sixteen hub lncRNAs may represent important new therapeutic drug targets for depression and anxiety disorders. 2. Materials and methods 2.1. Animals Animals used for collecting the brain samples represented three different behavioral phenotypes: (1) CMS-induced depression-susceptible (assessed in the sucrose preference test (SPT) and the forced swimming test (FST)); (2) CMS-induced anxiety-susceptible (assessed in the elevated plus maze test (EMT)); and (3) CMS-introduced but stress-insusceptible. Additional non-CMS-introduced group served as a negative control. Finally the four groups were obtained with five rats in each group. For a more detailed description, please refer to our previous work ([82]Tang et al., 2019). The animal protocols of this study were approved by the Ethics Committee of Chongqing Medical University, and animals were cared for in accordance with the National Institutes of Health protocols for the use and care of laboratory animals (2017013). 2.2. Collection of tissues and extraction of RNA After the tests, animals were decapitated, and brains were removed and dissected on ice. Hippocampal tissues were separated, weighed, rapidly frozen in liquid N[2] and stored at a temperature of −80 °C. From each group, the tissues of four rats were randomly selected for the following microarray analysis. Total RNA was extracted from tissues with TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions. Moreover, we measured RNA quality and quantity with a NanoDrop ND-1000. Then, we evaluated RNA integrity using standard electrophoresis. 2.3. Microarray analysis Arraystar profiling (Rockville, USA) was employed for expression profiling of protein-coding mRNAs and lncRNAs. We selected approximately 9000 lncRNAs from the NCBI RefSeq databases. We performed microarray analysis and sample labeling in accordance with the Agilent One-Color Microarray-Based Gene Expression Analysis protocol. Briefly, after removing rRNA, we performed purification of mRNA using mRNA-ONLY™. Afterward, we transcribed and amplified the samples into fluorescent cDNAs (Arraystar Flash RNA Labeling Kit, Arraystar). A random priming method was used to evaluate transcript length without 3’ bias. Using the Qiagen RNeasy Mini Kit, we purified the labeled cDNAs. With the NanoDrop ND-1000 tool, we examined the specific activity and the concentration of labeled cDNAs (pmol Cy3/μg cDNA). By adding 1 μl of 25 × fragmentation buffer, as well as 5 μl of 10 × blocking agents, we fragmented 1 μg of each labeled cDNA. Moreover, heating of the mixture was performed at 60 °C for 30 min. Afterward, dilution was performed by mixing the labeled cDNAs with 25 μl H[2]O, which was then mixed with GE hybridization buffer. Then, we dispensed the hybridization solution (50 μl) into the gasket slide, which was later assembled onto lncRNA expression microarray slides. The slides were subjected to 17-h incubation at 65 °C in a hybridization oven (Agilent Technologies). Using the Agilent DNA Microarray Scanner (G2505C), hybridized arrays were scanned, washed, and fixed. 2.4. Microarray data analysis Extraction of raw data and analysis of the array images was performed using Feature Extraction software (Agilent version 11.0.1.1). Moreover, we performed quantile normalization using the GeneSpring GX v12.1 software package (Agilent Technologies). Differential expression analysis of normalized signal data was performed using the Limma software package. A log10-fold-change of unity and a false discovery rate (FDR) of less than 0.05 was used as a threshold. For differentially expressed genes, we performed a hierarchical cluster analysis to visualize the data (TreeView). 2.5. Correlation analysis of lncRNAs and mRNAs Weighted correlation network analysis (WGCNA) assists in building a gene network based on coexpression relationships that helps to identify gene modules that are tightly coordinated across entire datasets ([83]Alaei et al., 2018; [84]Ren et al., 2015). We computed pairwise correlations between mRNA and lncRNA. The power function was used to analyze the topological overlap between lncRNAs and mRNAs, serving as a measure to evaluate neighborhood sharing or connection strength. Meanwhile, ceRNA analysis was employed to identify mRNAs and lncRNAs that share biologically meaningful correlations ([85]Zhang et al., 2018). Moreover, we searched for possible miRNA response elements in the mRNAs and lncRNAs. To predict an interaction between lncRNA, miRNA, and mRNA, we overlapped miRNA binding sites on both mRNAs and lncRNAs. Interactions of lncRNAs and mRNAs with miRNAs were predicted using miRNA target prediction software based on RNAhybird and miRanda. Finally, data obtained from the two methods were combined, and the overlapping sets revealed the relatively stronger correlations between lncRNAs and mRNAs. 2.6. Analysis of anti-depression/anxiety drug targets To identify anti-depression and anti-anxiety drugs and their target genes, we searched the DrugBank database (Version 5.1.4) ([86]https://www.drugbank.ca/). Information, including drug name, DrugBank ID and target genes, was collected in this study. The identified genes were considered primary targets of the drugs. For further analysis, species type for these genes was converted from homo sapiens to rattus norvegicus. Subsequently, interactions of the primary target genes were determined from six interaction databases: I2D, APID, Mentha, IntAct, BioGRID and MatrixDB. For each gene, interactions from all six databases were then combined to derive a universal gene set. To obtain reliably interacting genes, the network link detection Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) was used to analyze interactions between primary target genes and their interacting genes in a further step. Network analysis was conducted under the guidance of the medium confidence (STRING score = 0.4) and seven linkage criteria. In the generated network, genes directly interacting with primary target genes were considered secondary targets of the drugs. 2.7. Functional and network analysis For bioinformatics analysis, the multiomics tool OmicsBean ([87]http://www.omicsbean.cn) was used to identify enrichment of cellular component (CC), biological process (BP), and molecular function (MF) based on Gene Ontology (GO) categories ([88]Tang et al., 2019; [89]Xie et al., 2018). KEGG analysis was used to identify significantly enriched pathways, where p < 0.01 was considered significant. Furthermore, lncRNAs and their strongly correlated mRNAs were selected to develop gene expression networks. Pearson's correlation coefficients of greater than 0.99 were adopted for network construction. Calculating the degree of connectivity as previously described ([90]Wang et al., 2019b), hub genes in the networks were identified using CytoHubba, a plugin in Cytoscape software (Version 3.7.1). According to node degree, the top 30 hub genes were displayed in the Cytoscape software. 3. Results 3.1. Identification of differentially expressed lncRNAs and mRNAs of the CMS model rats In the present study, samples used were from the same batch of CMS rats in our previous study ([91]Tang et al., 2019). Based on the behavioral assessments of SPT, FST and EMT, a subset of control, depression-susceptible, anxiety-susceptible, and insusceptible mice was obtained. The CMS paradigm may represent an effective tool for analyzing unique and common neural characteristics of noncomorbid depression and anxiety. Next, we conducted the microarray screening to identify the lncRNA and mRNA profiles in hippocampal tissues from control, depression-susceptible, anxiety-susceptible and insusceptible rats, as shown in [92]Fig. 1A. Using volcano plot filtering, differentially expressed lncRNAs and mRNAs with an FDR of less than 0.05 and greater than 1.5-fold change in both pairwise comparisons were identified ([93]Fig. 1B–D; [94]Supplementary Table S1). In the depression-susceptible group, 454 lncRNAs and 420 mRNAs were upregulated, and 94 lncRNAs and 394 mRNAs were downregulated ([95]Fig. 1B). In the anxiety-susceptible group, 1051 lncRNAs and 945 mRNAs were upregulated, and 570 lncRNAs and 1463 mRNAs were downregulated ([96]Fig. 1C). In the insusceptible group, 848 lncRNAs and 991 mRNAs were upregulated, and 473 lncRNAs and 1304 mRNAs were downregulated ([97]Fig. 1D). Further analysis by Venn diagrams indicated that, among the insusceptible and the two susceptible groups, a set of lncRNAs (n = 1014) and a set of mRNAs (n = 1665) were regulated similarly as a result of exposure to CMS or potential intermediate variants ([98]Fig. 2A). Interestingly, we also observed 465 lncRNAs and 576 mRNAs with similar regulation between the two susceptible (depression and anxiety) groups, which may indicate some common abnormal components in these two CMS-induced disorders. Despite the similarity in the expression changes in the three groups, a substantial percentage (47%) of lncRNAs and mRNAs were unique to the three CMS groups, suggesting that these groups have different molecular mechanisms underlying responses to stress. In addition, according to unsupervised hierarchical cluster analysis, gene expression profiles were categorized into three distinct groups representing the three different responses to CMS ([99]Fig. 2B–C). Fig. 1. [100]Fig. 1 [101]Open in a new tab Analysis of differentially expressed long noncoding RNAs (lncRNAs) and messenger RNAs (mRNAs). (A) Overview of the workflow of the intersectional network analysis. (B, C and D) Identification of differentially expressed lncRNAs and mRNAs in the hippocampus of depression-susceptible (B), anxiety-susceptible (C) and insusceptible (D) groups. Volcano plot showing variations in gene expression. The fold change log (base 2) is on the x-axis, and the negative false log discovery rate (FDR) (base 10) is on the y-axis. Red signifies high relative expression, and blue signifies low relative expression. Dep-Sus, depression-susceptible; Anx-Sus, anxiety-susceptible; Insus, insusceptible; Cont, control. (For interpretation of the references to