Identification of key genes, pathways, and miRNA/mRNA regulatory networks of CUMS-induced depression in nucleus accumbens by integrated bioinformatics analysis

Abstract

Introduction

Major depressive disorder (MDD) is a recurrent, devastating mental disorder, which affects >350 million people worldwide, and exerts substantial public health and financial costs to society. Thus, there is a significant need to discover innovative therapeutics to treat depression efficiently. Stress-induced dysfunction in the subtype of neuronal cells and the change of synaptic plasticity and structural plasticity of nucleus accumbens (NAc) are implicated in depression symptomology. However, the molecular and epigenetic mechanisms and stresses to the NAc pathological changes in depression remain elusive.

Materials and methods

In this study, treatment group mice were treated continually with the chronic unpredictable mild stress (CUMS) until expression of depression-like behaviors were found. Depression was confirmed with sucrose preference, novelty-suppressed feeding, forced swimming, and tail suspension tests. We applied high-throughput RNA sequencing to assess microRNA expression and transcriptional profiles in the NAc tissue from depression-like behaviors mice and control mice. The regulatory network of miRNAs/mRNAs was constructed based on the high-throughput RNA sequence and bioinformatics software predictions.

Results

A total of 17 miRNAs and 10 mRNAs were significantly upregulated in the NAc of CUMS-induced mice with depression-like behaviors, and 12 miRNAs and 29 mRNAs were downregulated. A series of bioinformatics analyses showed that these altered miRNAs predicted target mRNA and differentially expressed mRNAs were significantly enriched in the MAPK signaling pathway, GABAergic synapse, dopaminergic synapse, cytokine–cytokine receptor interaction, axon guidance, regulation of autophagy, and so on. Furthermore, dual luciferase report assay and qRT-PCR results validated the miRNA/mRNA regulatory network.

Conclusion

The deteriorations of GABAergic synapses, dopaminergic synapses, neurotransmitter synthesis, as well as autophagy-associated apoptotic pathway are associated with the molecular pathological mechanism of CUMS-induced depression.

Keywords:

• stress
• depression
• nucleus accumbens
• miRNA
• mRNA

Introduction

Major depressive disorder (MDD) is a common, diverse, widespread, and intermittent neuropsychiatric burden that includes the wide range of symptoms, such as inability to feel pleasure, lack of motivation, retraction from social interaction, cognitive difficulties, and changes in appetite.^1,2 Depression affects a significant portion of the world’s population annually; according to the World Health Organization, it is the leading cause of disability in the world.³ Extensive research efforts have been done in past decades; however, the etiology of depression is still indefinable, its diagnosis unclear, and the pharmacotherapy ineffective. This may be due to poor understanding of the molecular pathophysiology of depression.⁴

It has been established that depression is a collection of diseases, with overlapping causal pathways, beginning with interaction among genetic and environmental factors.^5,6 The persistent stress to the genetically susceptible persons directs to the deficits of hypothalamus–pituitary–adrenal axis, neurotrophic factors, monoamine factor, and cytokine, which induce neuron atrophy and dysfunction in brain reward circuits, including the nucleus accumbens (NAc), prefrontal cortex (PFC), ventral tegmental area (VTA), and basolateral amygdala, in the depression patients and animal models with depression-like behaviors.^7,8 A number of abnormally regulated genes have been reported in genetically susceptible persons.^9–11 However, existing data cannot explain mechanism that mediates pathological changes in depression. Thus, it is important to investigate the potential changes in gene expressions and the regulation of these changes, specifically in brain regions that mediate emotional memory and anxiety.

It has been suggested that NAc is essential for emotion processing and receives dense innervations from mood-related structures including PFC, VTA, hippocampus, amygdala, and hypothalamus.^12,13 It is regarded as a neural interface between action and motivation, and has an important role in food intake, sexual behavior, reward-motivated behavior, and stress-related behavior.¹⁴ The dysfunction in this region may be correlated to low motivation, interest loss, and anhedonia in depression disorder. Recently, various reports have suggested that synaptic plasticity, structural plasticity changes in NAc medium spiny neurons, and reduced activity are implicated in depression symptomology.^15–18 Recently, neuromodulation, an intervention targeting the NAc, has been used in patients diagnosed with treatment-resistant depression.^19,20 Furthermore, NAc deep brain motivation has also been linked with anxiolytic and antidepressant effect and enhancement of quality of life in patients diagnosed with severe treatment-resistant depression.²¹

Gene expression disturbance in the brain is most probably linked with environmental factors which lead to depression.^22,23 This assumption has been established by analyzing miRNA in depressive patients^24,25 and in mice with depression-like behaviors.²⁶ A reaction chain from stress to miRNA dysregulation, alternations in mRNA/protein expression, and neuron atrophy and function deficit are proposed to be associated with MDD. However, the molecules acknowledged for major depression are exhibited differently in different reports. The inconsistency may be due to different experimental models in different studies.

This study was designed to characterize the NAc transcriptional and post-transcriptional landscape in mice with depression-like behaviors encouraged by chronic unpredictable mild stress (CUMS). We performed an unbiased, genome-wide bioinformatics analysis of mRNA transcriptome profiles and miRNA in the NAc from CUMS-induced depression mice with the use of high-throughput RNA sequencing. The miRNA–mRNA networks were created to probe and illustrate the complexity of stress-induced miRNA regulation. By using different associated analyses and comparisons, we found signal pathways in the NAc associated with stress-induced depression, in order to present the guidelines for understanding the molecular mechanisms of depression and exploring novel therapeutic strategies for this complex disorder.

Materials and methods

Animals

Male C57BL/6J mice (4 weeks, 16–18 g) were purchased from Beijing Vital River Laboratory Animal Technology Co., Ltd (Beijing, China). In this study, all procedures were carried out under the Guideline of National Institutes of Health, USA, for the Care and Use of Laboratory Animals and approved by Animal Care and Use Committee of Shandong University of Traditional Chinese Medicine Institutional (SDUTCM201805311223).

CUMS paradigm

The CUMS paradigms were performed with a minor modification according to the previous studies.^7,27 All experimental mice were acclimated to environment and sucrose intake for 1 week. We divided the experimental mice into two groups: control and CUMS. The 3–4 mice in the control groups were housed in one cage, and each mouse in the CUMS group was housed in a single cage. The control group was kept in undisturbed cages. A variate-stressor paradigm was used in the CUMS group.

The CUMS protocol consisted of 13 stressors including social isolation, food and water deprivation, empty bottles, soiled cage, restraint space, circadian disturbance, tilted cage, stroboscope, cold (4°C), rotating cage, wet cage, white noise, and space reduction (Table S1). This experiment was maintained for 4 weeks until the mice expressed depression-like behavior phenotype.

Behavioral assessments

Behavioral tests were conducted to assess whether the CUMS induced mice expressed depressed phenotype. The novelty-suppressed feeding test (NST) and sucrose preference test (SPT) were used to assess the anhedonia, and the tail suspension test (TST) and forced swimming test (FST) were used to estimate the behavioral despair. Behavioral assessment experiments were performed in the light cycle with a sound-proof behavioral facility.

TST and FST

TST and FST were carried out as previously reported.²⁸

SPT

SPT was conducted with a minor modification according to previous studies.^29,30

NST

Mice were observed for NST after 24-hour food deprivation (water was provided ad libitum). Details of the experiment have been previously reported.³⁰

RNA purification and sequencing

The protocol of RNA purification and sequencing has already been reported in our previously published study.⁷ In brief, if there were depression-like behaviors 24 hours after treatment (CUMS), both control and CUMS mice were anesthetized by using isoflurane and decapitated by guillotine. NAc samples were dissected and immediately placed on the ice-cold glass slides. These tissues were placed into the frozen vials that contained RNA and stored at -80°C for subsequent analyses. Total RNA was isolated from the colon with the RNAiso Plus Kit (TaKaRa, Dalian, China). After extracting RNA, the UV absorbance was used to measure RNA purity and quality. Absorbance at 260 nm was used to measure the amount of nucleic acid present in the sample. Absorbance at 280 nm were used to estimate the amount of protein in the sample. Absorbance at 230 nm was used to determine the amount of other contaminants that may be present in the samples, such as guanidine thiocyanate. The A260/A280 ratios ranging from 1.9 to 2.1 were acceptable ratios for purity in this experiment, and the requirement for A260/A230 ratios was >2.0. RNA concentration was calculated using the 260 nm reading. The samples with total RNA amount >10 µg, the concentration >200 ng/µL, were selected for the construction of transcriptome and small RNA libraries, respectively. The sequencings were done with the use of BGISEQ-500 platform (BGI, Shenzhen, China). The average reading length of two libraries was about 100 bp (pair-end) and 50 bp (single-end).

Bioinformatics

Bioinformatics analysis of transcriptome, miRNA expression profile, and integrated miRNA/mRNA network have been reported in our previous report.^7,30

Quantitative RT-PCR for the validations of miRNA and mRNA

To validate the RNA sequencing data, qRT-PCR was used to detect five miRNAs and five mRNAs which were significantly different between control and mice with depression-like behaviors (n=12/group), in which the tissues were included as those samples for sequencing (all gene primers are listed in Table S2). Briefly, total RNA was harvested from the NAc with the RNAiso Plus Kit (TaKaRa). Then, RNA was reverse-transcribed into cDNA with PrimeScript™ RT reagent Kit (TaKaRa). qRT-PCR was carried out using SYBR^® Premix Ex Taq™ kit (TaKaRa) on biosystems QuantStudio 7 Flex (Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer’s protocol. The relative expression level of mRNAs in the NAc was normalized to GAPDH. The relative miRNA expression in the NAc was normalized to U6. The 2^-ΔΔCt method was used to calculate the qRT-PCR results. Each sample of qRT-PCR was repeated in three replications.

Dual luciferase reporter assay

The 3′-untranslated region (UTR) sequence of Gad1 mRNA was amplified and fused into the NotI and XhoI sites of psiCHECK2. The site-directed mutation of the detected miR-144-3p targeting site of Gad1 mRNA 3′-UTR was constructed with the guideline of Site-Directed Mutagenesis Kit (Stratagene, La Jolla, CA, USA). The miRNA mimics and the negative control (primer sequences are listed in Table S3) were synthesized in Guangzhou Rui Bo Biological Technology Co., Ltd. HEK 293 T cells (24-well plates) were transfected with reporter vectors together with miR-144-3p mimics and inhibitor or the corresponding miR-NC. The renilla luciferase activity was measured and normalized to firefly luciferase activity at 48 hours post-transfection by using Dual-Glo^® Luciferase Assay System (E2920; Promega Corporation, Fitchburg, WI, USA). Each experiment was performed in triplicate.

Statistical analyses

The data of the behavior tests, gene analyses, and luciferase activity are presented as mean ± SEM. The association between miRNAs and prediction was investigated by Pearson’s correlation coefficients. Unpaired Student’s t-test was used to evaluate the statistical difference between the two groups. P<0.05 was considered statistically significant.

Results

CUMS induces depression-like behaviors

The mood states of CUMS-treated or control group mice were tested by NST, SPT, TST, and FST. Mice exposed to CUMS showed extensively improved immobility in FST (Figure 1A) and TST (Figure 1B) compared to control group, which demonstrated that those mice were in the status of helplessness and hopelessness. Stressed mice also displayed significantly reduced sucrose preference (Figure 1C) and increased latency to feed in the NST (Figure 1D), which showed that those mice were in the mood states of anhedonia and anxiety. Our data indicated that the CUMS paradigm induced depression-like behaviors.

Figure 1 CUMS leads the mice to express depression-like behaviors. Mice were subjected to the adaptation for a week, the CUMS for 4 weeks. (A) The values of immobile time in the FST were 170.68±1.28 seconds in CUMS-treated mice (n=12) and 141.60±1.8 seconds in controls (P<0.01). (B) The values of immobile time in the TST were 158.18±1.79 seconds in CUMS-treated mice and 124.94±0.86 seconds in controls (P<0.001). (C) The SPT values were 57.22%±1.9% in CUMS-treated mice (n=12) and 84.23%±0.92% in control mice (P<0.001). (D) The latency to eat food is 560.55±12.25 seconds in CUMS-treated mice and 344.58±7.22 in controls (P<0.001). The results are expressed as mean ± SEM. n=12 per group, ***P<0.001 compared with control.

Abbreviations: CUMS, chronic unpredictable mild stress; FST, forced swimming test; TST, tail suspension test; SPT, sucrose preference test; NST, novelty-suppressed feeding test; SEM, standard error of the mean.

Overall qualities of RNA-sequencing data set

In this study, we used RNA-sequencing and assessed miRNA expression and transcriptional profiles in NAc sample from mice with depression-like behaviors and control mice (n=4 per group). About 32,489,852–35,179,878 raw sequence reads were produced in the small RNA library. After trimming and filtering of reads with low quality and adaptor, a total of 29,759,815–32,525,311 clean small RNA reads were obtained (Table S4). The high qualities of small RNA sequencing and transcriptome data were further used for analysis. In parallel, RNAs from 45,123,294 to 45,175,652 raw sequence reads of about 100 bp were obtained from eight mRNA library Illumina sequencing. After filtering the reads that contained N with adaptor sequence and low quality, about 1.2 GB clean reads from each library were generated and mapped, which were equivalently about 75.82%–80.62% of the total reads from the University of California Santa Cruz mm10 (Table S5). In conclusion, these results from the high qualities of small RNA sequencing and transcriptome data can be used for further analysis.

miRNA expression profiles in NAc from CUMS-induced depression mice

The miRNA expression profile has been given in Table 1 with their expression change >1.5 fold and P-value <0.05. Seventeen miRNAs were significantly upregulated in the NAc of mice with CUMS-induced depression-like behaviors and 12 miRNAs were downregulated. These upregulated miRNAs included certain known miRNAs (mmu-miR-378d, mmu-miR-34b-5p, mmu-miR-34c-3p, mmu-miR-194-5p, mmu-miR-206-3p, mmu-miR-10a-5p, mmu-miR-144-3p, mmu-let-7a-1-3p, mmu-miR-34c-5p, mmu-miR-34b-3p, mmu-miR-21a-5p, mmu-miR-10b-5p, mmu-miR-31-5p, mmu-miR-219b-5p, and mmu-miR-338-3p), as well as two novel miRNAs (novel_mir42 and novel_mir39). Whereas, the downregulated miRNAs included mmu-miR-7a-2-3p, mmu-miR-124-5p, mmu-miR-351-5p, mmu-miR-127-5p, mmu-miR-671-5p, mmu-miR-335-3p, mmu-miR-377-3p, mmu-miR-574-5p, mmu-miR-200b-5p mmu-miR-486a-5p, mmu-miR-532-3p, and novel_mir29.

Table 1 miRNAs with quantitative change over 1.5-folds and their characteristics

miRNA id	Accession no	Control mean	CUMS mean	Fold change CUMS/control	P-value	Chromosomal location (mouse)	Seed sequence
mmu-miR-378d	MIMAT0025167	12.64	62.23	4.92	0.0072	chr10:126710282–126710391 [-]	5′-cuggccu-3′
mmu-miR-34b-5p	MIMAT0000382	5.89	15.09	2.56	0.0021	chr9:51103562–51103645 [-]	5′-ggcagug-3′
mmu-miR-34c-3p	MIMAT0004580	31.01	65.93	2.13	0.0049	chr9:51103034–51103110 [-]	5′-aucacua-3′
mmu-miR-194-5p	MIMAT0000224	106.99	227.42	2.13	0.0049	chr19:6264643–6264728 [+]	5′-guaacag-3′
mmu-miR-206-3p	MIMAT0000239	8.17	16.94	2.07	0.0105	chr1:20679010–20679082 [+]	5′-ggaaugu-3′
mmu-miR-10a-5p	MIMAT0000648	9.02	18.49	2.05	0.0003	chr11:96317165–96317274 [+]	5′-acccugu-3′
mmu-miR-144-3p	MIMAT0000156	9.83	19.94	2.03	0.0003	chr11:78073005–78073070 [+]	5′-acaguau-3′
mmu-let-7a-1-3p	MIMAT0004620	50.50	99.74	1.98	0.0169	chr13:48538179–48538272 [-]	5′-uauacaa-3′
mmu-miR-34c-5p	MIMAT0000381	25.35	44.39	1.75	0.0088	chr9:51103034–51103110 [-]	5′-ggcagug-3′
mmu-miR-34b-3p	MIMAT0004581	60.21	103.07	1.71	0.0436	chr9:51103562–51103645 [-]	5′-aucacua-3′
mmu-miR-21a-5p	MIMAT0000530	220.57	365.21	1.66	0.0011	chr11:86584067–86584158 [-]	5′-agcuuau-3′
mmu-miR-10b-5p	MIMAT0000208	8.05	13.02	1.62	0.0208	chr2:74726070–74726137 [+]	5′-acccugu-3′
mmu-miR-31-5p	MIMAT0000538	176.58	274.66	1.56	0.0269	chr4:88910557–88910662 [-]	5′-ggcaaga-3′
mmu-miR-219b-5p	MIMAT0029806	46.50	72.22	1.55	0.0001	chr2:29845647–29845711 [+]	5′-gaugucc-3′
mmu-miR-338-3p	MIMAT0000582	304.70	467.24	1.53	0.0262	chr11:120014765–120014862 [-]	5′-ccagcau-3′
novel_mir42	–	24	61	2.54	0.0230	chr11:109343222–109343303 [-]	5′-guggugc-3′
novel_mir39	–	26	65	2.50	0.0023	chr14:65530947–65530994 [-]	5′-gagguag-3′
mmu-miR-7a-2-3p	MIMAT0017070	73.66	49.38	0.67	0.0049	chr7:78888277–78888373 [+]	5′-aacaagu-3′
mmu-miR-124-5p	MIMAT0004527	39.79	25.94	0.65	0.0105	chr3:17795685–17795741 [-]	5′-guguuca-3′
mmu-miR-351-5p	MIMAT0000609	11.36	7.22	0.64	0.0001	chrX:53053255–53053353 [-]	5′-cccugag-3′
mmu-miR-127-5p	MIMAT0004530	23.02	14.08	0.61	0.0011	chr12:109592846–109592915 [+]	5′-ugaagcu-3′
mmu-miR-671-5p	MIMAT0003731	50.40	30.05	0.60	0.0003	chr5:24592114–24592211 [+]	5′-ggaagcc-3′
mmu-miR-335-3p	MIMAT0004704	54.34	32.28	0.59	0.0269	chr6:30741299–30741396 [+]	5′-uuuucau-3′
mmu-miR-377-3p	MIMAT0000741	53.54	30.86	0.58	0.0162	chr12:109740510–109740577 [+]	5′-ucacaca-3′
mmu-miR-574-5p	MIMAT0004893	47.17	27.33	0.58	0.0416	chr5:64970318–64970395 [+]	5′-gagugug-3′
mmu-miR-200b-5p	MIMAT0004545	49.92	27.02	0.54	0.0021	chr4:156055681–156055750 [-]	5′-aucuuac-3′
mmu-miR-486a-5p	MIMAT0003130	36.33	16.99	0.47	0.0049	chr8:23142555–23142682 [+]	5′-ccuguac-3′
mmu-miR-532-3p	MIMAT0004781	32.99	9.61	0.29	0.0088	chrX:7248402–7248497 [-]	5′-cucccac-3′
novel_mir29		198	108	0.55	0.0010	chr17:57117642–57117705 [-]	5′-cccugag-3′

Abbreviation: CUMS, chronic unpredictable mild stress.

In order to confirm the results of miRNA sequencing analysis, three upregulated miRNAs (mmu-miR-10a-5p, mmu-miR-144-3p, and mmu-miR-219b-5p) and two downregulated miRNAs (mmu-miR-124-5p and mmu-miR-127-5p) were selected for performing qRT-PCR. Consistent with high-throughput sequencing, these miR-NAs from mice with CUMS-induced depression-like behaviors are significantly changed in qRT-PCR compared with those from control mice (Figure 2A). These consistent results from the miRNA sequencing (Figure 2B) and qRT-PCR supported the validation of our study.

Figure 2 The validation of differentially expressed miRNAs in the NAc from mice with CUMS-induced depression-like behaviors and controls. Three upregulated miRNAs and two downregulated miRNAs were involved in different cellular functions and were selected for qRT-PCR analysis. (A) qRT-PCR was used to analyze the relative values of mmu-miR-10a-5p, mmu-miR-144-3p, mmu-miR-219b-5p, mmu-miR-124-5p, and mmu-miR-127-5p from mice with CUMS-induced depression-like behaviors and controls (n=12 per group), in which the samples were used from those tissues for high-throughput sequencing. (B) The relative level of mmu-miR-10a-5p, mmu-miR-144-3p, mmu-miR-219b-5p, mmu-miR-124-5p, and mmu-miR-127-5p from mice with CUMS-induced depression-like behaviors (n=4) and controls (n=4), which were analyzed by high-throughput miRNA sequencing. U6 was set as the internal control. The relative values for control mice were normalized to be 1. The data are expressed as mean ± SEM. *P<0.1, **P<0.01, ***P<0.001.

Abbreviations: CUMS, chronic unpredictable mild stress; NAc, nucleus accumbens; SEM, standard error of the mean.

Transcriptional profiles in the NAc from mice with CUMS-induced depression behaviors

mRNAs in the NAc from mice with depression-like behaviors and controls were quantified by sequencing total RNAs. We computed Reads Per Kilobase Million (RPKM) values for those genes. Genes with lower expression level (RPKM <5) were removed, thus 16,918 mRNAs were left for differential expression. The mRNA expression profile has been given in Table 1 in which alterations in their expressions were >1.5-fold in mice with depression-like behaviors and controls with P-value <0.05, which was a criterion to ensure differentially expressed genes. Twenty-nine mRNAs were significantly downregulated in the NAc from mice with CUMS-induced depression-like behaviors and 10 mRNAs were upregulated (Table 2). The decreased expression of mRNAs in the NAc of mice with CUMS-induced depression-like behaviors included Bdnf, Gad1, Atg13, Slc32a1 (VGAT), Gad2, Th, Mbp, and so on. Whereas, compared to control mice, Aspa, Atg2b, Hba-a1, Hba-a2, Hbb-b2, and Tnfsf12 mRNAs were significantly increased in the NAc of mice with CUMS-induced depression-like behaviors.

Table 2 mRNAs with quantitative change over 1.5-folds and their characteristics

Symbol	Gene ID	Length (bp)	Chromosomal map	Fold change CUMS/control	P-value	Description
Adra2a	11551	3,801	19 D2	0.62	0.0089	Adrenergic receptor, alpha 2a
Angpt2	11601	3,560	8 A1.3	0.52	0.0091	Angiopoietin 2
Atg13	51897	3,573	2 E1	0.63	0.0173	Autophagy related protein 13
Bdnf	12064	4,122	2 E3	0.57	0.0021	Brain-derived neurotrophic factor
Cdk5	12568	2,058	5 A3	0.54	0.0067	Cyclin-dependent kinase 5
Drd3	13490	1,455	16 B4	0.56	0.0012	Dopamine receptor D3
Fas	14102	1,486	19 C1	0.38	0.0052	Fas (TNF receptor superfamily member 6)
Fgl2	14190	3,769	5 A3	0.67	0.0021	Fibrinogen-like protein 2
Gh	14599	903	11 E1	0.02	0.0078	Growth hormone
Gad1	14415	40,249	2 C2	0.63	0.0038	Glutamic acid decarboxylase 1
Gad2	14417	71,550	2 A3	0.58	0.0029	Glutamic acid decarboxylase 2
Hspb1	15507	913	5 G2	0.44	0.0129	Heat shock protein 1
Icam1	15894	2,540	9 A3	0.46	0.0375	Intercellular adhesion molecule 1
Il1r1	16177	4,759	1 B	0.47	0.0129	Interleukin 1 receptor, type I
Map3k6	53608	4,333	4 D2	0.63	0.0044	Mitogen-activated protein kinase 6
Map3k8	26410	2,507	18 A1	0.62	0.0015	Mitogen-activated protein kinase 8
Mpzl2	14012	3,216	9 A5	0.65	0.0069	Myelin protein zero-like 2
Mgp	17313	599	14q11	0.51	0.0057	Matrix Gla protein
Mbp	17196	10,515	18 E3	0.36	0.0021	Myelin basic protein
Ngfr	18053	19,880	11 D	0.66	0.0025	Nerve growth factor receptor
Nr3c1	14815	6,329	18 B3	0.36	0.0148	Nuclear receptor subfamily 3, group C member 1
Ntn3	18209	4,911	17 A3	0.54	0.0051	Netrin 3
Pik3cd	18707	4,933	4 E2	0.62	0.0066	Phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit delta
Slc32a1	22348	4,090	2 H1	0.49	0.0029	GABA vesicular transporter
Socs3	12702	2,732	11 E2	0.64	0.0066	Suppressor of cytokine signaling 3
Th	21823	7,191	7 F5	0.65	0.0067	Tyrosine hydroxylase
Tlr7	170743	3,812	X F5	0.64	0.0074	Toll-like receptor 7
Tuba1c	22146	1,977	15 F1	0.61	0.0021	Tubulin, alpha 1C
Tyrp1	22178	2,741.3	4 C3	0.48	0.0011	Tyrosinase-related protein 1
Aspa	11484	1,537	11 B4	1.80	0.0003	Aspartoacylase
Atg2b	76559	10,173	12 E	1.87	0.0002	Autophagy-related protein 2B
Hba-a1	15122	559	11 A4	2.06	0.0039	Hemoglobin alpha, adult chain 1
Hba-a2	110257	563	11 A4	2.24	0.0006	Hemoglobin alpha, adult chain 2
Hbb-b2	15130	630	7 E3	2.52	0.0002	Hemoglobin, beta adult minor chain
Irf7	54123	1,809.78	7 F5	1.79	0.0154	Interferon regulatory factor 7
Mknk2	17347	3,448	10 C1	1.89	0.0015	MAP kinase-interacting serine/threonine kinase 2
Rasgrp2	19395	2,193	19 A	2.24	0.0025	RAS, guanyl releasing protein 2
Sema4g	26456	4,291	19 C3	1.67	0.0032	Sema domain, immunoglobulin domain (Ig), transmembrane domain (TM), and short cytoplasmic domain, (semaphorin) 4G
Tnfsf12	21944	1,694	11 B3	1.59	0.0087	Tumor necrosis factor (ligand) superfamily, member 12

Abbreviations: TNF, tumor necrosis factor; GABA, γ-aminobutyric acid.

Based on the bioinformatics analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) database, these differentially expressed genes (DEGs) were significantly enriched in the following signaling pathways and processes, such as the MAPK signaling pathway, GABAergic synapse, dopaminergic synapse, neuroactive ligand–receptor interaction, glutamate metabolism, cytokine–cytokine receptor interaction, axon guidance, regulation of autophagy, apoptosis, and so on (Table 3). Thus, the GABAergic synapses, dopaminergic synapses, axon guidance, neurotrophin signaling pathway, neurotransmitter synthesis, autophagy-associated apoptotic pathway, and neural-immune process in the NAc were associated with the pathogenesis of depression. It was noteworthy that the unregulated genes were clustered in hemoglobin genes (Hba-a1, Hba-a2, and Hbb-b2).

Table 3 Signaling pathways identified by KEGG function analysis based on DEG data

Pathway	DEGs with pathway annotation (39)	All genes with pathway annotation (32,152)	Contributed genes	P-valuea	Pathway ID
African trypanosomiasis	5 (12.82%)	65 (0.2%)	Fas, Icam1, Hba-a1, Hba-a2, Hbb-b2	0.00657	ko05143
MAPK signaling pathway	8 (20.51%)	684 (2.13%)	Bdnf, Fas, Hspb1, Il1r1, Map3k6, Map3k8, Mknk2, Rasgrp2	0.00033	ko04010
Alanine, aspartate, and glutamate metabolism	3 (7.69%)	68 (0.21%)	Gad1, Gad2, Aspa	0.00069	ko00250
TNF signaling pathway	5 (12.82%)	280 (0.87%)	Fas, Icam1, Map3k8, Pik3cd, Socs3	0.00072	ko04668
Toll-like receptor signaling pathway	4 (10.25%)	233 (0.72%)	Map3k8, Pik3cd, Tlr7, Irf7	0.00289	ko04620
HIF-1 signaling pathway	4 (10.25%)	234 (0.73%)	Angpt2, Fgl2, Pik3cd, Mknk2	0.00294	ko04066
Cocaine addiction	3 (7.69%)	126 (0.39%)	Bdnf, Cdk5, Th	0.00401	ko05030
Tyrosine metabolism	2 (5.12%)	65 (0.2%)	Th, Tyrp1	0.01177	ko00350
GABAergic synapse	3 (7.69%)	224 (0.7%)	Gad1, Gad2, Slc32a1 (VGAT)	0.01912	ko04727
Ras signaling pathway	5 (12.82%)	664 (2.07%)	Angpt2	0.02626	ko04014
Natural killer cell mediated cytotoxicity	3 (7.69%)	263 (0.82%)	Fas, Icam1, Pik3cd	0.02894	ko04650
Cytokine-cytokine receptor interaction	4 (10.25%)	508 (1.58%)	Fas, Il1r1, Ngfr, Tnfsf12	0.03957	ko04060
Axon guidance	4 (10.25%)	559 (1.74%)	Cdk5, Ntn3, Pik3cd, Sema4g	0.04284	ko04360
Regulation of autophagy	3 (7.69%)	364 (1.13%)	Atg13, Pik3cd	0.02451	ko04140
Jak-STAT signaling pathway	3 (7.69%)	374 (1.16%)	Gh, Pik3cd, Socs3	0.0335	ko04630
Apoptosis	3 (7.69%)	384 (1.19%)	Fas, Pik3cd, Tuba1c	0.03251	ko0421
Neurotrophin signaling pathway	3 (7.69%)	398 (1.24%)	Bdnf, Ngfr, Pik3cd	0.0312	ko04722
Neuroactive ligand-receptor interaction	4 (10.25%)	649 (2.02%)	Adra2a, Drd3, Gh, Nr3c1	0.01482	ko04080
Dopaminergic synapse	2 (5.12%)	401 (1.25%)	Drd3, Th	0.01259	ko04728

Note:

a P-values from the hypergeometric tests were adjusted by Benjamini-Hochberg method.

Abbreviations: DEG, differentially expressed gene; KEGG, Kyoto Encyclopedia of Genes and Genomes; TNF, tumor necrosis factor.

In order to verify the differentially expressed genes data, five significantly changed mRNAs involved in different cellular functions were selected for performing qRT-PCR from the NAc tissues that had been used for the mRNA sequencing. The level of Hbb-b2 (P<0.001) gene expression was increased, and the expressions of Bdnf, Slc32a1(VGAT), Mbp, and Gad1 were decreased in mice with CUMS-induced depression-like behaviors compared with that in the control mice (all P<0.001, Figure 3A). The consistent results obtained by gene sequencing (Figure 3B) and qRT-PCR supported our research validation.

Figure 3 The validation of differentially expressed mRNAs in the NAc from mice with CUMS-induced depression-like behaviors and controls.

Notes: Five significantly changed mRNAs involved in different cellular functions were selected for performing qRT-PCR. (A) qRT-PCR was used to analyze the relative values of Hbb-b2, Bdnf, Slc32a1(VGAT), Mbp, and Gad1 from mice with CUMS-induced depression-like behaviors and controls (n=12 per group). The samples contained those tissues for high-throughput sequencing. (B) The relative level of Hbb-b2, Bdnf, Slc32a1(VGAT), Mbp, and Gad1 gene expression from mice with CUMS-induced depression-like behaviors (n=4) and controls (n=4), which were analyzed by high-throughput mRNA sequencing. GAPDH was set as the internal control. The relative values for control mice were normalized to be one. The data are expressed as mean ± SEM. **P<0.01, ***P<0.001.
Abbreviations: CUMS, chronic unpredictable mild stress; NAc, nucleus accumbens; SEM, standard error of the mean.

Integrated miRNA/mRNAregulatory networks

The level of mRNA expression in the cells was affected by miRNAs, through which the bindings of miRNAs with their dicers degrade mRNAs and weaken their translations. If the downregulated mRNAs in the NAc from mice with depression-like behaviors were caused by miRNAs, their corresponding miRNAs will be upregulated. To test this hypothesis and validate our data about mRNA changes, we performed a series of bioinformatics analyses.

The significantly changed miRNAs predicted target mRNAs well match the measures by mRNA sequencing based on the three databases (Targetscan, RNA22, and miRDB) about the complex interactions between miRNAs and mRNAs. Table 4 shows the alternated miRNAs and their predicted-target mRNAs. For instance, the unregulated mmu-miR-144-3p predicted the downregulated Gad1 and Nr3c1 as the target mRNAs. The altered mRNAs and their corresponding miRNAs are shown in Table 5. Among those predicted altered mRNAs, Bdnf and mRNAs could be predicted as the four changed miRNAs (mmu-let-7a-1-3p, mmu-miR-10a-5p, mmu-miR-206-3p, and mmu-miR-10b-5p). From Tables 1 and 3–5, we can find that the altered mRNAs and miRNAs were simultaneously and reversely changed. The consistent results by associatively sequencing mRNAs and miRNAs validated our analyses and strengthen our conclusion.

Table 4 The changed miRNAs predict target mRNAs

miRNAs	The predicted target mRNAs that match DEGs in transcriptomea
mmu-let-7a-1-3p ↑	Bdnf ↓, Gad1 ↓
mmu-miR-10a-5p ↑	Bdnf ↓
mmu-miR-378d ↑	Zic1 ↓
mmu-miR-34b-5p ↑	Ntn3 ↓
mmu-miR-206-3p ↑	Atg13 ↓, Bdnf ↓
mmu-miR-144-3p ↑	Gad1 ↓, Nr3c1 ↓
mmu-miR-34c-5p ↑	Ntn3 ↓
mmu-miR-34b-3p ↑	Angpt2 ↓
mmu-miR-10b-5p ↑	Bdnf ↓
mmu-miR-31-5p ↑	Gad2 ↓
mmu-miR-219b-5p ↑	Th ↓, Tyrp1 ↓
mmu-miR-7a-2-3p ↓	Atg2b ↑
mmu-miR-124-5p ↓	Aspa ↑, Irf7 ↑
mmu-miR-351-5p ↓	Mknk2 ↑
mmu-miR-127-5p ↓	Mknk2 ↑
mmu-miR-200b-5p ↓	Atg2b ↑

Notes:

a The target mRNAs should be predicted by Targetscan, RNA22, and miRDB, and then overlapped to DEGs in transcriptome. ↑ indicates up-regulation in the tissue of NAc from depression-like mice vs control mice, whereas ↓ represents down-regulation.

Abbreviation: DEGs, differentially expressed genes.

miRNA-144-3p can bind to the GAD1 mRNApredicted region

To confirm the silico prediction, we selected miRNA-144-3p to examine whether they targeted Gad1 by qRT-PCR and dual luciferase reporter assay. Two areas of Gad1 3′-UTRs (769–775 and 902–908) were predicted as the miRNA-144-3p binding site according to the three databases (Figure 4A). The results of qRT-PCR analysis showed that there were negative correlations between miRNA-144-3p and Gad1 mRNA expression in NAc tissue from mice with depression-like behaviors and controls (r=-0.946; P<0.001; Figure 4B). Furthermore, in dual luciferase report assay, we constructed luciferase reporter plasmids, which contained the two wild-type or corresponding mutant 3′-UTRs of the predicted binding sites. These reporter constructs were transfected into HEK293T cells. After miRNA-144-3p mimics or its negative controls were incubated, the relative activities of luciferase reporter for the 3′-UTR 769–775 and 902–908 of Gad1 mRNA were significantly lowered by the mimics of miRNA-144-3p, which were reversed by mutating the binding sites of miRNA-144-3p. Interestingly, this phenomenon can also be eliminated by mutating both the binding sites of Gad1, which appeared to have the synergistic effect (Figure 4C). These results suggested that miR-144-3p might act as negative regulators of translation through binding to the two sites of Gad1 3′-UTR, which also supported our analyses for the prediction of miRNA target genes.

Figure 4 miRNA-144-3p can bind to the Gad1 mRNA predicted region.

Notes: (A) miRNA-144-3p targeted to Gad1 was predicted by three miRNA target prediction databases (Targetscan, RNA22, and miRDB). The sequences of these seeds referring to the nucleotides in miRNA positions 2–8 are shown in red. Watson–Crick matches in the seed sequence are shown in blue. (B) Correlation between miRNA-144-3p and its prediction target Gad1 mRNA expression in NAc tissue from mice with depression-like behaviors and controls (r=-0.946; P<0.001). (C) Luciferase reporter assay was performed by the co-transfection of luciferase reporter containing the wild or mutant 3′-UTR of GAD1 mRNA with miR-144-3p mimic or their NC into HEK293T cells. Luciferase activity was determined 48 hours after co-transfection. The data are expressed as mean ± SEM. **P<0.01, ***P<0.001.
Abbreviations: NC, negative control; SEM, standard error of the mean.

Discussion

MDD is a recurrent, devastating mental disorder, which affects >350 million people worldwide, and exerts substantial public health and financial costs to society. Thus, there is a significant need to discover innovative therapeutics to treat depression efficiently. Growing evidence indicates that dysfunctions in brain reward circuitry are major causes of mood disorders including depression. In the present study, we performed transcriptome analysis in CUMS model and reported differentially expressed genes and miRNAs in the NAc, a key brain region of reward circuitry.

The differentially expressed mRNAs were significantly enriched in GABAergic synapses, dopaminergic synapses, axon guidance, neurotrophin signaling pathway, neurotransmitter synthesis, autophagy-associated apoptotic pathway, and neural-immune process. It is noteworthy that the upregulated mRNAs were clustered in hemoglobin genes. In addition, target miRNAs and mRNAs were upregulated. The inline results from sequencing miRNAs and mRNAs (Tables 1 and 3) as well as qRT-PCR (Figures 2 and 3) fortify our results. Moreover, the consistent results of qRT-PCR and dual luciferase reporter assays also strengthen the bio-informatics analysis of miRNA target prediction (Figure 4). To the best of our knowledge, here for the first time, we reveal the molecular mechanism of the NAc pathophysiology in CUMS-induced depression.

Various research reports have highlighted the significance of NAc in MDD.^12,31–35 In humans, deep brain stimulation of the NAc has been shown to alleviate anhedonic symptoms of depression.³⁶ However, the NAc is a heterogeneous structure comprising different cell types, various dopamine receptors, and other neuromodulatory signaling. These characteristics made it difficult to elucidate the precise role of the NAc circuitry in depressive behaviors. Increasing evidence supports the reduction of synaptic transmission and GABAergic neurons in the NAc associated with depression.^37–39 Vulnerable mice exhibit an elevated stubby spine NAc that is associated with high occurrence of miniature excitatory post-synaptic current.⁴⁰ Electrophysiological record of NAc excitatory diffusion followed stress expose were unbalanced reduction.⁴¹ Also, in our previous study, we found that GABAergic neurons were attributed as inhibitory synapse outputs, excitatory, and excitability synapse reception in CUMS-induced depression mice.⁸ Such indications are completely in-line with our observations (Tables 1–5). As GABAergic neurons regulated the coordination and encoding ability of excitatory neurons,^42,43 the attenuation of GABAergic neurons led to the dysfunction of neuronal networks in the NAc. Furthermore, our analyses indicated the functional defect in dopaminergic synapses, axon guidance, neurotrophin signaling pathway, and neuroactive ligand–receptor interaction. As was known to us, depression was characterized as low self-esteem and anhedonia, which might have occurred due to dysfunction of the synapses and neurons in the reward circuits, such as the synaptic transmission deficit and the degeneration of the axons in the NAc and medial PFC.⁴⁴ Taken together, these results provided insight into mechanisms underlying the physiological absence of the reward circuits in psychiatric disorders for example, bipolar disorder, MDD, attention schizophrenia and post-traumatic stress disorder.

Earlier studies indicated downregulation in the genes of MAPK signaling pathway, neurotransmitter synthesis and myelin-associated protein in rodent animal model of depression-like and depressive patients with the use of qRT-PCR, microarray or in situ hybridization.^45–49 These indications were also consistent with the current study in CUMS-induced depression-like mice. In addition, our analysis further implies autophagy-associated apoptotic pathway and neural-immune process, which may bring new insight into molecular mechanism essential depression. These latest findings might be due to the use of advanced mRNA sequencing method.

In this study, the leading observations were that the upregulated mRNAs were clustered in hemoglobin genes. Hemoglobin expression is high in erythrocytes. It has a key role in oxidative stress, response to injury, and neuronal respiration.^50,51 Cells can use neuronal hemoglobin to produce hemoglobin derived peptides action on opioid and cannabin receptors.⁵² Chronic peripheral inflammation was contributed to increase the frontal expression of hemoglobin genes in ischemic and intracerebral hemorrhage patients.^53,54 Upregulation of hemoglobin genes (Hbb-b1, Hba-a1, and Hba-a2) has also been found in the medial PFC of mice with CUMS-induced depression-like behaviors in our previous study.⁷ Human HbF and HbA2 were associated with disease progression in bipolar disorder, which had a shielding function of HbA2 against postpartum episodes.⁵⁵ Thus, the elevation in hemoglobin genes was found in diverse pathological states. Consistent results of related studies reveal hemoglobin genes as the potential markers of chronic mild stress-induced depression.

miRNAs are a group of single-stranded small non-coding RNA molecules with 19–22 nucleotides, which could post-transcriptionally degrade mRNA or inhibit the translation of mRNA via binding the 3′-UTR section of mRNAs and furthermore influence the expression of target genes.⁵⁶ In some cases, miRNAs might trigger translation or even control transcription level by binding to precise gene promoters. The role of miRNAs in depression is mainly explained from preclinical studies which showed the role of miRNAs in neurogenesis, neurodevelopment, and synaptic plasticity.^57,58 Accumulated facts imply that stress is a significant risk factor in depression and can differentially regulate miRNA-mediated expression of mRNA transcripts in the brain of rodent.^59–61 Additionally, mRNAs down-regulation in the mice NAc are rooted by miRNAs, and their related miRNAs will also be upregulated. To validate the hypothesis and strengthen our findings about mRNA changes, we investigated the miRNA expression in NAc by their sequencings from mice with depression-like behaviors.

In our current study, we found that some known miRNAs were significantly changed in mice with CUMS-induced depression-like behaviors (Table 1), which degraded mRNAs listed in Table 5. In other words, the analysis from miRNA sequencing was consistent with the analysis from mRNA sequencing. Furthermore, these dysregulation miRNAs were reported in previous studies by qRT-PCR, microarray, or situ hybridization in rodent depression animal model or depressive patients.^62–68 We next performed a series of bioinformatics analyses for integrated miRNA/mRNA regulatory networks (Tables 4 and 5). For instance, the unregulated mmu-miR-144-3p predicted the downregulated Gad1 and Nr3c1 as the target mRNAs. Whereas, the downregulated mmu-miR-124-5p predicted the unregulated Aspa and Irf7 as the target mRNAs. Among those predicted altered mRNAs, Bdnf mRNA could be predicted as the four changed miRNAs (mmu-let-7a-1-3p, mmu-miR-10a-5p, mmu-miR-206-3p, and mmu-miR-10b-5p).

Table 5 The changed mRNAs are regulated by miRNAs

mRNAs	miRNAs that are predicted to mRNAs by Targetscan, RNA22 and miRDB
Bdnf ↓	mmu-let-7a-1-3p ↑, mmu-miR-10a-5p ↑, mmu-miR-206-3p ↑, mmu-miR-10b-5p ↑
Gad1 ↓	mmu-miR-144-3p ↑, mmu-let-7a-1-3p ↑
Zic1 ↓	mmu-miR-378d ↑
Ntn3 ↓	mmu-miR-34b-5p ↑, mmu-miR-34c-5p ↑
Atg13 ↓	mmu-miR-206-3p ↑
Nr3c1 ↓	mmu-miR-144-3p ↑
Angpt2 ↓	mmu-miR-34b-3p ↑
Gad2 ↓	mmu-miR-31-5p ↑
Th ↓	mmu-miR-219b-5p
Tyrp1 ↓	mmu-miR-219b-5p ↑
Atg2b ↑	mmu-miR-7a-2-3p ↓, mmu-miR-200b-5p ↓
Aspa ↑	mmu-miR-124-5p ↓
Irf7 ↑	mmu-miR-124-5p ↓
Mknk2 ↑	mmu-miR-351-5p ↓, mmu-miR-127-5p ↓

Note: ↑ indicates up-regulation in the tissue of NAc from depression-like mice versus control mice, whereas ↓ represents down-regulation.

To verify the silico prediction, we selected miRNA-144-3p to test whether they targeted Gad1 by dual luciferase reporter assay. Our results demonstrated that miR-144-3p acted as negative regulators of translation through attaching with the two sites of Gad1 3′-UTR. miR-144-3p role still remained unclear in depressive disorders. But, some biological mechanisms can explain our judgments. miR-144-3p has extensive expression nature which is enriched in human brain, additionally, also enriched in malignant hematopoietic and normal cells and tissues.⁶⁹ miR-144-3p is very conserved and has various envisaged targets in both mice and humans. Numerous investigations have revealed that miR-144-3p was implicated against response to stress, aging diseases, and mood stabilizer treatment.^7,70,71 miR-144-3p targeting signaling pathways include Nrf2, Wnt/β-catenin, and MAKP pathways.^72–74 In addition, miR-144-3p may suppress ataxin 1 (ATXN1) expression in human cells, and interestingly, the Genetic Association Database illustrated that ATXN1 is linked with psychological disorders.⁷⁵ In terms of response to stress, miR-144-3p level is notably elevated in depressed patients compared with healthy young adults. These studies support that miR-144-3p-mediated target gene expressions were involved in the processes of the pathogenesis of depression.

The quantitative alterations of miRNAs in some reports without the alterations of their targeted mRNAs in others may be caused by the following reasons. The changes in miRNAs may not reach the threshold to regulate their targeted mRNAs.⁷⁶ circRNAs act as miRNA sponges and positive regulators of miRNA-targeted genes.⁷⁷ On the other hand, the quantitative changes of mRNAs without the alteration of their corresponding miRNAs imply that the altered mRNAs are likely regulated by other epigenetic mechanisms, such as DNA methylation and repressive histone modification in the promoters.^78–80 Therefore, the associated analysis of miRNA/mRNAs in the brain regions with depression-related dysfunction may validate the data and strengthen the conclusion.

Conclusion

In summary, we performed a series of bioinformatics analyses for RNA sequencing data in the NAc tissue from the CUMS-induced depression mice. The worsening of dopaminergic synapses, GABAergic synapses, and neurotransmitter syntheses and autophagy-associated apoptotic pathway are associated with the molecular pathological mechanism of CUMS-induced depression. Our analyses support the string from stress to neuron atrophy through miRNA/mRNA regulatory network and provide the guidelines for developing novel therapeutic strategies for this complex disorder.

Author contributions

All authors contributed toward data analysis, drafting and critically revising the paper, approved the final version to be published, and agree to be accountable for all aspects of the work.

Acknowledgments

This work was supported by the Shandong Co-Innovation Center of Classic TCM Formula, Scientific Innovation Team of Shandong University of Traditional Chinese Medicine, the National Natural Science Foundation of China (81673852 and 81603419), Shandong Province University Scientific Research Project (J18KZ014), Shandong Medical Health Technology Development Plan Project (2018WS203), and Student Research Training Program of Shandong University of TCM (2018090 and 2018093).

Supplementary materials

Table S1 Chronic unpredictable mild stress procedure

Stressor	Mon	Tue	Wed	Thu	Fri	Sat	Sun
Social isolation	09:00→09:00	09:00→09:00	09:00→09:00	09:00→09:00	09:00→09:00	09:00→09:00	09:00→09:00
Food and water deprivation	09:00→09:00
Exposure to empty bottles		9:00→11:00
Soiled cage		11:00→11:00			9:00→21:00
Restraint			11:00→12:00
Light/dark succession every 2 h			12:00→22:00
45° cage tilt				9:00→21:00
Stroboscope				21:00→9:00
Cold (4° for 1 h)					21:00→22:00
Cage rotation	11:00→12:00					9:00→10:00
Wet cage						10:00→22:00
White noise							9:00→21:00
Space reduction							21:00→9:00

Table S2 qRT-PCR prime information

Gene ID	Symbol	Prime sequence	Length (bp)	Tm (°C)
15130	Hbb-b2	Forward 5′-CTGATTCTGTTGTGTTGACTTG-3′	188	60
		Reverse 5′-AGGTCTCCAAAGCTATCAAAGT-3′
12064	Bdnf	Forward 5′-GAGACAAGAACACAGGAGGAAA-3′	176	60
		Reverse 5′-GACTAGGGAAATGGGCTTAACA-3′
22348	Slc32a1	Forward 5′-TGGTCATCGCTTACTGTCTC-3′	157	60
		Reverse 5′-TGCTGCATGTTGCCTTCG-3′
17196	Mbp	Forward 5′-ACCATCCAAGAAGACCCCA-3′	192	60
		Reverse 5′-ACCCCTGTCACCGCTAAAG-3′
14415	Gad1	Forward 5′-GGGCTATGTTCCCCTTTATGT-3′	184	60
		Reverse 5′-CCTTTCTATGCCGCTGAGT-3′
		Reverse 5′-TTTGATGTCACGCACGATTT-3′
11461	GAPDH	Forward 5′-CTACGAGGGCTATGCTCTCC-3′	145	60
		Reverse 5′- TTTGATGTCACGCACGATTT-3′

Table S3 3′-Untranslated region (UTR) and site-directed mutation prime sequence Gad1

Gene ID	Symbol	Accession	Prime sequence
14415	Gad 67	AF326547.1	Forward 5′-CCGCTCGAGTGTCTATCTTTGGGCAGGGG-3′
			Reverse 5′-ATAAGAATGCGGCCGCAACACTTGTGGGACTGGTCA-3′
Gad1 (769–775)			Forward 5′-TGGAGTAAGTTCTGGTATAAAATTATGGTATTTTCGT-3′
			Reverse 5′-ACGAAAATACCATAATTTTATACCAGAACTTACTCCA-3′
Gad1 (902–908)			Forward 5′-TGGAAAAGAATTCTGTATAATTACATAGAGTCATGTT-3′
			Reverse 5′-AACATGACTCTATGTAATTATACAGAATTCTTTTCCA-3′

Table S4 Filtering small RNA library raw data and quality control

Sample name	Sequence type	Raw tag count	Clean tag count	%
Control-1	SE50	32,489,852	29,759,815	91.6
Control-2	SE50	33,436,523	31,008,054	92.74
Control-3	SE50	34,718,657	32,042,912	92.29
Control-4	SE50	35,179,878	32,525,311	92.45
CUMS-1	SE50	34,771,677	32,514,665	93.51
CUMS-2	SE50	33,954,639	31,391,222	92.45
CUMS-3	SE50	32,524,414	30,252,851	93.02
CUMS-4	SE50	34,232,788	32,024,156	93.55

Table S5 Filtering transcriptome raw data and alignment

Sample name	Clean reads	Genome map rate (%)	Gene map rate (%)	Expressed gene	Expressed transcripts	Expressed exon	Extend gene
Control-1	126,728,056	78.83	72.17	16,651	20,371	189,111	2,631
Control-2	127,291,954	75.82	71.41	16,534	20,128	188,520	2,589
Control-3	127,371,504	78.95	72.16	16,681	20,321	189,181	2,647
Control-4	127,218,882	74.49	72.41	16,519	20,111	188,326	2,569
CUMS-1	127,219,734	79.62	74.53	16,508	20,182	186,646	2,581
CUMS-2	126,569,138	79.58	72.58	16,636	20,315	188,256	2,597
CUMS-3	126,395,876	80.62	76.26	16,513	20,126	186,636	2,578
CUMS-4	126,021,758	78.12	74.51	16,649	20,352	188,217	2,557

Abbreviation: CUMS, chronic unpredictable mild stress.

Disclosure

The authors report no conflicts of interest in this work.