microRNA and mRNA profiles in nucleus accumbens underlying depression versus resilience in response to chronic stress.

Major depression in negative mood is presumably induced by chronic stress with lack of reward. However, most individuals who experience chronic stress demonstrate resilience. Molecular mechanisms underlying stress- induced depression versus resilience remain unknown, which are investigated in brain reward circuits. Mice were treated by chronic unpredictable mild stress (CUMS) for 4 weeks. The tests of sucrose preference, Y-maze, and forced swimming were used to identify depression-like emotion behavior or resilience. High-throughput sequencing was used to analyze mRNA and miRNA quantity in the nucleus accumbens (NAc) harvested from the mice in the groups of control, CUMS-induced depression (CUMS-MDD), and CUMS-resistance to identify molecular profiles of CUMS-MDD versus CUMS-resilience. In data analyses and comparison among three groups, 1.5-fold ratio in reads per kilo-base per million reads (RPKM) was set to judge involvements of mRNA and miRNA in CUMS, MDD, or resilience. The downregulations of serotonergic/dopaminergic synapses, MAPK/calcium signaling pathways, and morphine addiction as well as the upregulations of cAMP/PI3K-Akt signaling pathways and amino acid metabolism are associated with CUMS-MDD. The downregulations of chemokine signaling pathway, synaptic vesicle cycle, and nicotine addiction as well as the upregulations of calcium signaling pathway and tyrosine metabolism are associated with CUMS-resilience. The impairments of serotonergic/dopaminergic synapses and PI3K-Akt/MAPK signaling pathways in the NAc are associated with depression. The upregulation of these entities is associated with resilience. Consistent results from analyzing mRNA/miRNA and using different methods validate our finding and conclusion.

INTRODUCTION

Major depression is presumably caused by chronic stress life plus genetic vulnerability, and is characterized by persistent negative mood, such as anhedonia, interest loss, and low self‐esteem (Camp & Cannon‐Albright, 2005; Hamilton, Chen, & Gotlib, 2013; Jabbi, Korf, Ormel, Kema, & den Boer, 2008; Keers & Uher, 2012; Klengel & Binder, 2013; Lohoff, 2010; Moylan, Maes, Wray, & Berk, 2013; Wilde, Mitchell, Meiser, & Schofield, 2013). Pathological changes in monoamine synapses, hypothalamus‐pituitary‐adrenal axis, and brain‐derived neurotrophic factor are believed to cause neuronal atrophy in the ventral tegmental area, nucleus accumbens, and prefrontal cortex from depressive patients and depression‐like animals (Banasr, Dwyer, & Duman, 2011; Bennett et al., 2008; Duman, 2010; Elizalde et al., 2008; Ma, Xu, et al., 2016; Pittenger & Duman, 2008; Sandi & Haller, 2015; Xu, Cui, & Wang, 2016). In addition, major depression may be caused by a lack of reward in life to lower the use in the function of brain reward circuits including the ventral tegmental area, nucleus accumbens, and prefrontal cortex (Wang & Cui, 2015; Zhu, Wang, Ma, Cui, & Wang, 2017). Conversely, numerous individuals do not suffer from major depression under the condition of chronic stress. The resistance to chronic stress is called as resilience (Southwick & Charney, 2012). The elucidation of endogenous mechanisms underlying resilience to chronic stress should shed light on developing the therapeutic strategies for major depression (Zhu et al., 2017). Although a few of molecules are presumably involved in major depression versus resilience (Bergstrom, Jayatissa, Thykjaer, & Wiborg, 2007; Christensen, Bisgaard, & Wiborg, 2011; Friedman et al., 2014; Manji et al., 2003; Vialou et al., 2010; Wang, Perova, Arenkiel, & Li, 2014). However, comprehensive molecular profiles in specific brain areas remain to be systemically figured in terms of resilience and susceptibility to chronic stress for major depression.

The nucleus accumbens is considered as an important region of brain reward circuits (Basar et al., 2010; Carlezon & Thomas, 2009; Day & Carelli, 2007; Floresco, 2015; Ikemoto, 2007; Salamone, Correa, Mingote, & Weber, 2005). Spiny GABAergic neurons in the nucleus accumbens are involved in major depression (Francis & Lobo, 2016; Lim, Huang, Grueter, Rothwell, & Malenka, 2012; Wang et al., 2010). Compared with depression‐like mice, the mice with resilience to chronic stress demonstrate normal functions in neuronal spiking capability and synaptic transmission (Zhu et al., 2017). These data suggest that the nucleus accumbens is correlated to major depression versus resilience. It is necessarily to study molecular profiles in the nucleus accumbens underlying resilience to chronic stress versus major depression.

In present study, we have analyzed molecule profiles relevant to major depression and resilience in the nucleus accumbens from control mice as well as resilience versus depression mice induced by chronic unpredictable mild stress (CUMS), since mice were considered as a desirable system for gene analysis in the central nervous system to reveal molecules relevant to this disorder (Crowley & Lucki, 2005; Urani, Chourbaji, & Gass, 2005), and the different changes were clearly detected in the nucleus accumbens from CUMS‐induced depression and resilience mice (Zhu et al., 2017). The high‐throughput sequencings of miRNA and mRNA were performed in the nucleus accumbens from CUMS‐induced depression, resilience, and control mice. Through associated analyses and comparisons, we expect to figure out comprehensive molecular profiles in the nucleus accumbens relevant to stress‐induced depression versus resilience, in order to provide the guidelines for addressing molecular mechanisms underlying major depression and endogenous anti‐depression as well as for developing new therapeutic strategies.

MATERIALS AND METHODS

Mice were bred under well‐being condition and housed in cages (32 × 16 × 16 cm) with free access to food and water. Illumination was given between 7:00 am and 7:00 pm. Ambient temperature was kept at 22 ± 2°C. Relative humidity was maintained at 55 ± 5%. Their living condition was specific pathogen free (SPF). All experiments were conducted in accordance with the guidelines and regulations by the Administration Office of Laboratory Animals at Beijing China. Experiment protocols were approved by Institutional Animal Care and Use Committee in this office (B10831).

Mouse model for major depression induced by chronic unpredictable mild stress

Strain C57 male mice were selected for our experiments starting at postnatal day 21. In Week 1 for their adaptation, body weight, locomotion, sucrose preference, and Y‐maze tests (YMTs) were measured to have self‐control data. Mice that show consistent values in mean ± 2 SD of these measurements at postnatal Day 28 were separated into the groups of CUMS and control, in order to reduce variation among these mice. Control mice lived without mild stresses. The use of juvenile mice for our study is based on a fact that young individuals have high prevalence to suffer from major depression in response to chronic stress (Xu et al., 2016; Zhu et al., 2017).

As weaknesses in cognition, emotion, interaction skill, circadian, and stress response were risk factors for major depression (Southwick & Charney, 2012), the CUMS was used to mice for inducing major depression in the following principles. Mice lived in stress environment, challenged stress condition, and experienced defeat outcomes. Memory to these negative events induced by the CUMS drove them to feel cognitive and emotional disability and to fall into anhedonia and low self‐esteem (Wang & Cui, 2017, 2018). The paradigms for the CUMS included social isolation, tilted cage, empty cage, damp sawdust cage, restraint space, white noise, strobe light, and circadian disturbance (Ma, Guo, Xu, Cui, & Wang, 2016; Ma, Xu, et al., 2016; Xu et al., 2016; Zhu et al., 2017). Except for social isolation, such conditions were randomly used to treat mice in separation or combination manner each day (see Table one in Xu et al., 2016 for detail). The CUMS was sustained for 4 weeks until some mice showed anhedonia and low self‐esteem. Extreme stresses in a single pattern including learnt hopelessness, electrical shock, social defeat, or tail clamp, which led to fear memory, were not used in our study, since these protocols might induce the outcome similar to specific anxiety and posttraumatic stressful disorder (Xu et al., 2016; Zhu et al., 2017).

Whether CUMS‐treated mice in 4 weeks fell into anhedonia and low self‐esteem was examined in Days 57∼59. Sucrose preference test (SPT) has been used to estimate anhedonia. YMT has been used to evaluate a loss of interest to their partners. Forced swimming test (FST) has been used to assess their self‐esteem (Dellu, Mayo, Cherkaoui, Le Moal, & Simon, 1992; Duman, 2010; Overstreet, 2012; Xu et al., 2016). Mouse ingestions of 1% sucrose water versus pure water in 4 hr were used in the SPT, whose value was presented as the ratio of the ingested sucrose water to the ingested sucrose water plus pure water. The YMT was conducted by monitoring mouse stay in a special arm and other two arms. The end of this special arm included a female mouse (named as M‐arm). In the monitoring about 3 min, M‐arm stay time was presented by the ratio of stay time in M‐arm to that in all three arms. The FST was merited by recording immobile time in a water cylinder (10 cm in diameter and 19 cm in depth at 25 ± 1°C). In these tests, the SPT, YMT, and FST were given before and after the CUMS. Before the SPT, mice were deprived from food and water for 3 hr to drive their enthusiasm of drinking water. In the YMT, the arms were cleaned by 70% ethanol and then water after each test to reduce the effect of odor on the test. Carefulness in these tests was taken by performing them in a quiet room, no addition stress, same circadian circle for all mice, and their adaptation in the test environment.

Depression‐like behaviors were accounted if the CUMS mice demonstrated decreases in sucrose preference (twice at the ends of Week 3 and 4) and M‐maze stay time, as well as the increase in immobile time, compared to these values during their self‐control period (Week 1 for adaption) and in control mice. These measurements in each of mice would be considered as significant change if the values of the SPT and YMT reduced above 20% of its self‐control value and the immobile time of the FST increased above 15% of its self‐control value. These standards set up based on the averaged values in previous studies (Ma, Xu, et al., 2016; Xu et al., 2016; Zhu et al., 2017). Mice with significant changes in three tests were called as CUMS‐induced depression‐like mice or depression‐like mice. Mice with less than 5% changes in all three tests were named as CUMS‐resistance mice, that is, resilience mice. As showed in Table 1, CUMS‐treated mice in 4 weeks met this criterion about 27.54% (i.e., their vulnerability to chronic stress), and 18.84% of them were resilience (i.e., their invulnerability to chronic stress). These depression and resilience mice as well as the control mice were used to analyze microRNA and mRNA. In other words, tissues of the nucleus accumbens from control, CUMS‐induced depression, and CUMS‐Resilience mice were harvested for the high‐throughput sequencing of microRNA and mRNA. Mice with significant change in one or two measures are atypical in major depression, which are not included in our study.

Table 1

CUMS‐induced behavioral changes in the mice

	Number	Percentage (%)
Depression‐like mice	19	27.54
Resilient mice	13	18.84
Atypical mice	37	53.62
Total CUMS‐treated mice	69	100

RNA purification from the nucleus accumbens

Mice in groups of CUMS‐MDD, CUMS‐resilience, and control were anesthetized by using isoflurane, perfused by normal saline (4°C) through left atrium and decapitated by the guillotine. Both sides of the nucleus accumbens (NAc) were quickly isolated and dissected on ice‐cold glass slide. Total RNAs from the NAc in each mouse were isolated by using TRIzol Reagent (Life Technologies, Carlsbad, CA) as described (Ma, Xu, et al., 2016), after which RNA samples in the dry ice were delivered to Beijing Genomics Institute (BGI) in China for high‐throughput sequencing analysis. The concentration of total RNA, the value of RNA integrity number (RIN), and the ratios of 28S to 18S ribosomal RNA were measured by using 2100 Bioanalyzer (Agilent Technology) with RNA 6000 nano Reagents Port 1 for quality control of RNA samples. The samples with total RNA amount larger than 10 μg, the concentration larger than 200 ng/μL, the RIN larger than 8, and the ratio of 28S to 18S larger than 1.5 were selected for the construction of transcriptome and small RNA libraries (Ma, Guo, et al., 2016; Ma, Xu, et al., 2016). NAc tissues from three mice of control, CUMS‐MDD, and CUMS‐resilience groups were used for high‐throughput sequencing with correlation coefficient larger than 0.8.

RNA sequencing

mRNAs were extracted from total RNA and randomly sheared into 200 bp fragment by using oligo (dT) beads. Fragments were reversely transcribed to cDNA by random oligonucleotides. These synthesized cDNAs were purified by QiaQuick PCR extraction kit and ligated by sequencing adaptors after end repair. To select and purify cDNAs by agarose gel electrophoresis, the amplifications were done with Illumina PCR Primer Cocktail in 15 PCR‐reaction cycles. Finally, cDNAs with sizes between 200 and 300 bp were used for library construction.

Low molecular weight RNAs (18–30 nt) in total RNAs were isolated by polyacrylamide gel electrophoresis, which were used for constructing microRNA sequencing library. RNAs ligated to 5′‐RNA adapter by T4 RNA ligase were size‐fractionated and 36–50 nucleotide fractions were excised. The precipitated RNAs were ligated to 3′‐RNA adapter by T4 RNA ligase and size‐fractionated, and 62–75 nucleotide fraction (small RNA + adaptors) was excised. To produce sufficient templates for the sequencing, small RNAs ligated with adaptors were subjected to RT‐PCR, in which the products were purified and collected by gel purification and ready for the high‐throughput sequencing.

The qualities of transcriptome and microRNA libraries were assessed by using 2100 Bioanalyzer (Agilent Technologies, CA). Their quantities were verified by using ABI StepOnePlus Real‐Time PCR System. Their sequencings were performed by using Illumina HiseqTM 2500 platform (Illumina Inc., San Diego, CA). In two libraries, the average reading length were about 100 bp (pair‐end) and 49 bp (single‐end), respectively.

Bioinformatics for mRNA

Original image data was transferred into raw data or raw reads by base calling. Dynamic Trim Perl script implemented in SolexaQA package was done to control the quality of raw sequencing data according to the following criteria. The reads with adapters, the reads with unknown bases more than 10% as well as the reads with 50% of the bases and low quality score (PHRED score 5) were removed. The remained reads were called as “clean reads” and mapped to mouse genome reference sequence (UCSC mm10) by using TopHat v1.0.12 which incorporated Bowtie v0.11.3 software to perform the alignments. The maximum of allowable mismatch was set to three for each read in the alignment and mapping. To calculate gene expression level, the sole reads uniquely aligned to genes were used. The reads per kilo‐base per million reads (RPKM) were used for gene expression and the genes with low expression level (RPKM < 0.5) were removed for further analysis.

The differentially expressed genes (DEGs) were screened based on NOIseq package method that affected DEGs between two groups with the biological replicates, such as control versus CUMS‐induced major depressive disorder (CUMS‐MDD), control versus CUMS‐resilience, and so on. The threshold used to identify DEGs was fold‐change larger than 1.5, since the threshold at 2.0 might not be sensitively to detect changes in RNA profiles (Ma, Guo, et al., 2016; Ma, Xu, et al., 2016). The pathway enrichment analysis in DEGs association with physiological or biochemical processes were conducted. Subsequently, the hypergeometric test implemented in tool WebGestalt (version 2) and the canonical pathways from the Kyoto Encyclopedia of Genes and the Genomes (KEGG) database were used in these enrichment analyses. Compared with whole genome background, the enriched metabolic pathways or signal transduction pathways in DEGs would be identified in these analyses. p‐Values from hypergeometric tests were adjusted by Benjamini–Hochberg method, and pathways with adjusted p‐values less than 0.05 were considered to be significant enrichments.

Bioinformatics for microRNA (miRNA)

Adaptor sequences, low‐quality reads and contaminants in 49nt‐ sequence tags from Hiseq sequencing were initially removed to obtain credible clean reads. The reads were aligned to Genbank database and Rfam database with blast or bowtie software to remove reads from noncoding RNA, such as ribosomal RNAs, transfer RNAs, small nuclear RNAs, small nucleolar RNAs, and repeat RNA. To obtain miRNA count, high‐quality clean reads ranging from 18 to 25 nt were subsequently matched to the known miRNA precursor of corresponding species in miRBase. miRNAs satisfied by the following criteria, the tags aligned to miRNA precursor in miRBase with no any mismatch and the tags aligned to the mature miRNA in miRBase with at least 16 nt overlap allowing offsets would be counted to get the expression of identified miRNAs. For the remained reads without any annotation, Miredp was used to predict the potential novel miRNAs and its stem loop structure based on references (Friedländer, Mackowiak, Li, Chen, & Rajewsky, 2012; Friedländer et al., 2014). To correct the biased results from low expression, miRNAs with read counts less than 5 were discarded in the differential expression analysis.

DESeq software algorithm based on negative binomial distribution and biology duplicate samples was used to compare the known or novel miRNA expression in control versus CUMS‐MDD, control versus CUMS‐resilience. Threshold used to identify different expression of miRNAs was fold‐change larger than 1.5 and p‐value less than 0.05. Three miRNA target prediction software, for example, RNAhybrid, Targetscan, and miRanda, were used to predict gene targets of differentially expressed miRNAs.

Integrated miRNA/mRNA network analysis

A series of bioinformatics analyses were conducted to find the correlations between different expression of miRNAs and their target mRNAs. miRNAs were usually negatively correlated with their targeted mRNAs, in spite of a few exceptions (He et al., 2016). The datasets of differentially expressed miRNAs and transcripts were integrated to identify the potential miRNA‐regulated target genes based on the following criteria: (a) miRNAs and mRNAs should be simultaneously and reversely changed in our analyses; (b) mRNAs should be predicted by miRNAs from software of RNAhybrid, Targetscan, and miRanda. The interactive networks from differentially expressed miRNAs and simultaneously expressed target mRNAs were visualized by using Cytoscape software (San Diego, CA).

Quantitative RT‐PCR for the validations of miRNA and mRNA

Quantitative real‐time RT‐PCR (qRT‐PCR) by SYBR Green technique were used to analyze 16 mRNAs and 14 miRNAs that were involved in different cellular functions as well as were significantly difference among control (n = 13), resilience (n = 13), and depression‐like mice (n = 13). Supporting Information Table S1 listed the primers used in this study. Briefly, Real‐time PCR was performed with a Bio‐Rad CFX96Touch. Total RNA was extracted from the nucleus accumbens with a TRIzol Kit. cDNA was synthesized using the PrimeScript RT Reagent Kit (TaKaRa, RR037A, Kusatsu, Japan) and Mir‐X miRNA First‐Strand Synthesis Kit (Clontech, CA) for mRNA and microRNA, respectively. Amplification of mRNA was performed in a 20 μL reaction with 1 μL sample cDNA, 0.5 μL of each primer (10 nmol/L), 10 μL 2×qPCR Mastermix (Green), and 8 μL ddH2O. The real‐time PCR was initiated at 94°C for 2min, followed by 45 cycles of denaturation for 15 s at 94°C, annealing and elongation for 30 s at 60°C, and melt curve 65 to 95°C increment 0.5°C for 5 s. For miRNAs, the qRT‐PCR was performed in a 20 μL reaction with 1.6 μL sample cDNA, 0.4 μL mRQ3'Primer, 0.4 μL miRNA‐specific Prime (10μM), 10 μL 2×SYBR Advantage Premix, and 7.6 μL ddH2O, in which the program was set to 95°C for 10 s, followed by 40 cycles of denaturation for 5 s at 95°C, annealing and elongation for 20 s at 60°C, and melt curve 65 to 95°C increment 0.5°C 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. All qRT‐PCR runs were repeated in three replications. The results were calculated with the 2−ΔΔCt method.

Dual luciferase reporter assay

The sequence containing targeted sites of targeted gene was amplified, digested by XhoI and NotI and fused into the luciferase vector psiCHECK2 (Ma, Guo, et al., 2016; Ma, Xu, et al., 2016). The site‐directed mutation of the detected miRNA targeting site was carried out with QuikChange Lighting Site‐Directed Mutagenesis Kit (Stratagene, La Jolla, CA) according to the manufacturer's instructions. Luciferase reporter detection assays were performed as previously described (Ma, Xu, et al., 2016). Briefly, HEK293T cells were planted at 5 × 104 cells per well in 24‐well plates and maintained in DMEM containing 10% FBS. After 24 hr, these cells were co‐transfected with 50 ng psiCHECK2‐wild‐type or mutant reporter plasmids, 50 nM miRNAs mimic or miR‐NC by using Lipofectamine 2000 transfection reagent (Invitrogen, Carlsbad, CA). After 48 hr, the activities of firefly and Renilla luciferase were assessed by using Dual‐Glo Luciferase Assay System (Promega, Cat. E2920) according to the manufacturer's instructions. Each treatment was carried out in the triplicates in three independent experiments.

Differential expression and biological indication

In the comparisons between control and CUMS‐MDD, between control and CUMS‐resilience as well as between CUMS‐MDD and CUMS‐resilience groups in the nucleus accumbens (NAc), the 1.5‐fold ratio of RPKM values in differential expressions was set as the involvement of CUMS, MDD or resilience. If some mRNAs were differentially expressed between control and CUMS‐MDD mice, these mRNAs were presumably associated to CUMS treatment or CUMS‐MDD. mRNAs that were differentially expressed between control and CUMS‐resilience were presumably associated with CUMS treatment or resilience. mRNAs that were differentially expressed between CUMS‐MDD and CUMS‐resilience were likely associated with CUMS‐MDD versus CUMS‐resilience. Moreover, if some mRNAs were differentially expressed between control and CUMS‐MDD as well as between control and CUMS‐resilience, but showed similar expression between CUMS‐MDD and CUMS‐resilience, these mRNAs were likely associated with CUMS treatment.

Statistical analyses

Using NOIseq and DESeq software algorithm, initial processing raw data of mRNA and miRNA expression profiles were done, respectively. The data of the behavior tests, luciferase activity, and gene analyses are presented as mean ± SEM. Relationships between miRNA and its target prediction were assessed by Pearson's correlation coefficients. The unpaired Student t test was used to make the statistic comparison between control and CUMS‐MDD, control and CUMS‐resilience, CUMS‐MDD and CUMS‐resilience, and so on. p < .05 is considered statistically significant (Ma, Guo, et al., 2016; Ma, Xu, et al., 2016).

RESULTS

Chronic unpredictable mild stress to mice leads to either depression‐like behaviors or resilience

Mice were treated by CUMS or control in 4 weeks. Their mood was assessed by the SPT, YMT, and FST. CUMS‐treated mice appear significant change in all three tests (depression mice), no change in anyone (resilience), or changes in one or two of these tests (atypical). Figure 1 shows statistical analysis for mice with major depression (blue bars), resilience (reds), and control (whites). Compared with control mice (n = 7), the SPT values in CUMS‐MDD mice are 68.678 ± 4.129% after the CUMS and 95.169 ± 1.791% before the CUMS (p < .001, n = 8, paired t test; Figure 1a). The SPT values in CUMS‐resilience mice are 92.04 ± 1.03% and 91.76 ± 1.56% after and before the CUMS (n = 12; Figure 1a), respectively. In comparison with the YMT values in control mice (n = 6), the ratios of stay time in the M‐arm to stay time in total arms from CUMS‐MDD mice are 48.260 ± 3.303% after the CUMS and 80.226 ± 2.005% before the CUMS (p < .001, n = 12, paired t test; Figure 1b). These ratios in CUMS‐resilience mice are 74.98 ± 2.67% after the CUMS and 74.59 ± 1.62% before the CUMS (n = 12, Figure 1b). Compared with control mice (n = 5), the values of immobile time in the FST in CUMS‐MDD mice are 266.00 ± 13.18% after the CUMS and 167.80 ± 13.36% before the CUMS (p < .001, n = 5, paired t test; Figure 1c). These values in CUMS‐resilience mice are 193.25 ± 10.54% after the CUMS and 166.750 ± 16.439% before the CUMS (n = 5, Figure 1c). CUMS leads to depression‐like behaviors in certain mice versus resilience in others.

Figure 1

Chronic unpredictable mild stress (CUMS) leads mice to express depression‐like behaviors or resilience. (a) shows the SPT values (%) in the mice from CUMS‐MDD (blue bar), CUMS‐Resilience (red bar), and control group (white bar). (b) illustrates the ratios of stay time in M‐arm to stay time in three arms by the YMT in the mice from CUMS‐MDD (blue bar), CUMS‐Resilience (red bar), and control group (white bar). (c) illustrate immobile time of staying in the water cylinder by the FST in the mice from CUMS‐MDD (blue bar), CUMS‐Resilience (red bar), and control group (white bar). Three asterisks show p < .001, in which one‐way ANOVA was used for the comparisons among CUMS‐MDD, CUMS‐Resilience, and control mice, while paired‐t test was for the comparisons before and after the CUMS [Color figure can be viewed at http://wileyonlinelibrary.com]

Taken our data with the fact that acute and chronic stresses influence the nucleus accumbens (NAc) to induce a depressed mood (Walsh & Han, 2014), we have examined molecular mechanisms underlying major depression and resilience by analyzing miRNAs and mRNAs with high‐throughput sequencing to quantify their expression levels in the NAc from CUMS‐MDD, CUMS‐resilience, and control mice.

Overall qualities of RNA‐sequencing dataset

High‐throughput RNA‐sequencing was applied for transcriptome and microRNA profiles in each sample of the NAc from CUMS‐MDD, CUMS‐resilience, and control mice. RNAs of 53.9 Mb to 60.43 Mb raw sequence reads about 100 bp were obtained from mRNA library Illumina sequencing (three control libraries, three CUMS‐MDD libraries and three CUMS‐resilience libraries). After filtering the reads that contained N with adaptor sequence and low quality, 44.15 Mb to 44.98 Mb clean reads from each library were generated and mapped, which were about 84.93∼90.38% of total reads from the mouse genome (UCSC mm10) equivalently for all samples (Supporting Information Table S2). In addition, the totals about 29,045,586∼30,896,262 raw tag counts were produced in small RNA library. After filtering reads with low quality and adaptor, the totals of clean tag counts about 24,883,900∼28,355,107 were obtained (Supporting Information Table S3). The distribution of clean small RNA reads varied in a range of 10 to 44 nucleotides each library. The most abundant lengths were 22 nucleotides (Supporting Information Figure S1). All high‐quality clean reads larger than 18 nucleotides were mapped to the mouse genome. Genome‐matched reads were divided into different categories of small RNAs based on their biogenesis and annotation (Supporting Information Figure S2). The most abundant RNA category from each library was miRNA. The high qualities of transcriptome and small RNA sequencing data were used for further analysis.

mRNA differential expression in NAc among CUMS‐MDD, CUMS‐resilience, and control mice

mRNAs in the NAc were quantified by sequencing total RNAs. Their RPKM values were calculated. Genes in lower expression level (RPKM < 0.5) were removed. Others were left for differential expression analyses by NOISeq. After mapping reads referred to the mouse genome, these 18,534 (18,548) mRNAs from clean read sequences with high quality included 9316 (9308) of upregulated mRNAs and 9218 (9240) of downregulated mRNAs by comparing control versus CUMS‐MDD and control versus CUMS‐resilience. The criterion to make sure differential expression of genes was the ratio of their RPKM values above 1.5‐fold (Biggar & Storey, 2011; Dwivedi et al., 2015; Ma, Guo, et al., 2016; Ma, Xu, et al., 2016). If the ratio of their expression altered above 1.5‐fold between control and CUMS‐MDD, between control and CUMS‐resilience as well as between CUMS‐MDD and CUMS‐resilience (probability ≥ 0.8), the differential expression of mRNAs was warranted.

Taken out the genes with low expression (i.e., the single digit of RPKM value), we found that 120 mRNAs had the 1.5‐fold ratio of CUMS‐MDD to control in RPKM value, where 61 mRNAs were downregulated and 59 mRNAs were upregulated (Supporting Information Table S4). Their involved signaling pathways are included in Table 2 and Supporting Information Table S5. The lowered expressions of mRNAs in the NAc from CUMS‐MDD mice mainly include Plin4, Cckar, Hif3a, Kcnj5, Per2, Hnf1b, Sec31b, Zbtb16, Svep1, and so on. Based on bioinformatics for mRNA‐guided protein translation (KEGG database), the downregulated signaling pathways in CUMS‐MDD mice are dopaminergic synapses, circadian entrainment, and calcium signaling pathway (Table 2). In contrast, the raised expressions of mRNAs in the NAc from CUMS‐MDD mice mainly include Alox12, Tshr, CD40, Tcaf2, Pcdhb4, Pcdh12, Opn1mw, and so on. According to bioinformatics for mRNA‐guided protein translation (KEGG), the upregulated signaling pathways in CUMS‐MDD include tyrosine metabolism, regulation of actin cytoskeleton, and inflammatory mediator regulation of TRP channels (Table 2). It is noteworthy that the expression of some genes in a given signaling pathway is upregulated and others in this pathway are downregulated from CUMS‐MDD mice, such as serotonergic synapse, neuroactive ligand‐receptor interaction, and MAPK signaling pathway (Table 2). This up‐ and down‐expression of genes in the given pathways indicates the imbalance of molecules in these pathways. In addition to the imbalance among genes in a single signaling pathway, the upregulated and downregulated multiple signaling pathways make entire molecular networks in the NAc to be imbalance, which lead to neuronal dysfunction in the NAc for major depressive disorder.

Table 2

Signaling pathways identified by KEGG function analysis based on DEGs data in control versus CUMS‐induced depression mice

Pathway	DEGs with pathway annotation (120)	All genes with pathway annotation (18534)	Symbol	Pathway ID
Neuroactive ligand‐receptor interaction	2（1.67%）	247(1.33%)	Cckar↓, Tshr↑	ko04080
Serotonergic synapse	2（1.67%）	109(0.59%)	Kcnj5↓, Alox12↑	ko04726
Dopaminergic synapse	1（0.83%）	123(0.66%)	Kcnj5↓	ko04728
Inflammatory mediator regulation of TRP channels	1（0.83%）	100(0.54%)	Alox12↑	ko04750
Regulation of actin cytoskeleton	2（1.67%）	202(1.09%)	Itgal↑, Myl2↑	ko04810
Morphine addiction	1（0.83%）	88(0.47%)	Kcnj5↓	ko05032
Glycine, serine, and threonine metabolism	1（0.83%）	37(0.20%)	Alas2↑	ko00260
Tyrosine metabolism	2（1.67%）	36(0.19%)	Gm21948↑, Adh1↑	ko00350
MAPK signaling pathway	3（2.50%）	246(1.33%)	Map3k6↓, Dusp9↑, Fas↓	ko04010
Calcium signaling pathway	1（0.83%）	176(0.95%)	Cckar↓	ko04020
cAMP signaling pathway	1（0.83%）	183(0.99%)	Tshr↑	ko04024
Cytokine‐cytokine receptor interaction	3（2.50%）	186(1.00%)	Figf↑, Cd40↑, Fas↓	ko04060
Chemokine signaling pathway	1（0.83%）	166(0.90%)	Ccl9↓	ko04062
PI3K‐Akt signaling pathway	1（0.83%）	313(1.69%)	Figf↑,	ko04151
Apoptosis	1（0.83%）	58(0.31%)	Fas↓	ko04210
Cell adhesion molecules (CAMs)	3（2.50%）	139(0.75%)	Itgal↑, Cd6↑, Cd40↑	ko04514
Parkinson's disease	1（0.83%）	125(0.67%)	Uba1y↑	ko05012
Autoimmune thyroid disease	3（2.50%）	42(0.23%)	Cd40↑, Tshr↑, Fas↓	ko05320
Circadian entrainment	2 (1.6667%）	89 (0.4803%)	Kcnj5↓, Per2↓	ko04713

Note.↑ indicates upregulation in the tissue of nucleus accumbens from CUMS‐MDD versus control mice, whereas ↓ represents downregulation.

Furthermore, our study showed that 57 mRNAs had the 1.5‐fold ratio of CUMS‐resilience to control in RPKM values, in which 30 mRNAs are downregulated and 27 mRNAs are upregulated (Supporting Information Table S6). The involved signaling pathways are listed in Table 3 and Supporting Information Table S7. The decreased expressions of mRNAs in the NAc of CUMS‐resilience mice include Rtn4r, Nrn1, Slc17a6, Zbtb16, and so on. Based on the bioinformatics of mRNA‐guided protein translation (KEGG), the downregulated signal pathways in CUMS‐Resilience include synaptic vesicle cycle, glutamatergic synapse, nicotine addiction, and thyroid hormone signaling pathway (Table 3 and Supporting Information Table S7). Conversely, the raised expressions of mRNAs in the NAc from CUMS‐resilience mice mainly include Edar2, Adh1, Gna14, Ttr, Xist, Twist1, and so on. Based on the bioinformatics of mRNA‐guided protein translation (KEGG), the upregulated signal pathways in CUMS‐resilience include cytokine‐cytokine receptor interaction, tyrosine metabolism, and calcium signaling pathway (Table 3 and Supporting Information Table S7). Interestingly, compared to CUMS‐MDD, there is only a metabolic pathway showing that some genes express increasingly and others express decreasingly in CUMS‐resilience. Thus, compared with CUMS‐MDD, less changes of signaling pathways in CUMS‐resilience.

Table 3

Signaling pathways identified by KEGG function analysis based on DEGs data in control versus CUMS‐Resilience mice

Pathway	DEGs with pathway annotation (57)	All genes with pathway annotation (18548)	Contributing Genes	Pathway ID
Tyrosine metabolism	1 (1.75%)	34 (0.18%)	Adh1↑	ko00350
ErbB signaling pathway	1 (1.75%)	86 (0.46%)	Areg↑	ko04012
Calcium signaling pathway	1 (1.75%)	176 (0.95%)	Gna14↑	ko04020
Cytokine‐cytokine receptor interaction	1 (1.75%)	192 (1.04%)	Eda2r↑	ko04060
Chemokine signaling pathway	1 (1.75%)	169 (0.91%)	Grk1↓	ko04062
Synaptic vesicle cycle	1 (1.75%)	61 (0.33%)	Slc17a6↓	ko04721
Retrograde endocannabinoid signaling	1 (1.75%)	98 (0.53%)	Slc17a6↓	ko04723
Glutamatergic synapse	1 (1.75%)	108 (0.58%)	Slc17a6↓	ko04724
Phototransduction	1 (1.75%)	20 (0.11%)	Grk1↓	ko04744
Protein digestion and absorption	1 (1.75%)	74 (0.40%)	Slc16a10↓	ko04974
Nicotine addiction	1 (1.75%)	39 (0.21%)	Slc17a6↓	ko05033

Note. ↑ indicates upregulation in the tissue of nucleus accumbens from CUMS‐Resilience versus control mice, whereas ↓ represents downregulation.

In order to validate our data above, we ran quantitative RT‐PCR (qRT‐PCR) from tissues that had been used for mRNA sequencing. The expressions of Myl2, Itgal, Alox12, Cd40, and Tshr are raised, as well as the expressions of Cckar, Fas, Hnf1b, Kcnj5, Sec31b, and Per2 are decreased in CUMS‐MDD mice, compared to control mice (Figure 2). The expressions of Twist and Xist are raised, as well as the expressions of Rtn4r, Nrn1, and Slc17a6 are decreased in CUMS‐resilience mice, compared with control mice (Figure 3). Consistent results achieved by mRNA sequencing and qRT‐PCR confirm the validation of our study.

Figure 2

The validation of differentially expressed mRNAs in the NAc from control mice versus CUMS‐MDD mice. Three asterisks show p < .001, two asterisks show p < .01, one asterisk show p < .05, in which unpaired‐t test was used for the comparisons between CUMS‐MDD and control mice [Color figure can be viewed at http://wileyonlinelibrary.com]

Figure 3

The validation of differentially expressed mRNAs in the NAc from control mice versus CUMS‐Resilience mice. Three asterisks show p < .001, two asterisks show p < .01, in which unpaired t test was used for the comparisons between CUMS‐Resilience and control mice [Color figure can be viewed at http://wileyonlinelibrary.com]

By comparing changed mRNAs above 1.5‐fold among groups of control, CUMS‐MDD and CUMS‐resilience in their separation and overlapping, we find that the expression changes of certain genes are specifically involved in CUMS‐MDD or CUMS‐resilience (Figure 4), except for Adh1, Ptgis, Tecrl, Ttr, Zbtb16, 4932441J04Rik, and B930018H19Rik in both groups. These results will help us to figure out those genes related to major depression or resilience. The genes changed in both CUMS‐MDD and CUMS‐resilience mice are likely involved in the CUMS treatment.

Figure 4

Venn diagram showed the differently expression genes (DEGs) between control versus CUMS‐MDD and control versus CUMS‐Resilience. MDD (blue) shows the certain genes specifically involved in CUMS‐induced depression. Resilience (red) shows the certain genes specifically involved in CUMS‐Resilience. CUMS (yellow) indicates the genes changed in both CUMS‐MDD and CUMS‐Resilience mice, which are likely involved in the CUMS treatment [Color figure can be viewed at http://wileyonlinelibrary.com]

In principle, the level of mRNAs in cells is affected by miRNAs, through which the bindings of miRNAs with their dicers degrade mRNAs and weaken their translations (Afonso‐Grunz & Muller, 2015; Beilharz, Humphreys, & Preiss, 2010; Dalmay, 2013; Valinezhad Orang, Safaralizadeh, & Kazemzadeh‐Bavili, 2014). If mRNAs downregulated in the NAc from CUMS‐MDD mice are caused by miRNAs, their correspondent miRNAs will be upregulated, or vice versa. If mRNAs upregulated in the NAc from CUMS‐resilience mice are caused by miRNAs, correspondent miRNAs will be downregulated, or vice versa. To test this hypothesis and validate our data about mRNA changes, we analyzed the changes of miRNAs by their sequencings in the NAc from control, CUMS‐MDD, and CUMS‐resilience mice.

Changes of miRNA expression in the NAc among control, CUMS‐MDD, and CUMS‐resilience mice

The expression profile of miRNAs is presented in Supporting Information Table S8 if their expressions change above 1.5‐fold in all of CUMS‐MDD mice versus control, in which some miRNAs are upregulated or downregulated. Their predicted target mRNAs match the measures by mRNA sequencing based on the databases (RNAhybrid, Targetscan, and miRanda) about complex interactions between miRNAs and mRNAs. Supporting Information Table S9 shows the altered miRNAs and their predicted‐target mRNAs. Supporting Information Table S10 shows the altered mRNAs and their correspondent miRNAs. Interactive networks about the miRNAs and overlapped mRNAs, which were based on transcriptome expression data and predicted target genes from three databases, were made in the Cytoscape (Figure 5). By reading Supporting Information Tables S9–S11, we can find that the regulations of mRNAs and the regulations of miRNAs are matched well.

Figure 5

MicroRNA‐mRNA network in control versus CUMS‐MDD mice. microRNA/mRNA networks were constructed between the 162 miRNAs and 75 overlapped mRNAs with using transcriptome expression data and predicted target genes from RNAhybrid, Targetscan, and miRanda databases. Blue symbols present the elevated expression of miRNAs or mRNAs. Red symbols present the down regulated miRNAs or mRNAs [Color figure can be viewed at http://wileyonlinelibrary.com]

The expression profile of miRNAs is demonstrated in Supporting Information Table S12 if their expressions change above 1.5‐fold in all CUMS‐resilience mice versus control, in which certain miRNAs are upregulated or downregulated. Their predicted target mRNAs match the measures by mRNA sequencing based on RNAhybrid, Targetscan, and miRanda about the complex interactions between miRNAs and mRNAs. Supporting Information Table S13 shows the altered miRNAs and their predicted‐target mRNAs. Supporting Information Table S14 shows the altered mRNAs and their correspondent miRNAs. Interactive networks about miRNAs and overlapped mRNAs, which were based on transcriptome expression data and predicted target genes from three databases, were made in the Cytoscape (Figure 6). Based on Supporting Information Tables S13–S15, the regulations of mRNAs and the regulations of miRNAs are well matched. Consistent results by jointly sequencing mRNAs and miRNAs validate our analyses and strengthen our conclusion.

Figure 6

MicroRNA‐mRNA network in control versus CUMS‐Resilience mice. microRNA/mRNA networks were constructed between the 82 miRNAs and 41 overlapped mRNAs with using transcriptome expression data and predicted target genes from RNAhybrid, Targetscan, and miRanda databases. Blue symbols present the elevated expression of miRNAs or mRNAs. Red symbols present the down regulated miRNAs or mRNAs [Color figure can be viewed at http://wileyonlinelibrary.com]

In order to validate the finding by sequencing miRNA analysis, we selected certain miRNAs to do qRT‐PCR, including miR‐96‐5p, miR‐142a‐3p, miR‐384‐5p, miR‐770‐3p, miR‐667‐3p, miR‐540‐5p, miR‐874‐3p, and miR‐151‐3p in CUMS‐MDD versus control mice as well as miR‐15b‐5p, miR‐384‐5p, miR‐96‐5p, miR‐874‐3p, miR‐540‐5p, and miR‐449a‐5p in CUMS‐resilience versus control mice. Consistent to high‐throughput sequencing, these miRNAs are significantly changed in qRT‐PCR from CUMS‐MDD versus control mice as well as CUMS‐ resilience versus control mice (Figure 7). Their predicted target mRNAs match the actually altered miRNAs (Supporting Information Tables S9–S11 and S13–S15). Consistent results from analyses by miRNA sequencing and qRT‐PCR also support the validation of our study. It is noteworthy that certain miRNAs reversely express in CUMS‐MDD mice versus CUMS‐resilience mice, such as let‐7c‐2–3p, miR‐199a‐3p, miR‐199b‐3p, miR‐202‐5p, miR‐224‐5p, miR‐3074‐5p, miR‐3105‐5p, miR‐34b‐3p, miR‐34b‐5p, miR‐669c‐5p, miR‐7019‐3p, miR‐466c‐3p, and miR‐551b‐3p. In other words, these miRNAs may influence mice in the CUMS to be either depression or resilience (Supporting Information Tables S8 and S12).

Figure 7

The validation of differentially expressed microRNAs in the NAc. (a) The validation of differentially expressed microRNAs in the NAc from Control mice versus CUMS‐MDD mice. (b) The validation of differentially expressed microRNAs in the NAc from Control mice versus CUMS‐Resilience mice [Color figure can be viewed at http://wileyonlinelibrary.com]

Cckar and Rtn4r mRNA are the targets of miRNA‐384‐5p

To validate silico prediction (Supporting Information Tables S9, S10, S13, and S14), we selected miRNA‐384‐5p to examine whether it targeted to Cckar and Rtn4r by qRT‐PCR and dual luciferase reporter assay. There are inverse correlations between them from qRT‐PCR analysis (Figure 8). In dual luciferase report assay, we constructed luciferase reporter plasmids, which contained the wild‐type or mutant of the predicted binding sites of the miRNAs in Cckar and Rtn4r. These reporter constructs were transfected into HEK293T cells. The relative activities of luciferase reporter for the Cckar and Rtn4r mRNA are significantly lowered by the mimics of miRNA‐384‐5p, but not negative control (Figure 8), which are reversed by mutating the binding sites of miRNA‐384‐5p. These results support that Cckar and Rtn4r mRNA are the direct targets of miR‐384‐5p, which is consistent with our bioinformatics analyses in the prediction of miRNA target genes.

Figure 8

The miRNAs targeted mRNAs were validated by qRT‐PCR and Luciferase reporter assay. Correlation between miRNAs and its prediction target expression by qRT‐PCR in NAc tissue. (a) The correlation between Cckar and miR‐384‐5p (r = −0.9893; p < .001). (b) Luciferase reporter assay performed by the co‐transfection of luciferase reporter containing wild or mutant of Cckar mRNA with miR‐384‐5p mimic or negative control (NC) into HEK293T cells. (c) The correlation between Rtn4r and miR‐384‐5p (r = −0.995; p < .001). (d) Luciferase reporter assay performed by the co‐transfection of luciferase reporter containing wild or mutant of Rtn4r mRNA with miR‐384‐5p mimic or negative control (NC) into HEK293T cells. Data are the mean ± SEM [Color figure can be viewed at http://wileyonlinelibrary.com]

DISCUSSION

By high‐throughput sequencings of mRNAs and miRNAs, we have analyzed quantitative changes of mRNAs and miRNAs in the nucleus accumbens from control, CUMS‐MDD, and CUMS‐resilience mice. The downregulated mRNAs in CUMS‐MDD mice are those genes that encode neural processes about dopaminergic synapse, circadian entrainment and calcium signaling pathway, whereas the upregulated mRNAs are those genes that encode tyrosine metabolism, regulation of actin cytoskeleton, and inflammatory mediator regulation of TRP channels. To the nucleus accumbens in CUMS‐resilience mice, the downregulated mRNAs are those genes that encode neural events about synaptic vesicle cycle, glutamatergic synapse, nicotine addiction and thyroid hormone signaling pathway, whereas the upregulated mRNAs are those genes that encode neural processes about cytokine–cytokine receptor interaction, tyrosine metabolism and calcium signaling pathway. Our data suggest that the impairment of dopaminergic synapse, circadian process, and calcium signaling in the nucleus accumbens may be involved in CUMS‐induced depression, whereas the upregulation of tyrosine metabolism, cytokine‐cytokine receptor interaction, and calcium signaling in the nucleus accumbens may be involved in mouse resilience to the CUMS.

In addition to the downregulation or the upregulation of all mRNAs in some signaling pathways, the imbalanced expression of mRNAs in other signaling pathways may play role in major depression versus resilience. As indicated in Table 2, some mRNAs in the given signaling pathways may increase, whereas others in these pathways decrease in the nucleus accumbens of CUMS‐induced depression mice, such as serotonergic synapse, neuroactive ligand‐receptor interaction, and MAPK signaling pathway, indicating the imbalance of molecular networks in these pathways. In comparison with CUMS‐induced depression mice, there is only a metabolic pathway that shows the increase of one gene and the decrease of another one in CUMS‐resilience. Therefore, the imbalanced expression among mRNAs in a single signaling pathway may be responsible for major depression, whereas less bidirectional change of genes in intra‐signaling pathways may help mice to be CUMS‐resilience, that is, the prevention of mice from CUMS‐MDD. Therefore, the upregulation and downregulation among multiple signaling pathways as well as among mRNAs in intra‐signaling pathway make molecular networks in the nucleus accumbens to be imbalance, which lead to neuronal dysfunction in the nucleus accumbens for major depressive disorder.mRNAs and their translated proteins are presumably related to the CUMS if their expressions alter in the same direction in CUMS‐MDD and CUMS‐resilience mice, compared to control. This hypothesis is based on the fact that mice receive the CUMS, but their consequences differ in emotion events. By comparing the altered mRNAs among control, CUMS‐resilience and CUMS‐MDD mice in their separation and overlapping, we find that some genes are specifically changed in CUMS‐MDD or CUMS‐resilience (Figure 4), However, Adh1, Ptgis, Tecrl, Ttr, Zbtb16, 4932441J04Rik, and B930018H19Rik are changed in both groups. Adh1/Ptgis and their translated proteins influence tyrosine metabolism. Ttr/Zbtb16 and their translated proteins are involved in nerve regeneration (Fleming, Saraiva, & Sousa, 2007) and mental retardation (Fischer et al., 2008) (Tables 2 and 3). These neural events may be influenced by the CUMS. Therefore, these data help us to figure out mRNAs related to major depression and resilience as well. mRNAs changed in the same direction at the nucleus accumbens of CUMS‐induced major depression and CUMS‐resilience mice are likely involved in CUMS treatment. It is noteworthy that our study is focused on expressional changes in the level of mRNAs, which are presumably influenced by the epigenetic regulation of miRNAs. In fact, the changes in mRNA may be regulated by changes in genetic level and DNA methylation induced by chronic stress. These possibilities will be taken into our future studies.

Relationships between signaling pathways versus neural impairment or resilience to the CUMS are explained below. The CUMS‐induced downregulation of dopaminergic synapse, circadian and calcium signaling pathway in the nucleus accumbens from CUMS‐induced depression mice may weaken the functions of reward circuit, circadian and neuronal responsiveness. As the nucleus accumbens is thought as one of targets of dopaminergic neurons in the reward circuit for positive emotion, the impaired neurons, and dopaminergic synapses in the nucleus accumbens lead to negative mood in major depressive disorder. Conversely, the upregulations of tyrosine metabolism, calcium signaling pathway and cytokine‐cytokine receptor interaction will benefit dopamine synthesis and neuron responses in the nucleus accumbens for positive emotion in resilience mice, which counterbalances CUMS‐induced depressed mood.

In terms of the validation of our study and result, we have done the high‐throughput sequencing of mRNAs and miRNAs, the quantitative RT‐PCR of some RNAs changed in major depression or resilience, as well as the analysis of the interaction between mRNAs and mRNA. Our results indicate that the changed expression of mRNAs matches the changed expression of miRNAs well in high‐throughput sequencing (Supporting Information Tables S9, S10, S13, and S14). Some genes with their altered expressions in the high‐throughput sequencing have been confirmed by quantitative RT‐PCR analysis (Figures 2 and 3). Moreover, interactions between Cckar/Rtn4r mRNAs and miRNA‐384‐5p are confirmed by dual luciferase report analysis. Taken these analyses together, we are confident to our results, which is better than previous analyses in either miRNAs or mRNAs, and so on.

A few advantages in our study are presented below. We have paid attention to analyze and compare the profiles of mRNA and miRNAs in the nucleus accumbens of CUMS‐induced depression and CUMS‐resilience mice, which helps to figure out the molecules that are involved in depression, resilience or stress treatment. Second, the nucleus accumbens is an area that is believed to be involved in the emotional reaction to reward. If major depression is due to lack of reward, the analysis of molecular profiles in the nucleus accumbens from CUMS‐induced depression and resilience should help to figure out the role of the nucleus accumbens in major depression and resilience based on molecules analyzed. The alternations of certain genes have not been observed in previous analyses (Bai et al., 2014; Elizalde et al., 2010; Li et al., 2013; Liu et al., 2015; Moreau, Bruse, David‐Rus, Buyske, & Brzustowicz, 2011; Rajkowska et al., 2015; Smalheiser et al., 2012, 2011, 2014), including Cckar, Tshr, Kcnj5, Alox12, Slc17a6, Per2, Rtn4r, and Nrn1. These genes in the nucleus accumbens may be newly working for major depression or resilience.

AUTHORS CONTRIBUTIONS

Y. Si, X. Sun, and Z. Song contributed to experiments and data lanalyses. Jin‐Hui Wang contributed to concept, project design, and paper writing.

COMPETING INTERESTS

All authors declare no competing interest. All authors have read and approved the final version of the manuscript.

Appendix S1.

Click here for additional data file.