Chronic stress induces significant gene expression changes in the prefrontal cortex alongside alterations in adult hippocampal neurogenesis.

Adult hippocampal neurogenesis is involved in stress-related disorders such as depression, posttraumatic stress disorders, as well as in the mechanism of antidepressant effects. However, the molecular mechanisms involved in these associations remain to be fully explored. In this study, unpredictable chronic mild stress in mice resulted in a deficit in neuronal dendritic tree development and neuroblast migration in the hippocampal neurogenic niche. To investigate molecular pathways underlying neurogenesis alteration, genome-wide gene expression changes were assessed in the prefrontal cortex, hippocampus and the hypothalamus alongside neurogenesis changes. Cluster analysis showed that the transcriptomic signature of chronic stress is much more prominent in the prefrontal cortex compared to the hippocampus and the hypothalamus. Pathway analyses suggested huntingtin, leptin, myelin regulatory factor, methyl-CpG binding protein and brain-derived neurotrophic factor as the top predicted upstream regulators of transcriptomic changes in the prefrontal cortex. Involvement of the satiety regulating pathways (leptin) was corroborated by behavioural data showing increased food reward motivation in stressed mice. Behavioural and gene expression data also suggested circadian rhythm disruption and activation of circadian clock genes such as Period 2. Interestingly, most of these pathways have been previously shown to be involved in the regulation of adult hippocampal neurogenesis. It is possible that activation of these pathways in the prefrontal cortex by chronic stress indirectly affects neuronal differentiation and migration in the hippocampal neurogenic niche via reciprocal connections between the two brain areas.

Introduction

Adult hippocampal neurogenesis is thought to play an important role in depression and antidepressant effects (Miller and Hen, 2015; Egeland ; Boldrini ), however its exact role in the course of disease neurobiology is unknown. Chronic stress has been identified as a significant risk factor for depression (Bernet and Stein, 1999; Fava and Kendler, 2000; Pawlby ). Unpredictable chronic mild stress (UCMS) is widely used as a valid model of depression. UCMS induces such behavioural endophenotypes as anhedonia, anxiety, decreased grooming and weight loss (Ducottet ; Mutlu ; Nollet ). However some groups report hyperactivity and anxiolytic effect of this paradigm (Pothion ; Schweizer ; Strekalova ). UCMS is also associated with a decline in hippocampal neurogenesis. While multiple hypotheses have been put forward to explain the mechanism of neurogenic decline (Egeland ), a full coherent clinically relevant picture has not yet been formed. One of the ways to study neurobiological mechanisms in their full complexity is to utilize hypothesis-free explorative approach of a genome-wide transcriptome changes in selected tissues and brain regions. This approach is currently being widely utilized in clinical depression research to discover novel genes and pathways involved in the neurobiology of depression (Hervé ; Cattaneo ).

The majority of the research on depression and chronic stress focuses on the brain areas of the cortico-limbic circuit, including but not limited to the hippocampus (HIP), hypothalamus and the prefrontal cortex (PFC). Human neuroimaging studies confirm the involvement of these regions in depression (Lorenzetti ; van Tol ). The role of PFC in depression is further supported by the evidence from the deep brain stimulation trials, showing that deep brain stimulation over PFC and anterior cingulate cortex produces an antidepressant effect (Levkovitz ; Berlim ). Investigation of PFC pathology in animal models of depression demonstrated that exposure to chronic stress and increased glucocorticoid levels cause damage and neurodegenerative changes in the PFC, such as appearance of reactive microglia (Hinwood ), atrophy of pyramidal neurons (Cerqueira, 2005), dendritic atrophy (Liston ; Dias-Ferreira ) and reduction of synaptic proteins expression and synaptic currents (Li ; Müller et al., 2011).

Yet the majority of transcriptomic studies of chronic stress effects in rodents focussed on the consequence of stress on gene expression in the HIP. The HIP has been shown to be particularly sensitive to stress, potentially due to a high expression of the glucocorticoid receptor activated by corticosteroid hormone released in the course of stress response. Both rat and mouse studies find activation of immune response-related gene networks and signalling pathways often centred on nuclear factor kappa B (Malki , 2015; Gray ). A number of other studies investigating the effect of chronic mild or restraint stress on hippocampal gene expression in rats and mice identified signalling pathways involved in cell fate, such as cell apoptosis (Bergström ), proliferation and cell cycle control (Liu ), axonal guidance and cell migration (Jungke ; Datson ; Zitnik ).

Another region of interest in the chronic stress research is the hypothalamus. As a central organ in the hypothalamo-pituitary-adrenal axis regulating the systemic level of glucocorticoids, hypothalamus plays a central role in the endocrinal stress response of the body. Hypothalamus also receives neuroendocrine input from the gut metabolic hormones ghrelin and leptin (Lep) and is connected with the reward-related ventral tegmental area (van Zessen ). Another hypothalamic area implicated in chronic stress response is the suprachiasmatic nucleus, a major regulator of the circadian clock, relevant for sleep disturbance and alteration of molecular circadian rhythms observed in depression (Li ).

For these reasons, the three limbic circuit brain areas described above were selected in this study to analyse genome-wide transcriptomic changes induced by UCMS exposure to find novel potential mechanisms of stress-induced neurogenic decline using a hypothesis-free approach.

Materials and methods

Experimental design

Seven weeks old male BALB/cAnNCrl mice were obtained from Charles River (Margate, Kent, UK). All mice were housed in the Biological Services Unit at the Institute of Psychiatry, Psychology and Neuroscience in standard conditions [19–22°C, humidity 55%, 12 h:12 h light:dark cycle with lights on at 7.30 am, food (Rat and Mouse No. 1 Diet; Special Diet Services, Essex, UK) and water ad libitum]. Mice were allowed to acclimatize to new housing conditions for 1 week before experimental procedures commenced. All housing and experimental procedures were carried out in compliance with the local ethical review panel of King’s College London under a UK Home Office project licence held in accordance with the Animals (Scientific Procedures) Act 1986 and the European Directive 2010/63/EU.

Two cohorts of 8 weeks old male BALB/c mice were exposed to UCMS or control (CNTRL) conditions. Mice were randomly allocated to experimental groups. The BALB/c strain was selected for this study due to its high sensitivity to stress and high immune system reactivity (Potter, 1985; Kim ; Surget and Belzung, 2009). The first cohort was sacrificed for tissue collection after 6 weeks of UCMS exposure to examine the behavioural changes and the state of the hippocampal neurogenic niche (n = 10 per group). As the behavioural data from Cohort 1 demonstrated the strongest behavioural phenotype at Week 4, Cohort 2 was sacrificed for tissue collection after 4 weeks of UCMS exposure to detect UCMS-induced gene expression events potentially preceding changes in the adult hippocampal neurogenesis (n = 8 per group). For graphical representation of the study design, see Supplementary Fig. 1.

UCMS protocol and behavioural assessments

The UCMS protocol used was based on the methodology described by Nollet . The UCMS group was exposed to four daily stressors in a random-like order repeated every week. In order to avoid habituation, the longevity of stressors varied from 30 min on Week 1 to 4 h towards the end of the UCMS protocol. For the schedule of stressors included in the protocol, see Supplementary Methods  Table 1.

Table 1

Canonical pathways deemed significant in the PFC dataset by IPA

Ingenuity canonical pathways	−log (P-value)	Ratio	Z-score	Associated differentially expressed genes
Axonal guidance signalling	2.61	0.05		Ephb2, Tuba4a, Gng13, Sema4f, Robo3, Gng7, Tubb2b, Ephb6, Sema6d, Mag, Wnt10a, Efna5, Arhgef6, Lingo1, Mras, Myl4, Fzd5, Sema7a, Adamts4
Glutamate receptor signalling	2.12	0.10		Slc17a7, Slc17a6, Homer1, Grm4, Gng7
Calcium signalling	2.03	0.06	2.45	Camk2a, Tnnc1, Myh3, Rcan3, Myl4, Mef2c, Itpr1, Camkk2, Camk2g
Wnt/β-catenin signalling	2.01	0.06	0.33	Csnk2a2, Nlk, Wnt10a, Dkk3, Sox10, Fzd5, Acvr2b, Sox11, Tcf7l2
GABA receptor signalling	1.94	0.09		Slc32a1, Gabrg1, Gad1, Mras, Ap2s1
Paxillin signalling	1.65	0.07	0.45	Actn2, Actb, Arhgef6, Mras, Itgb4, Mapk11
Thrombin signalling	1.64	0.05	1.89	Camk2a, Arhgef6, Mras, Myl4, Gng13, Itpr1, Mapk11, Gng7, Camk2g
Gαq signalling	1.63	0.06	1.13	Napepld, Adrbk1, Mras, Gng13, Arhgef25, Itpr1, Chrm3, Gng7
Ephrin B signalling	1.50	0.07		Ephb6, Ephb2, Mras, Gng13, Gng7
RhoGDI signalling	1.49	0.05	−0.82	Dgkz, Arhgdig, Actb, Arhgef6, Mras, Myl4, Gng13, Gng7
B-cell receptor signalling	1.40	0.05	2.12	Synj2, Map3k10, Camk2a, Inpp5f, Mras, Mef2c, Mapk11, Camk2g
CCR5 signalling in macrophages	1.35	0.08		Mras, Gng13, Mapk11, Gng7
G protein signalling mediated by tubby	1.34	0.09		Mras, Gng13, Gng7
Triacylglycerol degradation	1.33	0.14		Faah, Mgll
α-Adrenergic signalling	1.31	0.06		Adra2a, Mras, Gng13, Itpr1, Gng7
Glutamate-dependent acid resistance	1.30	0.50		GAD1

Canonical pathways deemed significant in the PFC dataset ordered by the log (P-value) derived from Fisher’s exact test, with dataset genes/pathway genes ratio and Z-score of predicted pathway up- or downregulation, calculated for pathways where differential gene expression showed consistent direction of change.

Animals’ weight and coat state were assessed during UCMS on a weekly basis. Coat score was determined on a scale of 0–7 according to Nollet by a researcher blinded to experimental conditions. Sucrose preference test was conducted during and at the end of the UCMS exposure. The following additional tests were conducted at the end of the UCMS exposure: open field, novelty suppressed feeding (NSF), splash test and Porsolt swim test (PST). For the details of experimental procedures, see Supplementary Methods.

Blood collection

For cytokine and corticosterone analysis, blood was collected by incision method from the lateral tail vein (Sadler and Bailey, 2013) 24 h before and 30 min after the PST. Whole blood (30–50 µl) was collected into ethylenediaminetetraacetic acid covered (EDTA) microvette CB 300 tubes (Sarstedt, Leicester, UK) and separated by centrifugation for 10 min at 3000 rpm (4°C). Blood was also collected by the cardiac puncture after terminal general anaesthesia has been induced. In total, 100–200 μl of whole blood was collected using syringe injected into the cardiac cavity and subsequently processed as described above. This blood draw was used for cytokine, C-reactive protein (CRP) and Lep measurements.

Corticosterone measurement

Plasma corticosterone levels were measured using a commercially available corticosterone enzyme-linked immunosorbent assay kit (Enzo Life Sciences, Lausen, Switzerland) according to the manufacturer’s instructions for small sample volume. All samples were analysed in technical duplicates. Obtained data were analysed using five-parameter logistic curve analysis using the My Assays online data tool, MyAssays Ltd, http://www.myassays.com/five-parameter-logistic-curve.assay (last accessed 5 October 2020).

Luminex

The level of cytokines in the plasma was determined using the multiplex screening assay based on magnetic Luminex® xMAP® technology as described in Hye . For this assay, custom-made pre-mixed multianalyte kit was purchased from RnD systems, Minneapolis, USA (catalogue N LXSAMSM). This kit contained multi-coloured magnetic microparticles pre-coated with antibodies to 10 selected analytes (granulocyte-monocyte colony stimulating factor, interleukins IL-1β, IL-2, IL-4, IL-5, IL-6, IL-10, tumour necrosis factor-alpha and interferon gamma, CRP and Lep). To measure the fluorescent signal, Luminex® 100/200TM system was used. All assay steps were conducted according to the manufacturer instructions. Due to the small volume of blood collected, technical replicates for plasma samples were not included. Obtained data were analysed using five-parameter logistic curve analysis using the My Assays online data tool, MyAssays Ltd, http://www.myassays.com/five-parameter-logistic-curve.assay (last accessed 5 October 2020).

Brain tissue collection

To collect fixed tissue, animals were anaesthetized with Euthatal i.p. (Merial Animal Health Ltd, Harlow, UK) at a dose of 40 mg/kg pentobarbital sodium. Deep anaesthesia was confirmed by the loss of righting and pain reflexes and slowing of the rate of respiration. Animals were subsequently transcardially perfused with 30 ml of saline and 50 ml of 4% paraformaldehyde (Paraformaldehyde, prills, 95%, 441244, Sigma-Aldrich, Poole, UK) in phosphate-buffered saline (pH = 7.4 from phosphate-buffered saline tablets, 18912-014, Gibco by Life Technologies, Paisley, UK) through the left cardiac ventricle. Brains were post-fixed overnight in 4% paraformaldehyde at 4°C for 24 h, and subsequently stored at 4°C in 30% sucrose (Sigma-Aldrich, Poole, UK) in phosphate-buffered saline until sectioning (1–3 weeks).

For fresh frozen tissue collection, following decapitation brains were removed from the skull and placed dorsal side down on a wetted filter paper on a petri dish kept on ice. For the dissection of the brain areas technique, see Supplementary Methods. Right and left dissected areas were pooled for each sample, while brain tissue from each mouse was analysed as individual samples. Dissected brain areas were immediately flash frozen on dry ice.

Immunohistochemistry

Serial coronal sections of 40-μM thickness were cut on a HM430 sliding freezing microtome (Thermo Scientific) and stored in Tris-buffered saline (pH = 7.4) with 30% v/v glycerol, 15% w/v sucrose and 0.05% w/v sodium azide (all Sigma-Aldrich, UK). For immunohistochemistry on free floating sections, every sixth section was used. Immunohistochemistry was conducted as previously described (Dias ). For neuroblast detection, anti-doublecortin (anti-DCX) antibody (Ab18723 raised in rabbit, Abcam, UK, 1:1000 dilution) was used; for microglia immunostaining, antibody for ionized calcium binding adaptor molecule 1 (Iba1; 019-19741, Wako, Japan, 1:500 dilution) was used.

Stereological analysis of immunopositive cell density

For microscopy, the Axioskop 2 MOT Plus Microscope (Zeiss) with automated stage connected to Stereoinvestigator v.7 (MBF Bioscience, USA) software was used. Immunopositive cell density was determined using the volume and cell population number estimated by semi-automated Optical Fractionator method applied by the stereological software (see Supplementary Methods). All quantifications were done by an experimenter blinded to the experimental conditions.

DCX cell classification based on dendrite morphology and their morphometric analysis

The DCX positive cells were visually classified according to the categorization of Plümpe . EF-type neurons from each section were selected for morphometric analysis of dendrites. For this, sections were imaged using a 40× objective, and micrographs were acquired using a Zeiss AxioCam MR Rev3 camera. Primary, secondary and tertiary dendrites were manually traced using the NeuronJ plugin for ImageJ developed by Meijering and measured using image pixel to millimetre calibration as has been done previously (Srivastava ; Dias ). For details, see Supplementary Methods.

Migration distance of DCX+ cells

To estimate relative distance of DCX+ cells migration through the granular zone (GZ), the height of the GZ has been divided into 10 bins, and a relative distance of migration was estimated by assigning a score from 0 to 1 to each DCX+ cell based on the position of its cell body relative to the borders of the GZ as has been done previously (Han ). The average percentage of cells positioned within each distance bin for each brain was used in statistical analysis.

Microarray and data analysis

Genome-wide gene expression was assessed in fresh frozen brain tissue using the Mouse WG-6 BeadChip Kit (Illumina, USA) (for details of RNA extraction, see Supplementary Methods). Eight brains were sampled for each treatment condition; no technical replicates were utilized. Sample labelling, hybridization and signal detection and imaging were performed according to manufacturer’s instructions. Images were analysed using the GenomeStudio software (Illumina, USA). Raw data were pre-processed using the R package Lumi (Du ). Statistical analysis of differential expression was performed using the statistical analysis of microarrays method and R-based software (Tusher ). The dataset for each brain region was analysed separately as a two-class unpaired data. For further details, see Supplementary Methods. Microarray data analysis was verified using real-time quantitative polymerase chain reaction. For gene selection, primer design and quantitative polymerase chain reaction, see Supplementary Methods. Pathway, upstream regulator and predicted functions analyses were conducted through the use of Ingenuity Pathway Analysis (IPA) (QIAGEN Inc., https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis (last accessed 5 October 2020).

Statistical analysis

Unpaired two-tailed t-test was used for between group comparisons after normal distribution of the data was confirmed with the Kolmogorov–Smirnov test for single measures. For the weekly coat state between group comparisons, Mann–Whitney U test used for each week as the data were not normally distributed. For related measures, such as DCX+ cells, dendritic morphology types and migration distance groups, two-way ANOVA was used with post hoc Bonferroni multiple comparisons within each type or distance group.

Data availability

Output of statistical analysis of microarrays for each brain region is available in the Supplementary materials for this article. Raw microarray output data are available on request.

Results

UCMS impaired coat state and grooming behaviour in exposed mice

UCMS induced significant coat state deterioration detectable from Week 4 of stress exposure (Mann-Whitney for Week 4 U = 2, P < 0.0001; for Week 5 U = 21, P = 0.023; for Week 6 U = 9, P = 0.004) (see Fig. 1A). Repeated measures two-way ANOVA showed that body weight of animals was significantly affected by week of exposure (Effect of Time F(6,108) = 127.2, P < 0.0001) and UCMS × Time interaction (F(6,108) = 5.026, P = 0.0001) (see Fig. 1B). Sucrose consumption, measured once every 2 weeks, was significantly decreased in UCMS group compared to CNTRL only on Week 4 when sucrose concentration was increased from 1% to 2% (see Fig. 1C). The CNTRL group responded to 2% sucrose with 79% preference, while the UCMS-exposed group did not prefer it over water (46%) (Effect of Time × Stress F(3,39) = 13.15, P < 0.0001). In the splash test, UCMS-exposed mice spent shorter time grooming (t(18) = 3.304, P = 0.004) (see Fig. 1D). UCMS also induced hyperactivity in the open field test (see Fig. 1E) (Mann–Whitney U = 16, P = 0.009). In the NSF test, UCMS-exposed group’s latency to feed was reduced compared to CNTRL group (unpaired two-tailed t-test t(18) = 3.227, P = 0.005) (see Fig. 1F). No differences were found in the time of immobility on the PST (see Fig. 1G).

Figure 1

Behavioural and serological parameters in the UCMS-exposed mice. Male BALB/cAnNCrl mice (n = 10/group) aged 7 weeks at the beginning of the experiment were subjected to UCMS or CNTRL conditions for 6 weeks. (A) Weekly coat state deterioration score measurements from Week 4, *P < 0.05, **P < 0.01 and ***P < 0.001 derived from Mann–Whitney U-test between CNTRL and UCMS, data represent median and interquartile range. (B) Weekly weight monitoring, *P < 0.05 derived from Bonferroni multiple comparison between CNTRL and UCMS group means at Week 1 of UCMS. (C) Sucrose consumption measured once every 2 weeks over 2 consecutive nights, **P < 0.01 derived from Bonferroni multiple comparisons. (D) Time spent grooming during the splash test. (E) Distance moved during a 5-min exposure to a dimly lit open field. (F) Latency to start eating the pellet in the NSF test. (G) Immobility in the PST. (H) Plasma corticosterone (CORT) response to the PST measured 24 h before (PRE-PST) and 30 min after (POST-PST) the test. (I) Plasma levels of the CRP in the blood collected by the cardiac puncture at the end of the behavioural testing battery (7 days after termination of the UCMS protocol). (J) Plasma levels of leptin in the blood collected by the cardiac puncture at the end of the behavioural testing battery (7 days after termination of the UCMS protocol). In D–J, *P < 0.05, **P < 0.01 and ***P < 0.001 derived from two-tailed unpaired t-test, data represent mean ± SEM.

UCMS did not affect plasma corticosterone response to acute swim stress but increased plasma CRP and reduced plasma leptin levels

To analyse the HPA axis regulation of glucocorticoid release in response to stressful stimulation, blood samples were collected 24 h before and 30 min after a 6-min long PST. Corticosterone levels were significantly elevated after PST in both groups but UCMS did not affect the baseline level of CORT or response to swim stress (Effect of PST F(1,16) = 108.6, P < 0.0001) (see Fig. 1H). None of the 10 cytokines, measured at the time of sacrifice, were detected in the plasma. However, CRP was significantly elevated in the UCMS group (t(16) = 2.508, P = 0.023) (see Fig. 1I). Lep levels, reflecting hormonal regulation of appetite, were significantly decreased in UCMS group compared to CNTRLs (t(14) = 2.477, P = 0.027) (see Fig. 1J).

UCMS reduces the density of mature DCX+ neuroblasts and the percentage of DCX+ cells residing in the subgranular zone

UCMS group presented with a reduced number of DCX+ neuroblasts in the hippocampal dentate gyrus (DG) (t(14) = 2.954, P = 0.01) (see Fig. 2D). Density reduction was specific to neuroblasts with mature-like dendritic trees reaching into the molecular layer of the DG [type ‘EF’ in Plümpe ]. See Fig. 2D for quantification and Fig. 2E for examples of each neuroblast type.

Figure 2

Effect of UCMS on the density, morphology and migration of DCX+ neuroblasts and microglial density. (A) Number of neuroblasts residing in the subgranular zone (0) and in 10 different layers of the GZ (0.1–1). (B) Examples of the DCX+ cells with a high relative migration distance in the hippocampal GZ of the UCMS-exposed mice. (C) Average relative migration distance of DCX+ cell bodies. (D) The density of all DCX+ neuroblasts (TOTAL) and of the ‘AB’, ‘CD’ and ‘EF’ types classified based on their dendritic tree morphology. (E) Examples of AB, CD and EF types of neuroblasts. (F) Representative microphotographs of the Iba1+ cells in the DG. (G) Density of the Iba1+ cells in the GZ of the DG. (H) Density of Iba1+ cells in the medial prefrontal cortex (mPFC). (I) Representative microphotographs of Iba1+ cells in the mPFC. Data presented as mean ± SEM, *P < 0.05 and **P < 0.01 derived from unpaired t-test CNTRL versus UCMS; #P < 0.05 derived from post hoc Bonferroni multiple comparisons CNTRL versus UCMS.

As adverse environment has been previously shown to affect migration of the adult-born neuroblasts (Belarbi ), the effect of the UCMS on migration distance of the DCX+ cells was analysed (Effect of Interaction F(10,170) = 1.99, P = 0.037; Effect of Distance F(10,170) = 622, P < 0.0001). UCMS reduced the percentage of DCX+ cells residing in the subgranular zone (P = 0.005), see Fig. 2B for examples and Fig. 2A for quantification. However, the average relative distance of migration of DCX+ cell bodies from the subgranular zone was not significantly different between the two groups (see Fig. 2C).

UCMS does not affect the number of Iba1+ microglia in the hippocampal GZ but increases its density in the medial PFC

Analysis of Iba1+ microglial density showed that no significant effect of UCMS on the number of microglial cells was found (see Fig. 2F and G). However, density of Iba1+ microglia was increased in the medial PFC of the UCMS-exposed group compared to CNTRL (t(16) =3.139, P = 0.006) (see Fig. 2H and I).

UCMS exerted the strongest effect on gene expression in the PFC

To elucidate molecular pathways underlying UCMS effects, bulk genome-wide gene expression analysis of the brain tissue was conducted using the gene expression microarray. For this, a separate cohort of mice (Cohort 2) was exposed to UCMS or CNTRL conditions (n = 8 per group). As behavioural monitoring during 6 weeks of UCMS exposure in Cohort 1 showed the strongest phenotype on Week 4 (see Fig. 1A and C), this time point was chosen for gene expression analysis. Behavioural assessment preceded tissue collection at Week 4. Behavioural data confirmed the presence of low grooming and hyperactivity phenotypes seen in Cohort 1, however some differences were observed in the dynamics of sucrose preference response (see Supplementary Fig. 2).

Hierarchical cluster analysis based on the top 500 variable genes showed that samples primarily clustered based on the brain region they derived from, while within the regions some treatment group-based clustering was observed. The PFC samples clustered based on the UCMS exposure, apart from two outlier samples, which were excluded from subsequent analysis as shown in Fig. 3. In the HIP, only few CNTRL and UCMS-exposed samples clustered together, and no clear distinct pattern of expression could be seen on the heatmap for the hippocampal (HIP) UCMS and CNTRL samples. In the hypothalamus, no clustering based on stress exposure was observed (see Fig. 3).

Figure 3

Heatmap of cluster analysis based on the top 500 variable genes. Adult male BALB/c mice were exposed to UCMS or control (CON) for 4 weeks (n = 8/group) after which fresh frozen brain tissue from selected brain regions [whole HIP, hypothalamus (HYP) and PFC] were subjected to a genome-wide transcriptomic analysis using Illumina microarray platform. The heatmap shows the clustering of samples (top connecting lines) and the heatmap pattern of expression based on intensity values of the top 500 variable genes. The gradient of green and red represents the deviation of each sample intensity value for this gene from the mean intensity across all samples. Individual sample IDs and their groups highlighted at the bottom of the heatmap.

Differential expression analysis using statistical analysis of microarrays package yielded the highest number of significant genes in the PFC

Significance analysis of microarrays method detected 467 differentially expressed genes (DEGs) between the UCMS and CNTRL groups in the PFC with false discovery rate of 0. For the full list of top up- and downregulated genes selected based on their fold changes, see Supplementary Table 4. As such a low level of false discovery rate was calling a very small number of genes in other regions, false discovery rate level was relaxed to a recommended 0.095 for the HIP and the hypothalamus. Such settings called 46 and 27 DEGs in the HIP and in the hypothalamus, respectively (see Supplementary Table 5).

Pathway analysis revealed meaningful pathways, functions and predicted upstream regulators in the PFC dataset

Canonical pathway analysis conducted by the IPA™ software mapped the DEGs on known canonical pathways. The software identified a total of 240 canonical pathways associated with the DEGs. Among these, 17 were statistically significant (P < 0.05) based on Fisher’s exact test analysis utilized by the IPA (see Table 1). The software also quantified the ratio of the number of DEGs involved in each pathway to the total number of genes comprising the pathway. Table 1 lists all DEGs involved in each of the significantly activated pathways. Based on the fold changes of the DEGs, IPA also made a prediction of whether the pathway would be up (positive Z-score) or down (negative Z-score) regulated and its magnitude (Z-score value). Due to conflicting direction of DEGs change in some pathways, it was not always possible to make such a prediction (missing Z-score). Interestingly a pathway directly related to axonal maturation, the axonal guidance pathway was identified as the most significant. The list also included pathways previously shown to be involved in neurogenesis regulation such as Wnt/β-catenin and Ephrin signalling. Common pathways of neuronal signalling such as glutamate, calcium and GABA receptor signalling were also identified. Interestingly, immune response-related B-cell receptor pathway also appeared to be significantly activated.

Next, the predicted upstream regulators analysis was conducted using the IPA software. Upstream regulator analysis identifies molecules, which activation or suppression can explain gene expression changes in the dataset, as well as regulators’ own interactions with each other. The software utilizes a broad definition of the upstream regulator for this analysis, so that any type of regulatory molecule, from microRNA to a drug, can be included in the results. IPA identified 63 potential upstream regulators of DEGs, with 59 being significant according to Fisher’s exact test. Table 2 lists the top 15 predicted upstream regulators with P-value ≤ 0.01, connected to two or more DEGs. In some cases, the software could also make a prediction of the direction and magnitude of up- or downregulation based on the fold changes data available for target genes (activation Z-score).

Table 2

Upstream regulators of differentially expressed genes in the PFC dataset predicted by the IPA software

	Upstream regulator	Gene name	Activation Z-score	P-value	No. of target genes
1	Htt	Huntingtin	−0.004	0.0000	25
2	Lep	Leptin	−0.970	0.0000	8
3	Myrf	Myelin regulatory factor		0.0001	3
4	Mecp2	Methyl-CpG binding protein 2	1.126	0.0001	7
5	Bdnf	Brain-derived neurotrophic factor	−0.181	0.0001	9
6	Hdac4	Histone deacetylase 4		0.0001	8
7	Kmt2a	Lysine methyltransferase 2A	1.342	0.0002	5
8	Ntrk2	Neurotrophic receptor tyrosine kinase 2		0.0011	3
9	Dcc	DCC netrin 1 receptor		0.0018	2
10	Pou4f1	POU class 4 homeobox 1		0.0022	4
11	Comt	Catechol-O-methyltransferase		0.0035	2
12	Mtor	Mechanistic target of rapamycin		0.0054	3
13	Arntl	Aryl hydrocarbon receptor nuclear translocator like		0.0117	2
14	Ascl1	Achaete-scute family bHLH transcription factor 1		0.0123	3
15	Bckdk	Branched chain ketoacid dehydrogenase kinase		0.0153	2

IPA predicted upstream regulators in the PFC dataset, P-value derived from Fisher’s exact test. Z-score reflects the direction of change where the software was able to make such a prediction.

Figure 4 highlights several key upstream regulators with the most number of target genes [huntingtin (Htt), Lep, Myelin Regulatory Factor (Myrf), methyl-CpG binding protein 2 (Mecp2) and brain-derived neurotrophic factor (Bdnf)] at the centre of the regulatory network, as well as the predicted complexity of their regulatory relationships with the DEGs and among themselves. According to the analysis, Htt was identified as the leading regulator, with 25 DEGs from the PFC dataset listed as its target genes. Lep, Mesp2, Bdnf and Hdac4 followed Htt with seven to nine target genes among the DEGs. Some of the top upstream regulators could be grouped into three main categories: those previously implicated in neurogenesis regulation [Bdnf, mechanistic target of rapamycin (Mtor), achaete-scute homolog 1 (Ascl1)], those involved in epigenetic regulation (Mecp2 and Hdac4) and those related to some of the observed behaviours (Htt, Lep). The relationships between the main upstream regulators and their differentially expressed target genes are shown in Fig. 4. Colour coding reflects direction of fold change of the DEGs (red for upregulation and green for downregulation). Interestingly, most upstream regulators were not themselves differentially expressed, with the exception of BAF chromatin remodelling complex subunit B (Bcl11b), which was downregulated. In the HIP and the hypothalamus, the number of genes was not sufficient to produce meaningful results in the pathway analysis.

Figure 4

The predicted network of upstream regulators and their target differentially expressed genes in the PFC, designed by the IPA software based on the upstream regulator analysis. Network of upstream regulators linked to genes differentially expressed in the PFC dataset. The relative size of the upstream regulator molecules reflects their significance ranking with Htt, Lep, Myrf, Mecp2 and Bdnf being the top five predicted upstream regulators.

Functional analysis identified functions related to dendritic morphology

Next, we assessed if DEGs in the PFC dataset were linked to particular functions in the IPA database. IPA identified 185 functions related to the DEGs. The most relevant functions, based on their relation to the brain, the number of associated DEGs and their activation Z-score are listed in Table 3. These functions fell in the main three groups: those related to synaptic plasticity, such as long-term potentiation or synaptic depression, those related to dendritic morphology, such as neurite formation or neuritogenesis and branching of neurites, and those related to cell fate such as cell proliferation, differentiation and migration. Overall, functions related to cellular and nervous system development dominated the list.

Table 3

Predicted functions associated with the differentially expressed genes in the PFC dataset

Categories	Functions annotation	P-value	Activation Z-score	No. of genes
Cell-to-cell signalling and interaction	Long-term potentiation of brain	0.005	−1.177	9
	Long-term potentiation of cerebral cortex	0.007	−0.923	8
	Neurotransmission	0.009	−0.164	8
	Long-term potentiation	0.013	−0.978	10
	Release of neurotransmitter	0.016	−1.718	4
	Synaptic depression	0.031	−0.527	6
Cell morphology, cellular assembly and organization, cellular development, cellular function and maintenance, cellular growtd and proliferation, embryonic development, nervous system development and function, tissue development	Branching of neurites	0.004	−0.250	13
	Neuritogenesis	0.004	−0.250	15
	Dendritic growtd/branching	0.008	−1.342	11
	Outgrowtd of neurites	0.050	1.898	6
Cellular assembly and organization, cellular development, cellular growtd and proliferation, nervous system development and function, tissue development	Microtubule dynamics	0.003	−0.346	17
	Growtd of neurites	0.039	2.129	8
Cellular development, cellular growtd and proliferation, nervous system development and function, tissue development	Proliferation of neuronal cells	0.013	0.995	11
	Development of neurons	0.001	0.017	20
Cellular growtd and proliferation	Proliferation of cells	0.002	0.512	27
	Generation of cells	0.000	0.101	22
Cellular development	Differentiation of cells	0.001	0.460	19
Lipid metabolism, molecular transport, small molecule biochemistry	Concentration of lipid	0.007	−2.400	7
Nervous system development and function, tissue morphology	Density of neurons	0.014	−2.213	7
Developmental disorder, neurological disease, organismal injury and abnormalities	Cerebral dysgenesis	0.024	0.931	4
Behaviour	Behaviour	0.037	0.687	8
Cellular movement, nervous system development and function	Migration of neurons	0.050	−0.124	6

Gene functions predicted by IPA software and selected based on their relevance for the brain, number of genes related to each function and the ability of the software to make a prediction regarding the direction of change (activation Z-score).

Discussion

In this study, adult male BALB/c mice were exposed to 6 weeks of UCMS conditions, which resulted in low grooming, hyperactivity, transient anhedonia and change in food reward behaviour in UCMS-exposed mice. Behavioural changes were accompanied by increase in CRP and decrease in Lep plasma levels, reactive microglia changes in the PFC area and decline in the number of neuroblasts and their aberrant migration in the GZ of the hippocampal DG.

UCMS is known to induce a complex behavioural response, with our study being no exception. As expected, UCMS consistently reduced grooming in mice. Hyperactivity was also very prominent, in line with many previous UCMS studies (Couch ; Dournes ), especially where the stressors are administered in the light phase and therefore likely disrupt the circadian rhythm of mice (Aslani ). Hyperactivity following UCMS has been previously linked to impulsivity (Couch ). Contrary to expectations, UCMS-exposed mice demonstrated decreased latency to feed in the NSF test. There are several possible explanations for this effect. First, it could be simply linked to hyperactivity in the UCMS group. Second, it is possible that UCMS stimulated food reward drive in mice. Such explanation goes in line with a decline in blood Lep in the UCMS group. Lep is a negative feedback inhibitor released by adipose tissue to regulate food consumption through Lep receptors in the hypothalamus (Ahima and Osei, 2004), thus its reduction could lead to a higher drive for food rewards. Such effect of chronic stress on Lep levels has been previously documented (Ge ; Liu ). Importantly, calorie intake and satiety regulating hormones are known to modulate adult neurogenesis (Garza ; Hornsby ; Morgan ). Change in Lep levels might be reflecting a shift in circadian rhythm, which could have affected the time of feeding drive and levels of hormones regulating it, as has been shown in previous studies of circadian rhythm gene knockouts or mice exposed to disrupted light–dark cycles (Kettner ).

Lastly, it has been shown previously that mouse response in the NSF test is dependent on the ambiguity of the previous experience of aversive stimuli and the state of adult hippocampal neurogenesis. Glover demonstrated that neurogenesis-deficient mice, which undergone ambiguous cue fear conditioning training, showed lower latency to feed in the NSF test compared to their neurogenesis-intact counterparts. This effect was attributed to the failure of neurogenesis-deficient mice to generalize the aversive experience to novel environments. Due to the unpredictable nature of UCMS, it could have exerted an effect similar to ambiguous cue and led to analogous deficiency of fear generalization in the NSF arena. While fear conditioning paradigm differs significantly from the UCMS and their underlying mechanisms might not overlap, it is still interesting to consider that anomalous NSF response in neurogenesis-deficient mice has been described previously.

To investigate molecular changes, which underlie this combination of endophenotypes associated with chronic stress exposure, we used a hypothesis-free approach of genome-wide gene expression microarray. As behavioural monitoring showed the strongest stress behavioural response at 4 weeks, this time point was selected for gene expression analysis. A number of differences were observed in behavioural response between the two cohorts, which could be attributed to the difference in the length of UCMS exposure and variable nature of UCMS protocol well described in the literature (Willner, 2005).

Potential expression changes were assessed in areas previously involved in chronic stress response, such as HIP, the hypothalamus and the PFC. Interestingly, genome-wide gene expression changes were most apparent in the PFC among the three regions investigated, with no overlap among the regions. This finding goes in line with previous studies, which compared the effect of chronic mild stress on gene expression across different brain regions. As such, Hervé et al. (2017) found the strongest effect on UCMS in anterior cingulate cortex compared to the DG, with only five overlapping DEGs between the two regions. Similarly, Surget reported less DEGs in the DG upon UCMS compared to amygdala and the cortex, with the cortex response being more broad in terms of cellular functions involved. Interestingly, the RNAseq study conducted by Nollet showed much more overlap between HIP and the PFC upon UCMS exposure. It is possible that the differences in the methodology of this paper (surgery was conducted on all mice prior to UCMS; older age of mice at the start of UCMS; higher sensitivity of the RNAseq compared to microarray approach) could explain variation in results. Nonetheless, a strong gene expression response of the PFC to the UCMS exposure is a common finding among all studies and has been confirmed by our results. The lack of hippocampal response can be explained by the fact that UCMS mostly induces affective, not cognitive symptoms in mice, which would be expected to correlate with hippocampal response.

Many target genes, upstream regulators and signalling pathways involved in the PFC response in our study could be linked to dendritic remodelling and spine atrophy, a well-described effect of chronic stress on the PFC. Based on the IPA pathway and upstream regulator analysis, glutamatergic and calcium signalling, as well as Htt and Bdnf-centred networks stood out as the most significantly involved. Indeed, repeated stress is known to cause suppressed glutamate receptor expression and signalling in the PFC, which is thought to be linked to dendritic atrophy (Yuen ). Many recent studies have explored the antidepressant potential of ketamine in chronic stress animal models of depression, strengthening the glutamatergic theory of depression (Zhu ; Sun, 2016). In addition, disruption of glutamatergic signalling has been previously linked to hyperactivity (Procaccini ). Calcium signalling has also been implicated in dendrite remodelling, via its effect on protein kinase A and protein kinase C. Calcium is a necessary cofactor in the activation of protein kinase A and protein kinase C, and it has been shown that pharmacological suppression of the protein kinase A and protein kinase C activation prevented spine atrophy following chronic stress exposure (Hains , 2015).

The limited number of DEGs in the HIP did not permit carrying out pathway analyses. It is possible that the lack of differential gene expression detection comes from the heterogeneity of the dorsal and ventral hippocampal tissues, which were pulled together in this study. Alternatively, this could be a limitation of a time point of tissue collection, as the mice displaying neurogenesis changes were exposed to the UCMS for 2 weeks longer than the mice from which tissue for gene expression analysis was derived. However, it is also conceivable that the activation of the signalling pathways affecting dendritic morphology in the PFC has a downstream indirect effect on hippocampal neurons via PFC–hippocampal connections. Indeed, PFC and the HIP are known to be connected via multiple direct and indirect pathways, mostly involved in memory functions (Eichenbaum, 2017). This notion is supported by the enrichment of cell proliferation, differentiation and migration functions in the functional analysis of the PFC dataset. Indeed UCMS reduced dendritic length in PFC as well as in the HIP in a previous study (Morais ). Importantly, in line with our findings, this study also showed that dendritic tree alterations were associated with a more pronounced decrease in neuroplasticity-related gene expression in the PFC compared to the HIP. Interestingly, PFC–HIP connections have been previously implicated in food reward-related behaviour (Hsu ). In addition, increased functional connectivity in the PFC areas and the HIP among other limbic structures was associated with sleep disruption in people reporting depressive problems (Cheng ). Dendritic tree and dendritic spine density measurements in the PFC and the HIP could show if these molecular changes led to dendritic morphology alteration in adult neurons in these two regions.

Similarly, microglial cell density was increased in the PFC but not the HIP of the UCMS-exposed mice. While appearance of reactive microglia upon chronic stress has been described to take place simultaneously in the HIP and the PFC (Wohleb ), previous studies also suggested that the PFC is more sensitive to chronic stress than the HIP (Yuen ). At the same time PFC has been shown to be less prone to microglial dystrophy following chronic stress, which could also explain the presence of microglial increase at the time of sacrifice only in this area (Kreisel ). Interestingly however, inflammatory pathways did not come up in the gene expression analysis. This could be related to the length of UCMS exposure prior to frozen tissue harvesting, as most previous studies collected tissue after 7–9 weeks of UCMS schedule as opposed to 4 weeks in our study.

Finally, our genome-wide gene expression analyses suggested the involvement of two mechanisms not predicted from the behavioural or immunohistochemical experiments carried out in this project, namely the modulation of myelination by Myrf and epigenetic pathways regulated by Mecp2. Future myelin measurements and epigenetic regulation assays could show if indeed UCMS caused any functional changes in these domains.

To conclude, UCMS induces major gene expression changes in the PFC, which could potentially underlie UCMS-associated deficiency in the adult hippocampal neurogenesis. Glutamatergic and calcium signalling, Lep signalling, as well as Htt and Bdnf-regulated networks were identified as the main pathways involved. These gene expression changes preceded changes in food reward behaviour, increase in microglia density in the PFC and a decline in the number of neuroblasts in the hippocampal DG, as well as their aberrant migration through the GZ. These findings highlight the heterogeneous neurobiological effects of the chronic stress exposure, which could lead to depression-like phenotype.

Supplementary material

Supplementary material is available at Brain Communications online.

Click here for additional data file.