Abstract Background Major depression is a prevalent mood disorder. Chronic stress is presumably main etiology that leads to the neuron and synapse atrophies in the limbic system. However, the intermediate molecules from stresses to neuronal atrophy remain elusive, which we have studied in the medial prefrontal cortices from depression mice. Methods and Results The mice were treated by the chronic unpredictable mild stress (CUMS) until they expressed depression-like behaviors confirmed by the tests of sucrose preference, forced swimming and Y-maze. High-throughput sequencings of microRNA and mRNA in the medial prefrontal cortices were performed in CUMS-induced depression mice versus control mice to demonstrate the molecular profiles of major depression. In the medial prefrontal cortices of depression-like mice, the levels of mRNAs that translated the proteins for the GABAergic synapses, dopaminergic synapses, myelination, synaptic vesicle cycle and neuronal growth were downregulated. miRNAs of regulating these mRNAs are upregulated. Conclusion The deteriorations of GABAergic and dopaminergic synapses as well as axonal growth are associated with CUMS-induced depression. Introduction Major depressive disorder is featured as anhedonia and low self-esteem. The chronic stress to the genetically susceptible subjects leads to the deficits of monoamine, brain-derived neurotrophic factor and hypothalamus-pituitary-adrenal pathway [[32]1–[33]6], which induce neuronal atrophy in brain reward circuits, such as prefrontal cortex and nucleus accumbens, in depression patients and depression-like animals [[34]7–[35]13]. Genetic analyses by blood and dermal cells indicate that the deficits of numerous genes are related to genetic susceptibility for major depression [[36]14–[37]18]. However, these data may not reflect how stress environments influence these molecules in the brain that mediate pathological changes in major depression, such that gene analyses in depression-related brain regions are needed [[38]19–[39]22]. The analyses of gene expression in the limbic system from depression subjects indicate that many genes are involved in major depression, in which some of them differ from the analyses in the peripherals. The changes of these genes may lead to the depressed cell metabolisms and synapse functions [[40]19–[41]23]. The change of gene expression in the brain is presumably caused by a process that environmental factors make epigenetic dysregulation and in turn lead to major depression [[42]24–[43]26]. This assumption is granted by analyzing microRNA (miRNA) in certain brain areas from the postmortem tissues of major depression patients [[44]27–[45]31] and from the brain tissues of depression-like mice [[46]32, [47]33]. The reaction chain from stress to miRNA dysregulation, messenger RNA (mRNA)/protein expression alternations and neuronal atrophy is presumably related to major depressive disorders. However, the molecules identified for major depression appear varied in these reports. The inconsistence may be due to the analyses applied in these studies from the different experiment models, brain areas, animal species and/or the approaches by analyzing mRNA or miRNA. For instance, the changes of miRNAs in the reports without the alternations of their targeted mRNAs in other reports imply that these miRNAs may not reach the thresholds to regulate their mRNAs. The changes of mRNAs without the alternations of their correspondent miRNAs imply that these changed mRNAs are regulated by other mechanisms. Thus, the associative assessments of miRNAs and mRNAs [[48]14] from the brain regions that are involved in depression-related dysfunctions [[49]34] may examine the data validation and strengthen the conclusion about the chain from stress to neuron atrophy via miRNA and mRNA/protein expression dysregulation. We have studied molecule profiles related to major depression in depression-like mice induced by chronic unpredictable mild stress (CUMS), since the animals were thought to be the desirable system for gene analyses in the brain to reveal the molecules relevant to this disorder [[50]6, [51]35] and the neuropathological changes were clearly detected in medial prefrontal cortices from CUMS-induced depression mice [[52]34]. The gene expression was analyzed by sequencing miRNA and mRNA in the medial prefrontal cortices from CUMS-induced depression mice versus control. By the associated analyses and comparisons, we expect to figure out the signaling pathways in the medial prefrontal cortex related to stress-induced depression, in order to provide the guidelines for addressing the molecular mechanisms of major depressive disorder and for developing its therapeutic strategies. Methods and Materials Ethic issue All experiments were conducted in accordance with the guidelines and regulations by Administration Office of Laboratory Animals at Beijing China. All experimental protocols were approved by Institutional Animal Care Unit Committee in Administration Office of Laboratory Animals at Beijing China ([53]B10831). In terms of the living conditions for the mice under the normal life and control group, they are housed in the cages (32×16×16 cm) with free access water and food pellets under the circadian of 12 hours in the light (7:00 am to 7:00 pm) and the rest of 12 hours in the dark. Ambient temperature and relative humidity are maintained at 22±2°C and 55±5%, respectively. The standards are maintained in the specific pathogen free (SPF). The procedures of chronic unpredictable mild stress (CUMS) C57BL/6J juvenile mice were used to analyze the molecular profiles associated with major depression. The male mice were used starting at postnatal days 21 for making depression model. In week one for their adaptation to the experiments, the body weight, locomotion, sucrose preference and Y-maze tests of these mice were measured to select the healthy ones for our study. The mice of showing consistent values in these measurements were separated into two groups, CUMS and control, in order to reduce the variations among these mice. The control mice lived without the following stresses. Based on depression risk factors, such as weaknesses in cognitive function, emotional regulation, social interaction skill, circadian and stress responses [[54]36], we used chronic stress to produce depression-like mice in the following principle. The mice lived in stressful environment, made efforts to challenge these conditions and experienced defeat outcomes, which drove them to feel cognitive and emotional inabilities and in turn to be anhedonia and low self-esteem. The procedures for the CUMS mice include their adaptation, the CUMS and the behavioral tests [[55]34]. The stressful environments included social isolation, tilted cage, empty cage, damp sawdust cage, restraint space, white noise, strobe light and circadian disturbance [[56]1, [57]5, [58]37–[59]39]. Except for the social isolation, these conditions were randomly selected to treat the mice in the manners of their separations or combinations every day. These treatments were applied about 1~14 hours in durations and 1~12 hours in intervals [[60]34]. The duration and interval were unpredictable to the mice. The CUMS was sustained for three weeks until some of these mice expressed anhedonia and low self-esteem. Behavior assays for major depressive disorders Whether the CUMS-treated mice in three weeks fell into anhedonia and low self-esteem was tested in day 29~31. The sucrose preference test (SPT) and Y-maze test (YMT) were used to assess the anhedonia, and the forced swimming test (FST) was used to estimate their self-esteem [[61]9, [62]40–[63]43]. The SPT was conducted by 1% sucrose water versus water for four hours. The SPT value was presented as a ratio of the ingested sucrose water to the ingested sucrose water plus water. The YMT was performed by monitoring mouse staying in a special arm and other two arms for 2 minutes. The end of this special arm included a female mouse (named as M-arm). M-arm stay time was presented by a ratio of stay time in M-arm to that in three arms. The FST was done by recording immobile time in the water cylinder (10 centimeters in diameters and 19 centimeters in water depth at 25±1°C) for 6 minutes. To quantify the FST, immobile time and latency (a period to mouse immobility in the first time) were presented. The SPT, YMT and FST was done before and after the CUMS. Before each SPT, the mice in the CUMS and control were deprived from food and water for 3 hours to drive their intension of drinking water. In the YMT, these 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 additional stresses, same circadian circle for all mice and their adaptation in the test environment [[64]34]. An expression of depression-like behaviors was accepted if the mice in the CUMS group showed the decreases in sucrose preference ratio and M-maze stay time as well as the increase in immobile time significantly, compared to these values during their self-control period (the first week for adaptation) and in the control group of the mice. The mice with significant changes in all of three tests were defined as CUMS-induced depression-like mice or depression-like mice. The CUMS-treated mice in 3 weeks met this criterion about 30%, implying their vulnerability to the stresses. These mice were used as depression-like mice to study the molecular profiles. Some mice without any change in the three tests are considered as resilience, i.e., their invulnerability to CUMS, which were not included in our study. The mechanisms underlying vulnerability and invulnerability to the stress are not our current topic, which will be studied in the future. As 30% of CUMS-treated mice met the depression criteria and all CUMS-treated mice did not show a change of the SPT at the end of week one, the stressful situations in our study were thought to be mild stress. RNA purification from medial prefrontal cortices Twenty-four hours after the mice were surely tested to demonstrate CUMS-induced depression-like behaviors, these depression-like mice and controls were anesthetized by using Isoflurane and decapitated by guillotine. Both sides of the medial prefrontal cortices were quickly isolated and dissected on ice-cold glass slides. These cortical tissues were placed into the frozen vials that contained RNAlater RNA Stabilization Reagent (QIAGEN, Germany) at 4°C and stored at -80°C for subsequent analyses. Total RNAs from the tissue of medial prefrontal cortex in each mouse were isolated with TRIzol Reagent (Life Technologies, Carlsbad, CA, USA) based on the manufacturer instruction. Total RNA samples placed in the dried ice were delivered to Beijing Genomics Institute (BGI), China for sequencing analysis. RNA samples were done with quality-control by BGI staff using Agilent 2100 Bioanalyzer (Agilent Technologies, USA) with RNA 6000 nano Reagents Port 1. The concentrations of the total RNA, the values of RNA integrity number (RIN) and the ratios of 28S to 18S ribosomal RNA were measured. 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.0 were selected for the construction of transcriptome and small RNA libraries, respectively. RNA sequencing mRNAs were extracted from the total RNA by oligo(dT) beads. mRNAs were randomly sheared into 200 bp fragments that were reversely transcribed into the complementary DNA (cDNA) by random oligonucleotides. Upon end repair, these synthesized cDNAs were purified by using QiaQuick PCR extraction kit and ligated by sequencing adaptors. After the amplification with Illumina PCR Primer Cocktail within 15 cycles of PCR reaction, these cDNAs were size-selected and purified by agarose gel electrophoresis. cDNAs with the sizes between 200 and 300 bp were selected for the library construction. In the meantime, the miRNA sequencing library was constructed from those total RNAs, in which low molecular weight RNAs (18–30 nt) were isolated by the polyacrylamide gel electrophoresis. 5ˊ-RNA adapter was ligated to RNAs with T4 RNA ligase. The ligated RNAs were size-fractionated and 36–50 nucleotide fractions were excised. 3ˊ-RNA adapter was subsequently ligated to the precipitated RNA by using T4 RNA ligase. The ligated RNAs were size-fractionated and the 62–75 nucleotide fraction (small RNA+adaptors) was excised. Small RNAs ligated with adaptors were subjected to RT-PCR to produce sufficient templates for the sequencing. PCR products were purified and collected by gel purification and ready for the high-throughput sequencing. The qualities of transcriptome and small RNA libraries were assessed by using the Agilent 2100 Bioanalyzer (Agilent Technologies, CA USA). The quantity of PCR products was verified by quantitative PCR in ABI StepOnePlus Real-Time PCR System. After the quality and quantity of the prepared libraries were confirmed to be qualified, their sequencings were done by using Illumina HiseqTM 2500 platform (Illumina Inc., San Diego, CA USA). The average reading length of two libraries were about 100 bp (pair-end) and 49 bp (single-end), respectively. Bioinformatics for mRNA The original image data was transferred into the sequence data by base calling, which were defined as raw data or raw reads. DynamicTrim Perl script implemented in SolexaQA package was performed to control the quality of raw sequencing data based on the following criteria. 1) Remove the reads with adapters. 2) Remove the reads in that the unknown bases were more than 10%. 3) Remove the reads with 50% of the bases with low quality score (PHRED score 5). After filtering, the remained reads were called as the "clean reads" and then 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. For the alignment and mapping, the maximum of allowable mismatch was set to 3 for each read. The sole reads uniquely aligned to the genes were used to calculate gene expression level. 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. We screened the differentially expressed genes (DEGs) based on NOIseq package method, which determined DEGs between two groups with the biological replicates. The threshold used to identify DEGs was fold-change larger than 1.5 and the diverge probability larger than 0.8. We subsequently conducted a pathway enrichment analysis of DEGs’ association with physiological or biochemical processes. In these enrichment analyses, we used the hypergeometric test implemented in the tool WebGestalt (version 2) and the canonical pathways from the Kyoto Encyclopedia of Genes and the Genomes (KEGG) database. This analysis would identify the enriched metabolic pathways or signal transduction pathways in DEGs that were compared with whole genome background. P-values from the hypergeometric tests were adjusted by Benjamini-Hochberg method. Those pathways with adjusted p-values less than 0.05 were considered to be significant enrichments. Bioinformatics for miRNA 49 nt-sequence tags from Hiseq sequencing were initially processed to remove adaptor sequences, low-quality reads and contaminants for the credible clean reads. To remove the reads from ncRNA (noncoding RNA), such as rRNAs (ribosomal RNAs), tRNAs (transfer RNAs), snRNAs (small nuclear RNAs), snoRNAs (small nucleolar RNAs) and repeat RNA, we aligned the reads to the Genbank database and Rfam database with blast or bowtie softwares. All of the high-quality clean reads ranging from 18–25 nt were matched to the known miRNA precursor of corresponding species in miRBase to obtain the miRNA count. The detailed criteria include: 1) align the tags to miRNA precursor in miRBase with no any mismatch, 2) based on the first criteria, the tags align to the mature miRNA in miRBase with at least 16 nt overlap allowing offsets. Those miRNAs satisfied with these criteria would be counted to get the expression of identified miRNAs. The remained reads without any annotation were used to predict the potential novel miRNAs and its stem loop structure by Miredp based on the references [[65]44, [66]45]. To