Gene expression profile in peripheral blood mononuclear cells of postpartum depression patients

Postpartum depression (PPD) is a common mental health problem that causes maternal suffering and various negative consequences for offspring. The pathogenesis of PPD and the causes of consequences for offspring remain largely unknown. Here, we applied RNA sequencing to sequence the whole transcriptomes of peripheral blood mononuclear cells (PBMCs) from PPD patients (Edinburgh Postnatal Depression Scale [EPDS] score ≥13) and control subjects (EPDS = 0). We found that PPD was positively correlated with multiple genes involved in energy metabolism, neurodegenerative diseases and immune response, while negatively correlated with multiple genes in mismatch repair and cancer-related pathways. Remarkably, genes associated with appetite regulation and nutrient response were differentially expressed between PPD and control subjects. Then, we employed a postnatal growth retardation model by repeated immobilization stress (IS) stimulation to maternal mice. The expression of appetite regulation and nutrient response-related genes in the PBMCs of IS mice and in the hypothalamus of their offspring were also affected. In conclusion, this study provides a comprehensive characterization of the PBMCs transcriptome in PPD and suggests that maternal stress may affect appetite regulation and nutrient response in the hypothalamus of offspring mice.

Multiple pathways were altered in PPD samples. Gene Set Enrichment Analysis (GSEA) was then performed to investigate functional associations of gene expression changes in the PBMCs samples (PPD and normal control). We generated a gene list with greatest changes using RNA-seq data, and the enrichment of Kyoto Encyclopedia of Genes and Genomes (KEGG) biological pathways was evaluated by GSEA. GSEA analysis indicated that 41 and 11 pathways were significantly altered in PPD and control samples, respectively, with P value less than 0.05 and FDR less than 0.25 (Table 1). Remarkably, PPD was positively correlated with multiple genes in energy metabolism (oxidative phosphorylation, pyruvate metabolism, glycolysis/gluconeogenesis, sphingolipid metabolism, galactose metabolism, ether lipid metabolism, tryptophan metabolism, citrate cycle/TCA cycle and adipocytokine signaling pathway), neurodegenerative diseases (Parkinson's, Huntington's, Alzheimer's and Prion disease) and immune response (antigen processing and presentation, chemokine signaling pathway, cytokine/ cytokine receptor interaction, Nod like receptor signaling pathway and Toll like Receptor signaling pathway), while negatively correlated with multiple genes in mismatch repair and cancer-related pathways. Eight (Table 2). To confirm the results of RNA-seq, real-time PCR was performed to detect the mRNA expression of DUSP1, INSR, RXRA, ADRB3, CNR1 and PPARG in PBMCs from control (n = 27), PPD patients with lower EPDS score (n = 28) and PPD patients with higher EPDS score (n = 28). The variation trend was consistent with RNA-seq results. As shown in Fig. 1, mRNA levels of DUSP1 and CNR1 had significant difference between higher-EPDS-score group and control group, while comparable mRNA levels were observed in lower-EPDS-score group comparing with control group. INSR mRNA levels were higher in PPD groups compared to control group, and no significant difference was observed between higher-EPDS-score group and lower-EPDS-score group. RXRA mRNA levels were progressively increased, and mRNA levels of ADRB3 and PPARG were gradually decreased with the increasing of EPDS score. No obvious difference was observed in the p-values between ADRB3 (RNA-seq FDR = 1.61E-23) and other genes (RNA-seq FDR > 0.2) ( Table 3).

Confirmation of expression measurements with real-time PCR and enzyme-linked immunosorbent assay (ELISA) analyses.
Genes encoded chemokines and cytokines, including IL1B, C-X-C Motif Chemokine Ligand 2/3/16 (CXCL2/3/16) and C-C Motif Chemokine Ligand 3/24/28 (CCL3/24/28) were identified up-regulated genes by RNA-seq (Table 4). The altered levels of cytokines/chemokines in the serum have been reported in various human diseases, such as hear failure 22 , cancers 23 and systemic lupus erythematosus 24 . We then performed ELISA assays to detect the protein levels of IL1β, CXCL2 and CXCL3 in the serum form control (n = 27), PPD patients with lower EPDS score (n = 28) and PPD patients with higher EPDS score (n = 28). Figure 2 revealed the serum concentrations of IL1B, CXCL2 and CXCL3 were significantly increased in both PPD groups compared to control group.
RNA sequencing analysis on a mouse model. To further explore how PPD affects weight gain of children, we established a postnatal growth retardation model by repeated immobilization stress (IS) stimulation to maternal mice. At 1 week and 3 weeks after birth, the offspring mice were weighed. Maternal IS significantly decreased the body weight of their offspring (1 week old: Control, 4.94 ± 0.64 g; IS, 3.90 ± 0.66 g; 3 weeks old: Control, 14.12 ± 2.82 g; IS, 10.44 ± 1.90 g, P < 0.0001) within 3 weeks. We then performed RNA-seq on PMBCs from three pairs of control mice and IS mice. A total of 3,327 significantly DEGs with 1,555 up-regulations (Supplementary Table 3) and 1,772 down-regulations (Supplementary Table 4) were identified in IS mice, when compared with control mice (P < 0.05 and fold change > 1.2). Of these DEGs, 1817 genes (840 up-regulations and 977 down-regulations) were with Fold change > 1.5. We found that 55 DEGs were up-regulated (including appetite regulation and nutrient response-related genes, IL1B, DUSP1 and RXRA), while 36 DEGs (including CCND1) were down-regulated in PBMCs samples both from PPD patients and IS mice (Fig. 3). Maternal IS altered gene expression levels in hypothalamus of offspring mice. The appetitesatiety centers in the hypothalamus are also influenced by stress 25 . We hypothesized that the expression of appetite regulation and nutrient response-related genes in the hypothalamus of offspring mice were also affected by maternal IS. As shown Fig. 4, Dusp1, Insr, Rxra and Il1b significantly increased in the IS group, while Adrb3, Cnr1 and Pgarg notably decreased in IS group as compared to the Control group. The change trends were consistent with the findings in the PPD patients.

Discussion
Maternal postpartum depression has become a significant public concern because of its increasing prevalence, high severity and great impact on both children and mothers. Although there are numerous studies examining the prevalence and correlates of PPD, the pathophysiology and the causes of child consequences of PPD remains largely obscure. In this study, we conducted the transcriptome sequencing for PBMCs from PPD patients and  Table 2. DEGs associated with appetite regulation and nutrition level.   controls. We analyzed the expression difference at gene levels between PPD and controls and identified 2,164 DEGs. Some of the DEGs were identified as differentially expressed in a previous report such as Cyclin B1 (CCNB1), Early Growth Response 2 (EGR2), Polo Like Kinase 1 (PLK1) and Spermatogenesis Associated 20 (SPATA20) 20 , while most were novel genes. Functional analyses indicated that 41 and 11 KEGG pathways were significantly altered in PPD and control samples, respectively. Remarkably, PPD was positively correlated with multiple genes in energy metabolism, neurodegenerative diseases and immune response, while negatively correlated with multiple genes in mismatch repair and cancer-related pathways. A previous microarray study has found a signature of 73 genes in PBMCs (fold change > 1.5) between PPD patients and controls 20 , and found differences in transcription and immune activation. Here, we have identified more DEGs (621 with fold change > 1.5) and more pathways associated with PPD, suggesting RNA-seq is a powerful and sensitive tool for expression profiling. We revealed the effects of maternal PPD on the weight gain of infants, which were in line with previous reports [6][7][8] . Stressful life events strongly correlated with the onset and progression of PPD 2 . Repeated immobilization stress (IS) stimulation has been widely used to study major depression 26,27 . Here, repeated maternal IS stimulation 28 significantly decreased the body weight of their offspring within 3 weeks. RNA-seq was then performed on PBMCs from maternal mice samples. A total of 55 up-regulated genes and 36 down-regulated genes were identified in both samples from PPD patients and IS mice. Ray S et al. had reported the gene expression changes in the maternal brain during pregnancy and the postpartum period by using RNA-seq. 29 . They identified differential expression genes during pregnancy and the postpartum period implicated in PPD and depressive disorder. Among the DEGs identified in our present study, 2 DEGs implicated in PPD and 28 DEGs implicated in depressive disorder have been reported by Ray S et al. 29 (Table S6). The appetite-satiety centers in the hypothalamus are influenced by stress 25 . Glucocorticoid treatment lowered the gain in body weight of rats 30 and Dusp1 is a glucocorticoid target gene in rat hypothalamus 31 . Mice with a neuron-specific disruption of Insr gene had increased body fat 32 . Il1b is an anorectic gene and Cnr1 is an orexigenic gene in the hypothalamus 33 . PPARγ activation in the brain results in increased food intake and weight gain 34 . Here, qRT-PCR results showed that maternal IS increased the expression of Dusp1, Insr and Il1b, but reduced the expression of Cnr1 and Pgarg in the hypothalamus of offspring mice. These findings suggested that maternal stress may affect the weight of offspring by regulating the appetite. Adrb3 is an energy expenditure gene in the hypothalamus 33 . Rxra, an important adipogenesis regulator, has been found expressed in the central nervous system of mice 35 and its function in the brain is unclear. Here, maternal IS treatment led to decreased Adrb3 expression and increased Rxra expression in the hypothalamus of offspring, which also suggested the regulatory role of maternal stress on energy balance and nutrient response of their offspring. Further investigation are needed to explore the detailed mechanisms how maternal stress affected the expression of these genes.
In summary, RNA-based approach provides a comprehensive gene expression profiling in PMBCs of PPD. Gene-set enrichment analysis further identified differences in biological pathways relating to energy metabolism, neurodegenerative diseases, immune response, mismatch repair and cancer-related pathways between the PPD  Table 4. DEGs associated with immune response. and control groups. The IS mouse model experiments suggest that maternal stress may regulate the energy balance and nutrient response of their offspring, thus reducing the body weight.

Materials and Methods
Patients and sample collection. The study protocol was approved by the Institutional Review Board of Jinshan Hospital of Fudan University, and written informed consent was obtained. The study was carried out according to the relevant guidelines. Depressive symptom scores were evaluated using the Edinburgh Postnatal Depression Scale (EPDS) instrument 36 on six weeks after delivery. Mothers with EPDS score ≥ 10 (n = 56) were included in this study. Mothers with EPDS score = 0 (n = 27) after delivery and showing no depressive symptoms were set as control subjects. The control subjects were matched for age, sex, and BMI with the PPD subjects. Subjects were not included if they had a history of past or present psychiatric diagnoses, had medical or neurological illness, or had used antidepressants or other psychotropic medications. Peripheral blood samples were collected by venipuncture from study subjects after delivery. Ten paired samples were randomly selected for RNAsequencing, and the remaining samples were subjected to real-time PCR and ELISA analyses.

Animal experiments. The animal study was approved by Animal Care Institutional Review Board of Fudan
University and performed in accordance with the relevant guidelines and regulations. ICR mice (CLEA Japan, Inc., Tokyo, Japan) were kept at 21 ± 1 °C under a 12/12 h light-dark cycle. After 1 week of acclimatization, one female was housed with two males in each cage for 4 days until a copulation plug was found. On postpartum day 1, immobilization stress (IS) was performed by transferring mice from their home cage to a small meshed cylinder (3 × 6 cm). IS was applied to the animals for 3 hours every day for 3 weeks. The mice were sacrificed by cervical dislocation, blood was collected from the right ventricle, and the hypothalamus was dissected out as previously described 37 and kept at −80 °C until use.
RNA extraction, sequencing and data analysis. PBMCs were separated from collected whole blood by centrifugation at 300 g for 30 minutes with mononuclear/polynuclear cell-resolving medium (Flow Laboratories, Rockville, MD, USA). Total RNA was isolated with the Trizol reagent following the manufacturer's instructions (Invitrogen, Carlsbad, CA, USA). The RNA yields were quantified by NanoDrop ND1000 (Thermo-Fisher , which captured polyA-containing mRNAs, with 2 μg of total RNA. RNA single-end sequencing was performed using Illumina Genome Analyzer II using the standard protocol. Sequenced reads were base-called by using the Illumina standard pipeline. The cleaning reads were mapped to the human genome (hg19) using TopHat v2.0.11 with the default options with a TopHat transcript index built from Ensembl GRCh37. Count files of the aligned sequencing reads were gene rated by the htseq-count script from the Python package HTSeq with union mode, using the GTF annotation file. The read counts from each sequenced sample were combined into a count file, which was subsequently used for the differential expression analysis. Differential analyses were performed to the count files using DESeq. 2 packages, following standard normalization procedures 38 . Genes with less than 5 total counts in both conditions were removed from further analysis. Pathway analysis. We identified DEGs between PPD and control samples based on the following criteria: P < 0.05 and fold change > 1.2. Gene Set Enrichment Analysis (GSEA) software was used to determine whether a total of 178 gene sets from KEGG showed statistically significant, concordant differences between PPD patients (n = 8) and normal controls (n = 10) as describe previously 39 . Differentially expressed gene sets/pathways were identified with threshold false discovery rate (FDR) less than 0.5 and P less than 0.05.
Quantitative real time polymerase chain reaction (RT-PCR). PBMC RNA was extracted from 27 control subjects, 28 PPD patients with lower EPDS score (10 ≤ EPDS < 13) and 28 PPD patients with higher EPDS score (≥13) and reverse transcribed with RevertAid First Strand cDNA Synthesis Kit (Thermo-Fischer Scientific). Gene expression was conducted on Applied Biosystems 7300 instrument (Applied Biosystems, Foster City, CA, USA) using SYBR-green PCR Master Mix (Thermo-Fischer Scientific). Gene expression values were calculated with the ΔΔ Ct method 40 and GAPDH was served as reference gene. ANOVA analysis was used to calculate the statistical significance of difference among groups. Primers were listed in Supplementary Table 5. ELISA analysis. Serum was separated from collected whole blood by centrifugation at 2,500 rpm for 15 min at 4 °C. Serum concentrations of IL-1β, CXCL2 and CXCL3 were determined with ELISA assay (Bio-Swamp life science, Shanghai, China) following the instructions of the manufacturer. Absorbance was read at 450 nm using a microplate reader (Bio-Rad Laboratories Inc., Hercules, CA, USA). ANOVA analysis was performed to calculate the statistical significance of difference among groups.