RNA-Sequencing and Bioinformatics Analysis of Long Noncoding RNAs and mRNAs in the Prefrontal Cortex of Mice Following Repeated Social Defeat Stress

Background Repeated or continuous chronic psychological stress may induce diverse neuropsychiatric disorders; however, the underlying mechanisms remain unclear. In this study, we explored the expression profiles of long noncoding RNAs (lncRNAs) and mRNAs, along with their biological function and regulatory network, in mice after repeated social defeat (RSD) stress to explore their potential involvement in the development of anxiety-like behaviors. Main Methods RNA-sequencing was used to screen all differentially expressed (DE) lncRNAs and mRNAs between the RSD and control groups. Quantitative real-time polymerase chain reaction (qRT-PCR) was used for confirmation of the RNA-sequencing results. The function of DE lncRNAs was predicted by Gene Ontology (GO) enrichment and pathway analyses of target mRNAs. In addition, the functional regulatory network of the target mRNAs was constructed to reveal potential relationships between lncRNAs and their target genes with bioinformatics approaches. Key Findings In mice experiencing RSD, 373 and 454 lncRNAs, along with 1142 and 654, mRNAs were significantly upregulated and downregulated, respectively. The detailed regulatory network included 126 eligible lncRNA-mRNA pairs. Among them, 14 genes such as Arhgef1, Chchd2, Fam107a, Dlg1, Nova2, Dpf1, and Shank3 involved in neurite growth, neural development, and synaptic plasticity were direct targets of the DE lncRNAs. qRT-PCR of four of the DE lncRNAs and mRNAs confirmed the reliability of RNA-sequencing. GO clustering analyses showed that the top enriched biological process, cellular component, and molecular function terms were synaptic transmission, neuron spine, and glutamate receptor binding, respectively. Further, the top three significant enriched pathways were synaptic adhesion-like molecule (SALM) protein interactions at the synapses, trafficking of α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA) receptors, as well as glutamate binding, activation of AMPA receptors, and synaptic plasticity. Significance Hundreds of lncRNAs and mRNAs are dysregulated after RSD, and many of these lncRNAs might participate in the development of anxiety-like behaviors via multiple complex mechanisms such as target regulation. Available informatics evidence highlighted the likely role of synapse dysfunction and abnormal synaptic neurotransmission in these behaviors. Thus, our findings provide potential candidate biomarkers or intervention targets for chronic psychological stress-induced neuropsychiatric disorders.


Introduction
Psychosocial stress is a significant contributor to the development of mental disorders or behavioral abnormalities such as anxiety, depression, schizophrenia, drug abuse, fear memory formation, or posttraumatic stress disorder [1]. We [2] and others [3] previously reported that social defeat stress, as a model of chronic psychological stress, may significantly induce anxiety-like behaviors. It is estimated that up to approximately 30% of the population is likely to experience various degrees of psychosocial stress in their lifetime, inevitably jeopardizing mental health and leading 2 BioMed Research International to substantial social or economic burdens [4]. Therefore, understanding the pathogenesis underlying the transition from chronic stress to a neuropsychiatric disorder is of great clinical significance and could help to identify possible intervention targets.
Long noncoding RNAs (lncRNAs) are a type of nonprotein coding RNA with a length ranging from 200 nucleotides to several hundred kilobases and mainly play a role in the response to harmful external stimuli [5]. In particular, lncRNAs are highly expressed in the central nervous system with large copy numbers and broad diversity. Accumulating evidence suggests that lncRNAs are involved in various levels of regulation of gene expression, including transcription or posttranscription regulation and epigenetic modification [6]. Notably, individual lncRNAs may have various functions, and such complex roles can orchestrate a combination of diverse biological processes.
Related to this high diversity in the central nervous system, animal model and human studies have suggested the participation of specific lncRNAs in multiple neuropsychiatric disorders, including anxiety, depression [7], and schizophrenia [8], via regulating neurodevelopment, neurodegeneration, and neuroimmunity. Certain lncRNAs have even emerged as candidate biomarkers for disease diagnosis, therapy, or prognosis [9]. However, the differential expression sites, binding sites, acting modes, and complex network of lncRNAs, which involve mRNAs, microRNAs, cis regulatory elements, and competing endogenous RNA networks, make it difficult to elucidate the specific functions of most lncRNAs [10].
Therefore, in the present study, we sought to investigate the roles of lncRNAs in the mechanism by which repeated or continuous chronic psychological stress induces diverse neuropsychiatric disorders using a mouse model of repeated social defeat (RSD) stress, based on a resident and intruder paradigm. Despite the involvement of several brain regions after repeated stress exposure, such as the amygdala or hippocampus, we chose to focus on the prefrontal cortex (PFC), which is considered to be the most sensitive brain region to control fear and stress responses [11]. Even slight acute stress can lead to rapid and prominent injury of prefrontal cognitive ability, while long-term chronic stress may even change the architecture of the PFC dendrites and thus alter its functional connectivity to the rest of the brain [12]. After inducing anxiety-like behavior, a marker of successful model establishment, we examined differentially expressed (DE) lncRNAs and mRNAs in the PFC of the mice using an RNA-sequencing approach. We further confirmed the RNA-sequencing results with four of the identified DE lncR-NAs and mRNAs using quantitative reverse transcriptionpolymerase chain reaction (qRT-PCR). Subsequently, we predicted the biological function of the DE lncRNAs target mRNAs with Gene Ontology (GO) and pathway analyses in the RSD mice, and further constructed systematic functional landscapes of the lncRNAs-mRNAs network with specific bioinformatics approaches. These results can provide new insights and potential intervention targets for the early intervention or treatment of social stress-induced neuropsychiatric disorders.

Materials and Methods
. . Animals. Male 6-8-week-old C57/BL6 mice (20-25 g) were used as the residents, and 12-month-old CD-1 mice (30-35 g, retired breeders) were used as the intruders for RSD model establishment; all mice were purchased from Charles River Laboratory (Beijing, China). The C57BL/6 mice were housed in cohorts of three and the CD-1 mice were maintained in individual cages under a 12-h light/dark cycle at 22-24 ∘ C with sufficient food and water. All mice were acclimated to the surroundings for 7 days before model establishment. All procedures were conducted after ethical approval of the Institutional Animal Care and Use Committee of Nanjing Medical University (Approval no. 1706017).
. . RSD Model Establishment. The RSD model was established according to a resident-intruder paradigm as described previously [13]. In detail, an intruder male CD-1 mouse was placed into the home cage of three resident C57BL/6 mice. During the exposure, we recorded the behaviors of the mice to ensure that the intruder CD-1 mouse was acting aggressively, such as exhibiting attack behaviors, within 5-10 minutes after contact, while the resident C57/BL6 mice were presenting typical submissive behaviors, including an upright posture, escape attempts, or crouching. If no attack behaviors were observed by the CD-1 mouse or if the resident C57/BL6 mice attacked the intruder, a different CD-1 mouse was introduced. Each bout of exposure lasted 2 hours and was conducted for six consecutive nights; all anxiety-like behaviors of the C57/BL6 mice were monitored throughout the experiment. As a control, mice in a separate cage and room were left undisturbed under the same condition (3 mice per cage) throughout the experimental period.
. . Elevated Plus Maze (EPM) Test. The RSD-induced behavioral changes in C57/BL6 mice were evaluated with an EPM test, as a well-proven animal model for anxietylike behaviors detection based on the aversion that rodents have for open spaces [14]. In brief, the apparatus contained four arms (70 × 5 cm) in the shape of a cross placed at a height of 75 cm above the floor. At the beginning of each test, the mouse was set into the center of the maze facing an open arm, and allowed to explore the maze ad libitum. The test was performed in a quiet room with fixed brightness and the movement trajectory was recorded and analyzed using SuperMaze animal behavior analysis software (Xinxin Ltd., Shanghai, China). The total time spent in the open and closed arms was recorded for each mouse with a 5-minute duration. Either an increase in time spent in the closed arm or a decrease in time spent in the open arm is considered to indicate anxiety-like behavior.
. . Tissue Collection and RNA Isolation. Only mice that presented with anxiety-like behavior were subjected to RNAsequencing and bioinformatics analysis. For PFC tissue collection, eligible C57/BL6 mice were immediately decapitated after the EPM test. The brains were removed and cut coronally at the caudal border of the olfactory tubercle to remove . . Library Preparation for lncRNA Sequencing. For RNA sample preparations, 3 g RNA served as the input material for each sample. Sequencing libraries were obtained using NEBNext5 Ultra6 Directional RNA Library Prep Kit for Illu-mina5 (NEB, USA) according to the manufacturer's protocol with index codes supplemented to generate sequences for each sample. In brief, ribosomal RNA was depleted with residues randomly cut into short fragments. First-strand cDNA was produced via a random hexamer primer, and the second strand was synthesized using DNA Polymerase I and RNase H. dNTPs with dTTP were substituted by dUTP in the reaction buffer. The left overhangs were changed into blunt ends with exonuclease and polymerase. Following adenylation of DNA 3 -ends and hairpin loop structure ligation, the library fragments were purified to select desired cDNA fragments of 150-200 bp. Next, 3 l USER Enzyme (NEB, USA) was mixed with the resultant cDNA and incubated at 37 ∘ C for 15 min followed by 95 ∘ C for 5 min. Thereafter, PCR was conducted with Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index (X) Primer. Finally, the products were purified with the AMPure XP system and the library quality was assessed (Agilent Bioanalyzer 2100 system). The generated libraries were sequenced by KeyGEN BioTECH (Nanjing, China) on an Illumina Hiseq 2500 platform, and 100-bp paired-end reads were produced. There are three samples obtained from three mice in each group, with data pooled for statistical analysis.
. . Clustering and Sequencing of lncRNA. Clustering of the index-coded samples was performed with a cBot Cluster Generation System according to the manufacturer's protocol. The library preparations were then sequenced with an Illumina Hiseq system and 150-bp paired-end reads were generated. We controlled the error rates of lncRNAs to be lower than 0.02% for each sample. The transcripts with splices for each sample were screened and combined as lncRNAs using Cuffmerge according to the following criteria: FPKM ≥ 0.5 (Cuffquant), exon number ≥ 2, and length > 200 bp. In addition, with Cuffcompare Software, we excluded overlapping and potential coding transcripts located in the exon regions [16]. All sequencing services were provided by KeyGEN BioTECH (Nanjing, China).

. . Identification of DE lncRNAs and mRNAs.
We determined the DE lncRNAs and mRNAs using the edgeR package [17] along with the Cufflinks algorithm [18] among the transcripts obtained from RNA-sequencing. A |log 2 fold change| ≥ 1.0 and p value ≤ 0.05 were used as the threshold values for identification of DE transcripts. Other than comparing candidate sequences with those of known lncR-NAs with Cuffcompare [19], novel lncRNAs were identified according to the following criteria: (1) transcript length ≥ 200 bp and exon number ≥ 2; (2) a putative open reading frame < 300 bp; (3) overlapping prediction derived from the coding potential calculator, protein family database [20], and codingnoncoding index (CNCI) [21] related to possible lncRNA; and (4) not matched with known lncRNAs.
. . qRT-PCR Validation of DE lncRNAs and mRNAs. Four of the DE lncRNAs and mRNAs were randomly chosen for validation of the RNA-sequencing results in RSD versus control mice using qRT-PCR. In brief, total RNA was obtained from liquid nitrogen-frozen PFC tissues. Aliquots of 1 mg RNA were reverse-transcribed and amplified with a realtime PCR machine (Opticon, MA, USA) according to the manufacturer's recommendations. Expression levels of each lncRNA and mRNA were quantitatively calculated as a fold change with the 2 −ΔΔCT method. The Gapdh gene served as an internal control. The primer sequences used in qPCR are shown in Table 1.
. . Functional Prediction of the Target mRNAs of DE lncRNAs in RSD Mice. Since the majority of lncRNAs identified to date have not yet been functionally described, information on their function largely depends on the function prediction of their target mRNAs. lncTar software was used to predict the target mRNAs of the DE lncRNAs. All of the candidate target genes were then screened one by one to explore their possible involvement in psychiatric disorders with Medline database (PUBMED). Further, GO (http://geneontology.org/) was used to interpret the possible functions by forming hierarchical categories based on the biological process, cellular component, and molecular function terms of the DE target genes. GO is a useful bioinformatics tool that unifies the genes and gene product characteristics across all species and serves a good predictor of gene function. In addition, pathway analyses were performed to identify the significant pathways of the DE targets genes through the Reactome (https://www.reactome.org), Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.jp/kegg/ pathway.html), Protein ANalysis THrough Evolutionary Relationships (PANTHER; http://www.pantherdb.org) and BioCyc (http://biocyc.org) databases. Enrichment factor, gene number, and p values were used to evaluate the correlation between target genes and their predicted functions or pathways. A threshold of p < 0.05 was considered to indicate a significant correlation. Other than KEGG, which stores the higher-order functional information for genes, the Reactome, PANTHER, and BioCyc databases were further screened to systematically elucidate the possible pathways involved in the response to RSD.
. . Analysis of lncRNA-mRNA Regulatory Networks. The coexpression regulatory network was constructed to determine the correlated expressed genes and establish lncRNA-mRNA coexpression pairs. Pearson's correlation coefficient was used with a default filtering threshold set at 0.99 [22]; the corrected p value of the hypothesis test was determined with the Holm calculation, with a default filtering threshold set at 0.05. Correlations among similar types of RNAs were not considered. lncRNA target genes were predicted with lncTar software.
. . Statistical Analysis. Data are expressed as mean ± standard error (sem). Unpaired t-test with Welch's correction was performed for comparison between two groups, and p < 0.05 was considered statistically significant.

. . RSD for Six Days Successfully Induced Anxiety-Like
Behaviors in C /BL Mice. As there is individual variation among mice in response to chronic psychological stress, we examined anxiety-like behaviors in each mouse with the EPM test after six nights of exposure to social defeat stress. All six tested C57/BL6 mice demonstrated a significant increase in the time spent in the closed arm and a decrease in the time spent in the open arm compared to the controls ( Figure 1), suggesting anxiety-like behaviors and successful establishment of the RSD model.

. . Hundreds of DE lncRNAs and mRNAs
Were Identified in RSD Mice. RNA-sequencing revealed that 373 lncRNAs were prominently upregulated and 454 lncRNAs were downregulated in the PFC of RSD mice as compared to that of control C57/BL6 mice. In addition, 1142 and 654 mRNAs were significantly upregulated and downregulated, respectively, in the RSD mice. Data of all differential expressed lncRNAs and mRNAs have been uploaded to GEO database and released since March 5, 2019 (Accession number GSE127812). Tables 2  and 3 list the top 20 upregulated and downregulated lncRNAs and mRNAs in RSD mice, and the volcano plot and heat maps for clustering analysis of DE lncRNAs and mRNAs are shown in Figures 2 and 3, respectively. All differentially expressed lncRNAs and mRNAs, along with fold changes and p value were shown in Supplementary Tables 1-4.

. . qRT-PCR Validated the RNA-Sequencing Results.
qRT-PCR showed that NONMMUT033604.2 and NONMMUT068776.2 were upregulated lncRNAs, and NONMMUT064397.2 and NONMMUT032162.2 were downregulated lncRNAs in RSD mice. In addition, ENSMUST00000106332 and ENSMUST00000088086 were upregulated mRNAs, and ENSMUST00000177637 and ENSMUST00000126073 were downregulated mRNAs in RSD mice. These results were consistent with the RNA-sequencing of these four genes (Figure 4),  (Table 2), which are more likely to play important roles in the pathogenesis in RSD mice. Indeed, a majority of the lncRNAs had no target mRNAs. Thereafter, we screened all possible target mRNAs according to the fold change in the expression level of two upregulated and two downregulated lncRNAs ( Table 4). As shown in Table 4, one lncRNA could potentially target three or more mRNAs, associated with three or more regulatory functions. Thus, as expected, this result highlighted the complex regulatory network between lncRNAs and their target genes, as illustrated in Figure 5. Furthermore, a total of 14 promising target genes of the DE lncRNAs were identified to participate in the pathogenesis of certain psychiatric disorders (Table 5).

. . Functional Prediction of Target mRNAs of DE lncR-NAs in RSD Mice.
As shown in Figure 6, the top three significantly changed biological processes after RSD were found to be synaptic transmission, regulation of postsynaptic membrane potential, and establishment of synaptic vesicle localization. The top three enriched cellular components were neuron spine, postsynaptic density, and excitatory synapse, and the top three enriched molecular   Color key

Discussion
Social defeat is a well-validated model to mimic social stress in rodents [23]. In addition, RSD reproduces the pivotal behavioral, physiological, and immunological changes observed in humans that are subjected to chronic psychosocial stress [24]. As expected, six cycles of social defeat stress successfully induced anxiety-like behaviors in mice, as evidenced by a decrease in time spent in the open arm and in increase of time spent in the closed arm in the EPM test. Moreover, we identified hundreds of lncRNAs and mRNAs that are dysregulated in the mouse PFC after six cycles of social defeat stress, providing a first glimpse into the lncRNAs that might participate in the development of anxiety-like behaviors via multiple complex mechanisms via target regulation. The identified DE lncRNAs were scattered widely across all chromosomes other than sex chromosomes, although the majority are not annotated with available databases. In addition, four of the DE lncRNAs and four mRNAs were verified with qRT-PCR to confirm the reliability of the RNA-sequencing technique. Further, many of the candidate target genes of the DE lncRNAs, such as Arhgef , Chchd , Fam a, Dlg , Nova , Dpf , and Shank , are involved in neurite growth, neural development, and synaptic plasticity. Furthermore, informatics analyses highlighted the contribution of synapse dysfunction and abnormal synaptic neurotransmission. Collectively, our findings provide potential candidate biomarkers or intervention targets   for chronic psychological stress-induced neuropsychiatric disorders.
Accumulating evidence indicates that lncRNAs play important roles in diverse neuropsychiatric disorders. For example, the etiological lncRNA changes associated with major depressive disorder include the cognitive disordersrelated lncRNAs RNON, HAR , PACER, and BACE -AS; synaptic plasticity-related lncRNAs BDNF-AS, DISC , and BC ; or lncRNAs related to other psychiatric diseases such as GOMAFU and MALAT- [23]. Among these, the lncRNA GOMAFU, which is related to schizophrenia, was previously shown to be significantly downregulated in mice with anxiety-like behaviors after fear conditioning [24]. The authors also suggested that GOMAFU interacts with the polycomb repressive complex1, BMI1, which regulates expression of the schizophrenia-related gene beta crystallin (Crybbl), thus playing a regulatory role in fear-induced anxiety-like behaviors.
Of note, most of the lncRNAs expressed in the mammalian brain are regionally enriched in the hippocampus and PFC, two brain regions with prominent roles in the pathogenesis of neuropsychiatric disorders [25]. Although we identified hundreds of DE lncRNAs and mRNAs in the PFC of RSD mice in our study, the majority could not be associated with a target gene. In fact, only nearly 70% of the lncRNAs in the PFC have been characterized in the mouse genome, and the functions of the majority of these lncRNAs are not well defined [25]. Thus, the unidentified lncRNAs warrant further exploration and might also have a relevant function. Overall, only 126 target mRNAs were identified for the DE lncRNAs identified in our study. Accumulating evidence suggests that lncRNAs not only function as enhancers or repressors to regulate gene expression but also act as molecular assemblers to participate in epigenetic modification [26]. Consistent with this idea, we found a complex regulatory network for lncRNAs and their target genes. For example, upregulated NONMMUG093964.1 can downregulate Sgpl and Unc b but can upregulate Chst ; simultaneously, Chst was found to be upregulated by another lncRNA, NONMMUG0956668.1.
Notably, when screening the function of the 126 target genes of the DE lncRNAs, 14 promising genes were found to participate in the pathogenesis of certain psychiatric disorders. For example, Fam a (Drr ), the target gene of the lncRNA MERGE.19055.3, is recognized as acting as a unique link between stress and actin dynamics in the brain. FAM107A binds to actin, stabilizes actin filaments, facilitates bundling, and promotes actin-dependent neurite outgrowth, collectively regulating long-term potentiation and cognitive performance [27]. Increased expression levels of Fam a were previously observed in the PFC, hypothalamus, hippocampus, and lateral septum in response to various types of acute social stress [28] or chronic stress such as social defeat or maternal deprivation [29]. Moreover, dysregulated FAM A expression has been reported in postmortem PFC samples of suicide victims with schizophrenia and bipolar disorder [30]. In addition, Shank , regulated by both MERGE.36317.1 and NONMMUG028155.2, encodes a widely expressed scaffolding protein (SHANK3) in the excitatory synapses, and its mutation leads to significant impairment of hyperpolarization activated cation channels, neuronal morphology, and synaptic connectivity [31]. Shank gene disruption could reproduce autism-like behaviors, including anxiety, social interaction deficits, and repetitive grooming behaviors, in mice [32], while Shank overexpression could induce manic-like behaviors and seizures [33]. SHANK polymorphisms have also been reported in individuals with autism [34], schizophrenia, and bipolar disorder [35]. In contrast to the commonly reported expression regions such as the hypothalamus, ventral tegmental area, and hippocampus, this is the first report demonstrating that the expression of Shank is also dysregulated after social defeat stress in the PFC. Other than these two genes, 18 other genes listed in Table 5 also appear to play potential roles in the development of behavioral abnormalities after RSD. Although their exact contribution is still far from clearly elucidated, our results highlight these genes as promising candidates for subsequent validation studies.
Moreover, GO enrichment analyses revealed that the significantly enriched biological process and cellular component terms related to the target genes of DE lncRNAs in  RSD were mainly concentrated within the realm of synaptic functioning such as synaptic transmission, postsynaptic membrane potential regulation, and postsynaptic density. Further, the significantly enriched molecular function terms were glutamate receptor binding, kinase activity, or transferase activity, which are also involved in regulating synaptic plasticity and transmission. These findings are consistent with previous studies suggesting that synapse dysfunction such as synaptic formation [36], synaptic remodeling [37], and synaptic plasticity [38] plays a pivotal role in the pathogenesis of psychiatric disorders. The significantly enriched pathways largely agreed with the GO results. The top three significantly enriched pathways of the DE genes were SALM protein interactions at the synapses, trafficking of AMPA receptors, as well as glutamate binding, activation of AMPA receptors, and synaptic plasticity. Indeed, SALMs, as a group of cell adhesion molecules, participate in regulation of neurite outgrowth, branching, synapse formation, and synaptic maturation [39]. In particular, SALM member 5 has been reported to be associated with severe progressive autism and familial schizophrenia, suggesting its clinical significance. In addition, dysfunction of AMPA receptors as well as the glutamatergic signaling pathway is well recognized to participate in the pathophysiology of psychiatric and neurological disorders [40]. Stein et al. [41] also suggested that synaptic connections damage due to microglia overactivation and the subsequent overproduction of proinflammatory mediators or cytotoxins significantly contributed to the development of social defeat stress-induced psychiatric disorders such as anxiety or depression. These data collectively indicate that maintenance of synapse function and normal neurotransmission may be an effective strategy for the intervention of psychiatric disorders. Nevertheless, several other GO or pathway terms could also be promising alternative directions for further research, and our findings only provide a preliminary exploration for detailed investigations in this regard.
Despite this potential, there are several limitations of this study that should be acknowledged. First, although use of an animal model is an advantageous tool to obtain valuable bioinformatics-based evidence, the results cannot be directly extrapolated to humans due to species differences. Second, other than the PFC, other brain regions such as the amygdala or hippocampus also play important roles in psychological stress-induced neuropsychiatric disorders, and deserve comprehensive exploration in the same regard. Finally, our results provide numerous candidate biomarkers or intervention targets; however, these candidates are only predicted using available informatics analyses or literature, and thus further research is required for validation of their specific roles.

Conclusion
We provide the first identification of DE lncRNAs and mRNAs in the mouse PFC after RSD stress using RNAsequencing. Furthermore, bioinformatics analysis identified promising candidate target genes and the target regulatory network was constructed. These findings provide valuable insights into the molecular mechanisms underlying psychiatric disorders and also open up new possibilities for the development of novel therapeutic strategies for effective intervention.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.