Molecular mechanism of reward treatment ameliorating chronic stress-induced depressive-like behavior assessed by sequencing miRNA and mRNA in medial prefrontal cortex

Chronic stress and lack of reward may reduce the function of the brain's reward circuits, leading to major depressive disorder. The effect of reward treatment on chronic stress-induced depression-like behaviors and its molecular mechanism in the brain remain unclear. In this study, companion communication was used as a reward to study the effect of reward on CUMS-induced depression-like behaviors, and mRNA and miRNA profiles in the medial prefrontal cortex harvested from mice with depression-like and resilient behaviors were established by high-throughput sequencing. The results showed that accompanying with companion ameliorated CUMS-induced depression-like behaviors in mice. Furthermore, 45 differentially expressed genes (DEGs) associated with depression-like behaviors, 8 DEGs associated with resilience and 59 DEGs associated with nature reward (companion) were identified, and 196 differentially expressed miRNAs were found to be associated with companion. Based on the differentially expressed miRNAs and DEGs data, miRNA-mRNA network was established to be associated with companion. Taken together, our data here provided a method to ameliorate depression-like behaviors, and numerous potential drug targets for the prevention or treatment of depression.

hydrate. After the mice are fully anesthetized, perfused with 4 °C normal saline through the left atrium, and decapitated. The mPFC was separated on ice-cold glass slide. Total RNAs from mPFC were isolated in TRIzol Reagent (Life Technology, Carlsbad, CA, USA) as previously described [42], and RNA samples were delivered to Beijing Genomics Institute (BGI) in China for high-throughput sequencing analysis under dry ice conditions. Library Preparation and mRNA-Sequencing. mRNAs were puri ed with Oligo (dT)-attached magnetic beads, and sequentially fragmented into small pieces. Then, rst-strand cDNA was synthesized using random hexamer-primed reverse transcription, followed by a second-strand cDNA synthesis. Obtained cDNA fragments were ampli ed by PCR, and products were puri ed by Ampure XP beads, then dissolved in EB solution. The product was validated on the Agilent Technologies 2100 bioanalyzer for quality control. Then, undergo DSN treatment. The DSN treated library was assessed quality to ensure the high quality of the sequencing data by two methods: checked the distribution of the fragments size using the Agilent 2100 bioanalyzer, and quanti ed the library using real-time quantitative PCR. The quali ed library was ampli ed on cBot to generate the cluster on the ow cell. And the ampli ed ow cell was sequenced single end on the HiSeq4000 platform.
Differentially expressed transcript analysis and data analysis.
The reads per kilo-base per million reads (RPKM) was used to compute the level of transcript expression. Differentially expressed transcripts of samples were determined with DEseq.2. For screening differentially expressed transcripts, the threshold used to identify DEGs were |log2 fold change|≥0.584936 and FDR ≤ 0.05. Subsequently, DEGs were annotated in the gene ontology (GO) database (http://www.geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg) database.
Library Preparation, miRNA-Sequencing and data analysis.
Library was prepared with 1 µg total RNA for each sample. Total RNA was puri ed by electrophoretic separation on a 15% urea denaturing polyacrylamide gel electrophoresis (GAGE) gel and small RNA regions corresponding to the 18-30 nt bands in the marker lane (14-30 ssRNA Ladder Marker, TAKARA) were excised and recovered. Then the18-30 nt small RNAs were ligated to a 5'-adaptor and a 3'-adaptor.
The adapter-ligated small RNAs were subsequently transcribed into cDNA by SuperScript II Reverse Transcriptase (Invitrogen, USA) and then several rounds of PCR ampli cation with PCR Primer Cocktail and PCR Mix were performed to enrich the cDNA fragments. The PCR products were selected by agarose gel electrophoresis with target fragments 100-120 bp, and then puri ed by QIAquick Gel Extraction Kit (QIAGEN, Valencia, CA). The library was quality and quantitated in two methods: check the distribution of the fragments size using the Agilent 2100 bioanalyzer, and quantify the library using real-time quantitative PCR (TaqMan Probe). The nal ligation PCR products were sequenced using the BGISEQ-500 platform.
Subsequently, contaminated reads, including adapter dimers, junk, low complexity, common RNA families (rRNA, tRNA, snRNA, and snoRNA), were disregarded to obtain clean reads. Cleaned tags were annotated with miRBase 21.0 to identify known miRNAs. miRNA expression levels were calculated with TPM values [45]. DESeq software algorithm based on negative binomial distribution and biology duplicate samples was used to compare the known or novel miRNA expression among different groups. The threshold used to identify the different expression of miRNAs was fold-change larger than 1.5 and P-value < 0.05. Three prediction approaches (RNAhybrid, Targetscan, and miRanda) were used to identify the miRNA binding sites.
Integrated miRNA/mRNA network analysis.
A series of bioinformatics analyses found that there is a correlationship between miRNAs and their target mRNAs. miRNAs were usually negatively correlated with their targeted mRNAs in theory, despite a few exceptions [46]. To identify potential miRNA-regulated target genes, the datasets of differentially expressed miRNAs and transcripts were integrated according to the following criteria: (1) In our analysis, miRNAs and mRNAs should undergo reverse changes simultaneously; (2) The correlationship between miRNAs and their target mRNAs should be predicted by the software of RNAhybrid, Targetscan or miRanda. Cytoscape software (San Diego, CA USA) was used to visualize the interactive network of differentially expressed miRNAs and concurrently expressed target mRNAs.
Quantitative RT-PCR for the validations of miRNA and mRNA.
Quantitative real-time RT-PCR (qRT-PCR) were performed with the Bio-Rad CFX 96Touch as previously described [47]. Brie y, total RNA was extracted from mPFC with a TRIzol Kit. cDNA was synthesized using the HiScript RT SuperMix Reagent Kit (Vazyme, R323-01, China) for mRNA, and Mir-X miRNA First-Strand Synthesis Kit (Clontech, 638315, CA, USA) for miRNA. mRNAs were ampli ed in a 20 µl reaction with 1 µl sample cDNA (500 ng/µl), 0.4 µl of each primer (10 nmol/l), 10 µl 2 × ChamQ Universal SYBR qPCR Master Mix and 8.2 µl ddH 2 O. Real-time PCR was initiated at 95℃ for 30 s, followed by 40 cycles of denaturation for 10 s at 95℃, annealing and elongation for 30 s at 60℃, then melt curve 95℃ for 15 s, 60℃ for 1 min, and 95℃ for 15 s. For miRNAs, qRT-PCR was performed in a 20 µl reaction with 1.6 µl sample cDNA, 0.4 µl mRQ3′Primer, 0.4 µl miRNA-speci c Prime (10 µM), 10 µl 2 × SYBR Advantage Premix, and 7.6 µl ddH 2 O. The program was set to 95℃ for 10 s, followed by 40 cycles of denaturation for 5 s at 95℃, annealing and elongation for 20 s at 60℃ and melt curve 65℃ to 95℃ increment 0.5℃ for 5 s. The relative expression level of mRNAs in the tissue was normalized to an internal reference gene GAPDH. The relative expression level of miRNAs in the tissue was normalized to U6 small nucleolar RNA.
Each ampli cation was performed in triplicate. The results were calculated with the 2 −ΔΔCt method.
Dual luciferase reporter assay.
According to the miRNA/mRNA prediction, targeted gene sequence containing targeted sites was ampli ed by PCR. Then, the products were digested by XhoI and NotI and fused into the luciferase vector psiCHEK2 [44]. Site-directed mutation of targeted site was performed with QuikChange Lighting Site-Directed Mutagenesis Kit (Stratagene, La Jolla, CA) according to the manufacturer's instructions.
Luciferase reporter detection assays were carried out as previously described [48]. Brie y, 50 ng psiCHECK2-wild-type or mutant reporter plasmids, 50 nM miRNAs mimic or miR-NC mimic were cotransfected to HEK293T cells by using Lipofectamine 2000 transfection reagent (Invitrogen, Carlsbad, CA, USA). After 24 h, the activities of re y and Renilla luciferase were assessed by using Dual-Glo® Luciferase Assay System (Promega, Cat. E2920, USA) according to the manufacturer's instructions. Each treatment was carried out in triplicates in three independent experiments. Statistical analyses.
The data of behavioral tests, luciferase activity and gene analyses are presented as mean ± SEM. Relationships between miRNA and its target prediction were assessed by Pearson's correlation coe cients. Two-way ANOVA was used to make the statistic comparison among control, CUMS and CUMS-companion groups before and after treatment. Paired t-test was used for statistics of before versus after values within groups. The unpaired Student t-test was used to make the statistic comparison between control and CUMS-MDD, control and CUMS-resilience, control and Reward-MDD, control and Reward-resilience, CUMS-resilience and Reward-resilience, and so on. P < 0.05 is considered statistically signi cant.

Results
Accompanying with companion ameliorated CUMS-induced depression-like behaviors in mice.
In order to identify the effect of rewards on ameliorating depression-like behaviors, as a classic depression model [49], CUMS model was used to explore the possible antidepressant effect of accompanying with companion. As shown in Fig. 1a, C57BL/6J male mice were treated with CUMS or CUMS-companion for 4 weeks, after which depression-like behaviors or resilience were assessed by sucrose preference test, Y-maze test and forced swimming test. According to the behavioral tests, only mice changed signi cantly in all three tests could they be de ned as MDD mice, and mice did not change signi cantly in all three tests could they be de ned as resilience mice (Fig. 1a). The results showed that three behavioral parameters of the CUMS and CUMS-companion groups varied so much after 4 weeks treatments ( Fig. 1b-d). In CUMS group mice, the SPT values signi cantly decreased after 4 weeks treatments (83.86 ± 0.9225% versus 71.51 ± 1.788%, p < 0.0001), n = 67)( Figure 1b and Table S1), the ratios of stay time in M-arm to stay time in total arms also showed signi cant difference after treatments (78.54 ± 1.177% versus 57.22 ± 2.05%, p < 0.0001, n = 67) ( Fig. 1c and Table S2), and the FST's immobile time were 206.6 ± 7.086 seconds after treatments and 137 ± 6.935 seconds before the CUMS treatments (p < 0.0001, n = 67) ( Fig. 1d and Table S3). These results indicated that CUMS treatments can induce depression-like behaviors in mice.
It is noteworthy that pairwise comparisons between the three groups after CUMS showed signi cant differences in SPT, YMT and FST, respectively ( Fig. 1b- Table S1-3). Compared with CUMS group, the SPT and YMT values signi cantly increased in CUMS-companion group, and the ratios of stay time in M-arm to stay time in total arms showed signi cantly decrease in CUMS-companion group ( Fig. 1b- Table  S1-3). In addition, as indicated in Table 1, mice treated by the CUMS or CUMS and companion in 4 weeks, the percentage of MDD in CUMS group was about 29.85%, while that in CUMS-companion group decreased to 12.12%. Meanwhile, the percentage of resilience in CUMS group was about 11.94%, while that in CUMS-companion group increased to 36.36%. These results suggested that accompanying with companion ameliorated CUMS-induced depression-like behaviors in mice. Accompanying with companion disturbed the mRNA expression of medial prefrontal cortex in CUMSinduced MDD and resilience mice.
To investigate the molecular mechanism of ameliorated effect of accompanying with companion in CUMS-induced depression-like behaviors, the miRNA and mRNA pro les were analyzed by highthroughput sequencing in medial prefrontal harvested from control, CUMS-MDD, Reward-MDD, CUMSresilience, Reward-resilience mice. 15 samples were sequenced using RNA-Seq technology, and the number of raw sequencing reads and clean reads are shown in Supplementary Table S4. The unique mapping ratio with reference gene and the average genome mapping ratio were more than 79.48% (Supplementary Table S4). The criterion to make sure differential expression of genes were a 1.5-fold or greater change in transcript level between any two group mice and a P-value < 0.05. As shown in Supplementary Table S5, 311 differentially expressed mRNAs were obtained in control versus CUMS-MDD mice, in which 195 mRNAs were up-regulated and 116 mRNAs were down-regulated in CUMS-MDD mice compared to control mice. In the control versus Reward-MDD mice, 127 differentially expressed mRNAs changed in Reward-MDD mice compared to control mice, where 21 mRNAs were down-regulated and 106 mRNAs were up-regulated. For resilience mice derived from two treatments, 45 differentially expressed mRNAs were obtained in control versus CUMS-resilience mice, in which 24 mRNAs are upregulated and 21 mRNAs are downregulated; 45 differentially expressed mRNAs were obtained in control versus Reward-resilience mice, in which 22 mRNAs are upregulated and 23 mRNAs are downregulated. These results indicate that the molecular changes in the medial prefrontal region are different in different treatments that achieve consistent behaviors.
In order to validate the sequencing data above, we ran quantitative RT-PCR (qRT-PCR) from tissues that were used for mRNA sequencing. As shown in Supplementary Figure S1-S3, the expressions of Ccl28, Ciart, Fkbp5, Fmo2, Gm7120, Gpr149, Hif3a, Kcnh5, Npas4, Plin4, Xdh were decreased, as well as the expressions of Slc6a13 was increased in Reward-MDD mice, compared to control mice ( Figure S1). Moreover, the expressions of Ciart, Dbp, Fmo2, Hif3a, Lrrc39, Xdh were decreased, as well as the expressions of Col3a1, Pcdha1, Slc6a13 and Stac were increased in Reward-resilience mice, compared to control mice ( Figure S2). Consistent results achieved by mRNA sequencing and qRT-PCR con rm the validation of our study.
Identi cation of differentially expressed genes are associated with depression-like behaviors.
Since the depression-like behaviors in MDD mice were consistent whether they were derived from CUMS treated group or CUMS-companion group, the common differentially expressed genes (DEGs) in control versus CUMS-MDD and control versus Reward-MDD comparisons may be strongly associated with depression-like behaviors. To identify both unique and common genes in control versus CUMS-MDD and control versus Reward-MDD comparisons, numbers were calculated and presented using a Venn diagram (Fig. 2). As shown in Fig. 2  To gain in-depth insights into the molecular function of these unique and common DEGs, we performed Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis. As shown in Fig. 2  Identi cation of differentially expressed genes are associated with resilience.
Similar to the depression-like behaviors of MDD mice, the behaviors of resilience mice derived from CUMS treated group and CUMS-companion group were also consistent. Therefore, the common differentially expressed genes (DEGs) are associated with resilience were identi ed, calculated and presented using a Venn diagram (Fig. 3). As shown in Fig. 3 and Supplementary Table S7, 37 DEGs (I) speci cally involved in control versus CUMS-resilience comparison, and 37 unique DEGs (III) in control versus Reward-resilience comparison. In addition, we obtained 8 common DEGs (II) in control versus CUMS-resilience and control versus Reward-resilience comparisons, such as Cpm, Fmo2, Grin1os, Hif3a, Hist2h2aa2, Plin4, Scn5a and Trex1. What is more, the expression pattern of each common gene was consistent in both CUMS-resilience and Reward-resilience mice compared to control mice, either the expression simultaneously increased or decreased. Compared with MDD mice, the resilience mice derived from CUMS treated group and CUMS-companion group seem to share less molecular changes.
Subsequently, Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis was performed for unique DEGs (I), unique DEGs (III) and common DEGs (II). As shown in Fig. 3  Identi cation of differentially expressed genes are associated with accompanying with companion.
Most individuals may not suffer from major depression by chronic stress, i.e., resilience [50]. As shown in Table 1, approximately 11.94% of CUMS-treated mice in 4 weeks were resilience, and this value increased to approximately 36.36% after accompanying with companion. To identify the DEGs are associated with the accompanying with companion treatments, differentially expressed mRNAs were screened in CUMSresilience versus Reward-resilience mice. As shown in Supplementary Table S5, 59 differentially expressed mRNAs were obtained in CUMS-resilience versus Reward-resilience mice, in which 26 mRNAs were up-regulated and 33 mRNAs were down-regulated in CUMS-MDD mice compared to control mice. To gain in-depth insights into the molecular function of these differentially expressed mRNAs, we performed KEGG pathway enrichment analysis. As shown in Table 2, differentially expressed mRNAs were found to be mainly enriched in KEGG pathways associated with Neuroactive ligand-receptor interaction, Dopaminergic synapse, Rap1 signaling pathway, Focal adhesion, Metabolic pathways and PI3K-Akt signaling pathway. Furthermore, to validate the sequencing data above, we ran quantitative RT-PCR (qRT-PCR) from tissues that were used for mRNA sequencing. 10 mRNAs were selected to conduct qRT-PCR, the expressions of Artn, Dio3, Klf4, Pcdha1, Pcdhgb4, and Slc6a13 are increased, as well as the expressions of Col6a3, Drd2, Mmp9, and Rac3 are decreased ( Figure S3). Consistent results achieved by mRNA sequencing and qRT-PCR con rm the validation of our study.
Accompanying with companion altered miRNA expression in resilience mice.
Studies have revealed that patients with psychiatric disorders have altered miRNA expression pro les in the brain and circulation, and multiple roles of miRNAs in psychiatric disorders have summarized [51]. To investigate the impact of accompanying with companion on miRNA, we performed differential expression analysis of miRNA by high-throughput sequencing. As shown in Table S8, compared with the CUMSresilience mice, a total of 61 differentially expressed miRNAs were up-regulated and 135 differentially expressed miRNAs were down-regulated in the Reward-resilience mice. In order to validate the nding by sequencing miRNA analysis, quantitative RT-PCR were performed for certain miRNAs. As shown in Figure   S4, qRT-PCR analysis results for miR-149-5p, miR-183-5p, miR-199b-3p, miR-323-5p, miR-337-3p, miR-384-5p, miR-451a, miR-484, miR-671-5p, and miR-702-3p were consistent with miRNA sequencing results in Reward-resilience versus CUMS-resilience mice.
miRNAs can act as an 'expression switch' and block the expression of their target genes [52]. In this study,if the mRNAs downregulated were caused by miRNAs, their correspondent miRNAs should be upregulated, or vice versa. Based on three databases (RNAhybrid, Targetscan and miRanda), the target genes of miRNAs were predicted and matched with the mRNA measured by mRNA sequencing to form the complex interactions between miRNAs and mRNAs (Tables 3 and 4). As shown in Fig. 4a and 4b, 63 miRNAs established regulatory relationships with 34 mRNAs in CUMS-resilience versus Reward-resilience mice. Furthermore, the miRNAs targeted mRNAs were validated by qRT-PCR and dual luciferase reporter assay. As shown in Fig. 4c and 4d, there were inverse correlations between Pcdha1 and miRNA-337-3p, as well as between Rac3 and miRNA-671-5p from qRT-PCR analysis. In dual luciferase report assay, the relative activities of luciferase reporter for Pcdha1 and Rac3 are signi cantly decreased by the mimics of miRNA-337-3p and miRNA-671-5p respectively, but not negative control ( Fig. 4e and 4f). These reductions can be reversed by changing the binding sites of miRNA-337-3p and miRNA-671-5p in Pcdha1 and Rac3. These results support that mRNA Pcdha1 is the direct target of miRNA-337-3p, and mRNA Rac3 is the direct target of miRNA-671-5p, which is consistent with our bioinformatics analyses in the prediction of miRNA target genes. Note: ↑ indicates up-regulation in the tissue of mPFC from Reward-resilience mice compared to CUMS-resilience mice, whereas ↓ represents down-regulation. Note: ↑ indicates up-regulation in the tissue of mPFC from Reward-resilience mice compared to CUMS-resilience mice, whereas ↓ represents down-regulation.

Gene Symbol
The predicted target miRNAs that match DEGs in transcriptome *

Discussion
Research suggests that major depressive disorder emerges out of a three-pronged series of interacting impairments in the domains of stress regulation, reward and mentalizing [53]. The motivation induced by reward can effectively promote individual's emotional regulation ability [54], and lack of rewards in life can cause the brain's reward circuits to be not su ciently activated or actively used, which may reduce the function of the brain's reward circuits [7,8]. Studies have shown that encouraging patients to participate in reward activities during treatment has been found to be effective in alleviating MDD [27]. In this study, we add a natural reward (companion) on the basic of CUMS model to explore the effect of reward treatment on CUMS-induced depression-like behaviors. Behavioral results showed that accompanying with companion ameliorated depression-like behaviors induced by CUMS treatments.
Subsequently, the high-throughput sequencing for both mRNAs and miRNAs were performed to gain indepth insights into the molecular function. Excitingly, we screened a large number of genes associated with depression-like behaviors, resilience and accompanying with companion.
Depression is a heterogeneous, highly prevalent and moderately inherited disease. In addition, environmental factors, such as sexual, physical or emotional abuse during childhood, are strongly associated with the risk of developing MDD [55]. Therefore, the etiology of MDD is multifactorial, and no established mechanism can explain all aspects of this disease [55]. In this study, two depression-like behavior mice derived from CUMS and reward treatments, although they are phenotypically consistent, there are still large differences in gene expression in the mPFC, which provides an evidence for the different mechanisms of depression in different environments. Nevertheless, we screened 45 coexpressed genes from the mPFC of two depression-like behavior mice, which were de ned as DEGs associated with depression-like behaviors. Subsequent KEGG pathway analysis found that these genes were mainly enriched to Metabolic pathways, cAMP signaling pathway, ECM-receptor interaction, PI3K-Akt signaling pathway, Serotonergic synapse and Neuroactive ligand-receptor interaction. How important these genes are to preventing and treating depression, more research is needed in the future. Generally, most individuals may not suffer from major depression by chronic stress, i.e., resilience [50], implying the presence of endogenous anti-depression in the brain. In this study, we also found resilience mice under the CUMS and CUMS with accompanying with companion, respectively. Compared to 45 DEGs associated with depression-like behaviors, only 8 common DEGs were obtained between control versus CUMS-resilience and control versus Reward-resilience comparisons, which indicated that although resilience derived from CUMS and reward treatments are consistent in behaviors, there are still signi cant differences in molecular changes in the prefrontal brain. Therefore, different antidepressant measures and inherent genetic factors, different antidepressant molecular mechanisms may appear.
Recent studies have shown that low perceived social support and loneliness were related to depression course and recovery [56,57], and only the social support characteristic of few negative experiences with the support from a partner or close friend, and limited feelings of loneliness proved to have unique predictive value for a favorable course of depression [56]. Therefore, it's a great interest to test whether a partner or close friend can contribute to better recovery and prevention of MDD. In our study, mice treated by CUMS and the companion used for accompanying are from the same litter, which is a good way to simulate a partner or close friend. What is more, accompanying with companion achieved a certain effect in ameliorating CUMS-induced depression-like behaviors, and the molecular mechanism were in-depth investigated. Associated with accompanying with companion, 59 mRNAs and 196 miRNAs were identi ed and used to establish the miRNA/mRNA network. Noteworthy, in CUMS-resilience versus Reward-resilience, accompanying with companion disrupted the expression of dopamine receptors (Drd2 and Drd4), which are enrich in dopaminergic synapse and neuroactive ligand-receptor interaction pathways.

Conclusions
Taken together, our present study by analyzing miRNA and mRNA pro les provides numerous potential drug targets for depression and comprehensive view about molecular mechanisms underlying reward treatment ameliorating depressive-like behaviors, which giving better theoretical support for the application of reward intervention in the treatment of depression.  Differently expression miRNAs targeted differently expression mRNAs in CUMS-resilience versus Rewardresilience mice. a-b MicroRNA-mRNA network in CUMS-resilience versus Reward-resilience mice were constructed between the 63 miRNAs and 34 overlapped mRNAs with using transcriptome expression data and predicted target genes from RNAhybrid, Targetscan, and miRanda databases. Red symbols present the elevated expression of miRNAs or mRNAs. Blue symbols present the down regulated miRNAs or