The impact of chronic stress on the PFC transcriptome: a bioinformatic meta-analysis of publicly available RNA-sequencing datasets

Abstract

The prefrontal cortex (PFC) is one of several brain structures that are sensitive to chronic stress exposure. There have been several studies which have examined the effects on chronic stress, using various protocols such as chronic unpredictable stress and chronic social defeat stress, on the PFC transcriptome. In this report, a bioinformatic meta-analysis of publicly available RNA sequencing datasets (101 samples) from seven chronic stress studies was carried out to identify core PFC transcriptional signatures that underpin behavioral phenotypes including resilience and susceptibility. The results showed 160 differentially expressed genes in chronic stress mice compared to controls with significant enrichment in mechanisms associated with translation and localization of membrane-bound proteins with a putative effect on synaptic plasticity in glutamatergic neurons. Moreover, the meta-analysis revealed no differentially expressed genes in resilient mice but 144 in susceptible mice compared to controls, of which 44 were not identified in the individual studies. Enrichment analysis revealed that susceptibility genes were most affected in oligodendrocytes and linked to mechanisms which mediate biochemical, bidirectional communication between this cell-type and myelinated axons. These results provide new avenues for further research into the neurobiology and treatment of chronic stress-induced disorders.

Keywords:

• Stress
• resilience
• susceptibility
• prefrontal cortex
• transcriptome
• myelination

Introduction

The prefrontal cortex (PFC) is one of several limbic structures which directs the autonomic and neuroendocrine response to stress exposure in addition to its array of cognitive and behavioral functions (Arnsten, 2009; Ulrich-Lai & Herman, 2009). Chronic stress exposure has several reported effects on the PFC structure, morphology, and connectivity with other structures (for review, see McKlveen et al., 2019; Negrón-Oyarzo et al., 2016); this includes a reduction in the dendritic complexity of pyramidal neurons (Goldwater et al., 2009), alterations in synaptic architecture (Shepard & Coutellier, 2018), an imbalance in (glutamatergic) excitatory vs (interneuron) inhibitory neurotransmission in favor of the latter (Shepard & Coutellier, 2018), and a decoupling of connections between the PFC and its output structures such as the amygdala (Jovanovic et al., 2011). These findings have led to the theory that chronic stress exposure takes the PFC “offline,” removing the “top-down” control and influence over functions, and other brain structures that it would otherwise have which are critical for return to homeostasis (Datta & Arnsten, 2019).

Complementing these findings, changes in the PFC transcriptome detected using RNA-sequencing strategies have revealed new insights and mechanisms which potentially underpin many of the above stress-induced phenotypic changes (Bagot et al., 2016; Bondar et al., 2018; Cheng et al., 2018; Labonté et al., 2017; Laine et al., 2018; Wang et al., 2019; Weger et al., 2020). One of the benefits of these RNA-studies is that authors have deposited the sequencing datasets into the NCBI Gene Expression Omnibus (GEO) database (Barrett et al., 2013). This is free to access and researchers can download and re-analyze the datasets to carry out meta-analyses, accounting for the variability within studies and study-specific effects (Rau et al., 2014; Rung & Brazma, 2013).

Accordingly, in this study the datasets from the above studies were downloaded and analyzed using the same bioinformatic pipeline to identify PFC transcriptional signatures with high signal-to-noise ratio which underpin the behavioral, physiological, and endocrinological responses to chronic stress, including resilience and susceptibility. Moreover, functional enrichment analyses were carried out to highlight potential mechanisms and cell-types implicated in the etiology of stress-induced psychopathologies.

Methods

The NCBI Geo Database was queried with the search term “chronic AND mouse AND stress” with the additional filter, “Expression profiling by high throughput sequencing” (Barrett et al., 2013). Projects were selected for analysis which met the following criteria adapted from a similar recent comparative analysis of studies using both acute and chronic stress protocols (Flati et al., 2020): (i) had an accompanying publication with metadata which could be referenced, (ii) at least six samples were sequenced with three biological replicates per group, (iii) exclude studies/samples in which mice received injections or were transgenics, (iv) exclude studies which focus on early life stress, (v) whole transcriptome sequencing. This resulted in data from 101 samples from seven studies (Table 1) which were imported into a Galaxy environment.

Table 1. List and details of studies used in the meta-analysis.

Bioproject	Reference	Study title	Stress paradigm	Housing	Duration	Age	Substrain	Sex	Phenotype	Time of sacrifice	Number of samples	Sequencing
PRJNA293822^a	Bagot et al. (2016)	Circuit-wide transcriptional profiling reveals brain region-specific gene networks regulating depression susceptibility.	Chronic social defeat	Single	10 days	8 weeks	C57BL/6J	Male	Control, resilient, susceptible	2 days, 28 days after last defeat	9	Single-end, 100 bp
PRJNA398031^b	Labonté et al. (2017)	Sex-specific transcriptional signatures inhuman depression.	Chronic variable stress	Grouped	3 weeks	8 weeks	C57BL/6J	Male/Female	Control, stress	Day after last stressor	19	Paired-end, 50 bp
PRJNA323485^c	Bondar et al. (2018)	Molecular adaptations to social defeat stress and induced depression in mice.	Chronic social defeat	Single	10 days, 30 days	10–12 weeks	C57BL/6J	Male	Control, stress	2 days after last defeat	13	Paired-end, 100 bp
PRJNA430409	Laine et al. (2018)	Genetic control of myelin plasticity after chronic psychosocial stress.	Chronic social defeat	Single	10 days	7 weeks	C57BL/6NCrl	Male	Control, resilient, susceptible	6–8 days after last defeat	18	Single-end, 96 bp
PRJNA415948	Cheng et al. (2018)	10–11 translocation proteins modulate the response to environmental stress in mice.	Chronic restraint stress	Grouped	2 h/day for 2 weeks	7–8 weeks	C57BL/6J	Male	Control, stress	After last restraint	6	Single-end, 100 bp
PRJNA525440	Wang et al. (2019)	RNA-sequencing and bioinformatics analysis of long noncoding RNAs and mRNAs in the prefrontal cortex of mice following repeated social defeat stress.	Chronic social defeat	Grouped	6 days	6–8 weeks	C57BL/6NCrl	Male	Control, stress	After last behavioral test	6	Paired-end, 100 bp
PRJNA625244	Weger et al. (2020)	Mitochondrial gene signature in the prefrontal cortex for differential susceptibility to chronic stress.	Multimodal chronic restraint	Grouped	3 weeks	6 weeks	C57BL/6J	Male	Control, stress	2 days after behavioral testing	30	Single-end, 70 bp

^a Dataset from mice culled 2 days after last defeat was used.

^b Sequencing carried out on samples from non-behavioral cohort. Data from males were analyzed.

^c Dataset from mice subjected to 10 day protocol was used.

Reads were aligned to the reference Mus musculus genome (GRCm38/mm10) using HISAT2, annotated (Gencode v.M25) and counted to gene exons using featureCounts. Default settings were used for all analyses. Biostatistics and visualization were run in R (version 3.6.3) with the Rstudio GUI (version 1.2.5033). Filtering was done to include only protein-coding genes with at least 10 counts in half the samples of each study.

Differential gene expression of individual datasets was analyzed using the DESeq2 package with default parameters (Love et al., 2014). For the meta-analysis of datasets, two approaches were used. The first was by combining sample data from all studies. There were several lab-specific experimental variables for each study (Table 1) which ranged from housing type, age of the mouse to the sequencing strategy. As such a new variable known as “lab,” was used and DESeq2 analysis was adjusted to control for the “lab effect” (design=∼lab + group). The second approach was using the metaRNASeq package to perform inverse normal p value combination tests as previously described (Fernández-Beltrán et al., 2021; Rau et al., 2014; Ziff et al., 2022).

Overlap analyses between datasets were carried out and visualized using the UpSetR package (Conway et al., 2017).

Enrichment analysis was carried out using EnrichR on significant differentially expressed genes (Chen et al., 2013) using the following libraries: GO_Biological Process_2021, GO_Molecular Function_2021, GO_Cellular Component_2021, KEGG_2019_Mouse, DSigDB, and Allen_Brain_Atlas_10x_scRNA_2021. Cluster identification and cell-type information were obtained from the Allen Institute transcriptomic atlas (https://taxonomy.shinyapps.io/ctx_hip_browser/). GSEA analysis was performed as previous described (Mootha et al., 2003; Subramanian et al., 2005) using the ontology gene set database (ms_gsea_c5 = msigdb(species=“Mus musculus”, category=“C5”).

Identification of predicted upstream regulators of differentially expressed genes was done using QIAGEN IPA (QIAGEN Inc., https://digitalinsights.qiagen.com/IPA).

The top five significant results of each analysis have been presented in bar plots. Unless otherwise stated, significance for all analyses was set at FDR <0.05.

Results

Because different sequencing platforms were used for the analysis in each of the studies, there was a need first correlate measures of average gene expression (normalized counts) between datasets. The higher the correlation, the greater the chance of finding similarities between datasets. Across all seven studies, there were 8767 genes that were commonly identified and there were significant Pearson correlations between all studies suggesting expression profiles are largely comparable (Figure 1(A)).

Figure 1. Analysis of RNA-seq datasets from studies which have used chronic stress protocols in male mice to examine the impact on the PFC transcriptome. (A) Pearson correlation plot of average gene expression across all the seven studies that were the included in this study. (B) Heatmaps indicating degree of fold-change and the number of differentially expressed genes in each study and in the combined analysis (bottom row) for stress vs control mice. (C) UpSet plot of overlaps between the various datasets that were analyzed. The vertical lines which connect black circles indicate intersecting sets of genes between studies. The y-axis of the histogram indicates the number of intersecting genes and the bar plot on the right-hand side indicates the number of genes in each set. (D–G) EnrichR analysis of differentially expressed genes in the combined analysis for stress vs control mice. (H) Predicted upstream regulators of the differentially expressed genes in the combined analysis as determined using IPA.

Stress vs control

The number of genes that were differentially expressed in stressed mice compared to controls in each of the seven studies varied from none to 969 (Figure 1(B), Supplementary Tables 1(A–F)) and this is reflective of the impact of methodological differences on differential gene expression in the PFC. In combining all the datasets there were 160 genes that were differentially expressed (Supplementary Table 1G). Additionally, a total of 259 differentially expressed genes were found to overlap across six of the seven studies and the combined analysis (Figure 1(C)). MetaRNAseq analysis yielded no significant findings. To avoid bias for any specific study, EnrichR analysis was carried out on the combined analysis (160 genes) and yielded significant findings across multiple libraries. Notably, there was significant enrichment in functional mechanisms in stressed mice associated with transport of proteins to target membranes (Figure 1(D)) and ribosomal function (Figure 1(E,F)). There was also significant enrichment for genes that were highly expressed in specific glutamatergic neurons (Mouse 355 V3d up). IPA analysis identified several putative upstream regulators of differentially expressed genes in stressed mice including luteinizing hormone (Lh), and CTNNB1 also known as β-catenin (Figure 1(E)). GSEA analysis revealed no significant findings.

Resilient/susceptible vs control

There were 216 differentially expressed genes in resilient mice but only from one study (Bagot et al., 2016). In contrast, there were a 388 differentially expressed genes in susceptible mice across both studies (Figure 2(A), Supplementary Tables 2(A–C)). Analysis of combined datasets revealed no differentially expressed genes in resilient mice compared to controls but 92 differentially expressed genes in susceptible mice (Figure 2(A), Supplementary Table 2D). MetaRNAseq analysis of all the datasets yielded 144 differentially expressed genes in susceptible mice of which 44 were not identified in the individual studies alone (Supplementary Tables 2(E,F)). Additionally, a total of 310 differentially expressed genes were found to overlap across the two studies, the combined and metaRNAseq analyses in the susceptible mice (Figure 2(B)). In particular, there were 76 differentially expressed genes that were common to the combined and metaRNAseq analyses (Supplementary Table 2G). As above, to avoid bias for any specific study, EnrichR analysis was carried out on this set of 76 genes revealed significant effects on channel activity (Figure 2(C)), gap junction (Figure 2(D)), and oligodendrocytes were the most affected cell-type (Figure 2(E)). IPA analysis of this gene set revealed the transcription factor Sox2 as one of the upstream regulators. GSEA analysis revealed no significant enrichment.

Figure 2. Analysis of RNA-seq datasets from studies which have used chronic social defeat stress to examine the impact on the PFC transcriptome in male mice classified as resilient or susceptible. (A) Heatmaps indicating degree of fold-change and the number of differentially expressed genes in each study and in the combined analysis (bottom row) for susceptible vs control and resilient vs control mice. (B) UpSet plot of overlaps between the various datasets that were analyzed. (C–E) EnrichR analysis for differentially expressed genes that were common in the combined analysis and metaRNAseq analysis for susceptible vs control mice. (F) Predicted upstream regulators of the differentially expressed genes that were common in the combined analysis and metaRNAseq analysis for susceptible vs control mice.

Discussion

The PFC is one of the several brain structures sensitive to stress exposure which can, if chronic or uncontrollable, lead to it being taken “off-line,” thus compromising its important top-down role in controlling cognition, mood, and behavior (Datta & Arnsten, 2019; Woo et al., 2021). Here, I have presented the results of a bioinformatic meta-analysis of several studies that have examined the impact of chronic stress on the PFC to identify common transcriptional signatures which underpin stress-induced phenotypes.

Using the same bioinformatic pipeline on each dataset revealed wide ranging effects on chronic stress on the PFC transcriptome. However, analysis of the combined datasets yielded 160 differentially expressed genes which suggests their general involvement in the response to chronic stress, regardless of paradigm specifics. Enrichment analysis of these genes revealed significant involvement in several processes including signal recognition particle (SRP)-co-translational targeting of proteins to membranes as well as ribosomal function. The SRP-pathway is initiated by the binding of SRP to a specific amino acid sequence on a nascent protein emerging from a translating ribosome. Translation is briefly halted while this complex is delivered to the plasma membrane via a cognate SRP receptor where it unloads the cargo. Translation resumes, SRP and SR dissociate and move on to other targets (Egea et al., 2005; Nagai et al., 2003; Saraogi & Shan, 2011). Combined with the findings that specific glutamatergic neurons are more affected than other cell types provide supporting evidence that chronic stress in the PFC disrupts excitatory neurotransmission and synaptic plasticity, going “off-line,” in part due to the dysregulated expression of membrane-bound proteins (Arnsten, 2009; Moghaddam, 2002).

Selected top predicted upstream regulators for this set of differentially expressed genes in chronic stress mice included Lh and β-catenin. Lh is released in a pulsatile manner by gonadotrophs in the anterior pituitary downstream of an activated hypothalamic–pituitary–gonadal axis. This is a neuroendocrine network which functions to control reproductive function but is also under the control of the hypothalamic–pituitary–adrenal axis. Indeed, evidence from cross-species experiments have reported an inverse relationship between stress-induced increase in corticosterone levels and Lh levels (Son et al., 2022). Notably, work in rats has shown that stress-induced effects on Lh levels are sensitive to glucocorticoid receptor antagonists (Briski et al., 1995), antibodies to corticotropin releasing factor (Nemoto et al., 2010), serotonin synthesis inhibitors, and serotonin receptor antagonists (Armario et al., 1993). Lh binds to Lh-receptors and activates a range of intracellular signaling cascades which may culminate in transcriptional responses, though this evidence mainly comes from work using non-neural tissue (Carvalho et al., 2003; Davis, 1994; Sasson et al., 2004). β-catenin is a component of the Wnt-signaling pathway which when activated facilitates translocation of β-catenin into the nucleus where it mediates transcription of target genes. Postmortem studies have reported dysregulation in PFC β-catenin expression and this has also been observed in animal models of aspects of depression (Teo et al., 2018). It is worth noting here that in the nucleus accumbens, β-catenin is implicated in stress resilience (Dias et al., 2014). β-catenin in the PFC has also been implicated in the mechanism of action of antidepressants (Chen et al., 2012; Mohamed et al., 2020).

Overall, susceptibility to chronic stress was associated with greater degree of transcriptional change than resilience. This transcriptional change had significant effects at the bidirectional, biochemical communication interface between oligodendrocytes and the axons which they myelinate as indicated by the enrichment analysis. Oligodendrocytes express glutamatergic and GABAergic receptors which influence their ability to facilitate axonal conduction and they, in turn, are also susceptible to depolarization by passing action potentials in the axons (Edgar & Sibille, 2012; Mazuir et al., 2021). Recent studies have reported the differential impact of chronic stress on aspects of PFC myelination in susceptible and resilient mice, with the former group being more affected (Antontseva et al., 2020; Bonnefil et al., 2019; Laine et al., 2018). The findings reported herein suggest that the susceptibility phenotype may be, in part, due a breakdown in oligodendrocyte-axonal information transfer in the PFC and consequently lead to its disengagement and decoupling from the neural networks over which it exerts control (Edgar & Sibille, 2012).

Sox2 was identified as a predicted upstream regulator of gene expression changes in susceptible mice. It is a conserved transcription factor with a significant role in neurodevelopment but also expressed in adulthood in a region- and cell-type specific manner; in astrocytes, deletion of the Sox2 gene induces resistance to the effects of traumatic brain injury in mice (Chen et al., 2019) while deletion in GABAergic neurons in the suprachiasmatic nucleus, the brain’s “internal clock,” results in increased anxiety-like behaviors and sex-specific effects on motivation (Boehler et al., 2021). Its functional role in the PFC and the stress-response remains undetermined. An important caveat with this discussion is that susceptibility was defined as reduced social interaction with a novel conspecific. But this could also be a manifestation of an adaptive, high-vigilance response developed as a protective strategy over the course of the social defeat protocol (Meshalkina & Kalueff, 2016). Examining the PFC transcriptional response in stressed mice with differing reward learning and processing abilities is one approach to be considered in future work to rule out this alternative explanation and more accurately discriminate between phenotypes (Gururajan et al., 2019).

The deposition of transcriptomics into public database provides opportunities for researchers to make new discoveries and generate testable hypotheses at low-cost through integrative meta-analyses using bioinformatic tools. There are pros and cons to using this approach (Feichtinger et al., 2012; Ramasamy et al., 2008). Pros include (i) the ability to increase statistical power and obtain more precise estimates of gene expression by increasing the sample size, (ii) the ability to reveal a set of generalizable, reproducible, core-over-represented genes – a so-called meta-signature – in analogous studies, (iii) and reducing the probability of false negative results. Cons are that (i) studies with larger sample sizes may have a larger influence on the meta-analysis outcome than a smaller ones, (ii) the “overlooking” of any bias in the publications from which data are extracted, and (iii) the inherent assumption that different experimental paradigms make similar demands on the organism. Other limitations of this bioinformatic study are the (i) absence of sufficient RNA-sequencing data from chronic stress studies in females which met the exclusion criteria (Deonaraine et al., 2020; Girgenti et al., 2019; Issler et al., 2020) and that (ii) all input studies using bulk RNA-sequencing techniques which does not reflect the known heterogeneity in PFC cell-type specific gene expression responses in stress-related disorders (Nagy et al., 2020). The inclusion of both male and female animals in the same study and cell-type analyses in the future should address this knowledge deficit to provide high-resolution insights into sex-specific PFC transcriptional signatures. With this proviso, the results lead to the hypothesis that chronic stress disrupts mechanisms that regulate synthesis and localization of membrane-bound proteins involved in neurotransmission and synaptic plasticity in the PFC. Moreover, the findings also indicate the specific sensitivity of oligodendrocytes in the PFC to chronic stress that may lead to stress-susceptible phenotypes in mice. Future validation will be required to rigorously test these hypotheses before translation into better treatments for stress-induced psychopathologies.

Acknowledgements

AG wishes to thank Dr Kandarp Joshi and the Sydney Informatics Hub, a Core Research Facility of the University of Sydney, for their assistance with data analysis and Prof. Avi Ma'ayan with interpretation of the EnrichR brain cell-type analysis.

Disclosure statement

No potential conflict of interest was reported by the author.

Data availability statement

All data generated or analyzed during this study are freely available on the NCBI GEO website. Accession details have been included in Table 1. Scripts for analysis of data are available on request.

Additional information

Funding

AG was supported by a University of Sydney Research Fellowship.

Notes on contributors

Anand Gururajan

Dr Anand Gururajan is a Research Fellow at the Brain & Mind Centre at the University of Sydney. His research is focused on understanding the impact of stress on the brain by using animal models and high-resolution, high-throughput sequencing techniques, with a view to identifying new drug targets for the treatment of stress disorders such as major depression and PTSD.