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.


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� on-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� e et al., 2017;Laine et al., 2018;Wang et al., 2019;Weger et al., 2020). One of the benefits of these RNAstudies 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 metaanalyses, 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.
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� andez-Beltr� an 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).
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)).

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 b-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  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. Dataset from mice subjected to 10 day protocol was used.

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 SRPpathway 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 b-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 (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.
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 nonneural tissue (Carvalho et al., 2003;Davis, 1994;Sasson et al., 2004). b-catenin is a component of the Wnt-signaling pathway which when activated facilitates translocation of b-catenin into the nucleus where it mediates transcription of target genes. Postmortem studies have reported dysregulation in PFC b-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, b-catenin is implicated in stress resilience (Dias et al., 2014). b-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 celltype 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 highresolution 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. 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.

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.