Stressful life events are major environmental risk factors for anxiety disorders, although not all individuals exposed to stress develop clinical anxiety. The molecular mechanisms underlying the influence of environmental effects on anxiety are largely unknown. To identify biological pathways mediating stress-related anxiety and resilience to it, we used the chronic social defeat stress (CSDS) paradigm in male mice of two inbred strains, C57BL/6NCrl (B6) and DBA/2NCrl (D2), that differ in their susceptibility to stress. Using a multi-omics approach, we identified differential mRNA, miRNA and protein expression changes in the bed nucleus of the stria terminalis (BNST) and blood cells after chronic stress. Integrative gene set enrichment analysis revealed enrichment of mitochondrial-related genes in the BNST and blood of stressed mice. To translate these results to human anxiety, we investigated blood gene expression changes associated with exposure-induced panic attacks. Remarkably, we found reduced expression of mitochondrial-related genes in D2 stress-susceptible mice and in exposure-induced panic attacks in humans, but increased expression of these genes in B6 stress-susceptible mice. Moreover, stress-susceptible vs. stress-resilient B6 mice displayed more mitochondrial cross-sections in the post-synaptic compartment after CSDS. Our findings demonstrate mitochondrial-related alterations in gene expression as an evolutionarily conserved response in stress-related behaviors and validate the use of cross-species approaches in investigating the biological mechanisms underlying anxiety disorders.

Genetic and environmental factors contribute to the etiology of psychiatric diseases but the underlying mechanisms are poorly understood. Chronic psychosocial stress is a well-known risk factor for anxiety disorders. To identify biological pathways involved in psychosocial stress-induced anxiety and resilience to it, we used a well-characterized mouse model of chronic social defeat stress (CSDS) in two inbred mouse strains, C57BL/6NCrl (B6) and DBA/2NCrl (D2), which differ in their susceptibility to stress. We focused on the bed nucleus of the stria terminalis, a key brain region behind stress-response and anxiety, and carried out genome-wide analysis of mRNA, and miRNA expression, and protein abundance. Bioinformatic integration of these data supported differences in mitochondrial pathways as a major stress response. To translate these findings to human anxiety, we investigated blood cell gene expression in mice and in panic disorder patients exposed to fearful situations and experiencing panic attacks. Concurring with our brain findings, expression of mitochondrial pathways was also affected in mouse and human blood cells, suggesting that the observed stress response mechanisms are evolutionarily conserved. Therefore, chronic stress may critically affect cellular energy metabolism, a finding that may offer new targets for therapeutic interventions of stress-related diseases.

Funding: This work was supported by the ERA-NET NEURON grant ( https://www.neuron-eranet.eu/ ; AnxBio; I.H., A.E., C.W.T., A.C.), the University of Helsinki (I.H.), Sigrid Jusélius Foundation ( https://sigridjuselius.fi/ ; I.H.), Doctoral Program Brain and Mind, University of Helsinki (M.L. and Z.M.), the Orion Research Foundation ( https://www.orion.fi/en/rd/orion-research-foundation/ , Z.M.) and the Max Planck Society (Z.M., L.R., G.M., C.R., B.N, D.I.P, C.W.T.). The EM-unit is supported by the Institute of Biotechnology (University of Helsinki), Helsinki Institute of Life Science and Biocenter Finland (E.J.). A.C. is the head of the Max Planck Society – Weizmann Institute of Science Laboratory for Experimental Neuropsychiatry and Behavioral Neurogenetics. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

To identify the core differentially expressed molecules and pathways underlying pathological anxiety and resilience to it, we employed a multi-omics approach in mice. We applied the CSDS model to induce anxiety-related phenotypes and to identify molecular markers for susceptibility and resilience in male mice. We used two inbred strains, the largely stress-resilient C57BL/6NCrl (B6) and stress-susceptible DBA/2NCrl (D2), due to strong genetic background effects in the mouse stress response [ 4 , 20 ]. We investigated gene (RNA-seq), microRNA (miRNA-seq), and protein (LC-MS/MS) expression differences in the BNST between the stress-resilient, stress-susceptible, and control mice. As a translational effort, we also carried out RNA-seq from mouse post-CSDS blood, and longitudinal microarray-based gene expression profiling from blood cells of panic disorder patients who experienced high degrees of stress and anxiety due to exposure to phobic situations. To identify converging anxiety-related gene sets and pathways, we conducted pathway and gene set enrichment analyses of the data sets from mouse BNST samples, and mouse and human blood cells. Altogether, our results indicate anxiety induces a genetically-controlled evolutionarily conserved response in mitochondrial pathways.

Omics-based approaches require the investigation of etiologically relevant tissue. The CSDS model has been used to study transcriptomic effects of chronic psychosocial stress in brain regions classically associated with anxiety and depression-related behaviors, including the medial prefrontal cortex, hippocampus, amygdala, and nucleus accumbens [ 4 , 16 ]. In addition, converging lines of anatomical, physiological, and mechanistic evidence suggest that the bed nucleus of the stria terminalis (BNST) functions as the center of integration for limbic information [ 17 ], monitoring the genesis of long-term responses to anxiety [ 18 ], such as anticipatory anxiety. The BNST is an especially relevant brain region for stress-related anxiety. It has been suggested that it plays a role in valence surveillance by processing salient information on physical and social contexts, collected through its numerous projections throughout the brain [ 17 ]. However, only few anxiety studies have generated large-scale genomics data from the BNST [ 19 ].

The development of much-needed novel targets for therapeutic intervention of anxiety disorders is limited by the ignorance of the molecular and cellular mechanisms associated with events that initiate and maintain pathological anxiety. The phenotypic heterogeneity of human populations and the high variability of environmental influences [ 13 ], along with a limited access to brain tissue samples, make it difficult to identify the biological basis of anxiety disorders. These challenges can, to some extent, be controlled in animal models. Cross-species approaches are therefore expected to reveal specific biological mechanisms underlying anxiety disorders [ 14 ]. They can be especially powerful for multi-omics studies allowing hypothesis-free identification of the most significant biological pathways associated with specific exposures, while at the same time probing the effects of different genetic backgrounds on the outcomes. Therefore, the integration of multiple levels of information [ 15 ] and the translation of the results from animal models to the human condition are critical for the success of cross-species approaches.

Anxiety disorders are common stress-associated psychiatric disorders [ 7 , 8 ]. They are characterized by an excessive physiological and emotional response in the absence of real threat or imminent danger. Among the most debilitating anxiety disorders is panic disorder, which involves sudden recurrent surges of intense fear and discomfort called panic attacks [ 9 ]. However, panic attacks are not exclusive to panic disorder, but also frequent in other anxiety disorders. They are typically associated with severe perceived physical and mental stress, feeling of loss of control and avoidance behavior. Epidemiological studies show a major impact of both cumulative and specific life events or stressors, such as threat or psychosocial interpersonal life events, on the development of all anxiety disorders [ 10 , 11 ], including panic disorder [ 12 ].

Chronic stress is a significant risk factor for human anxiety disorders [ 1 ], yet not all individuals exposed to stress develop a clinically relevant anxiety symptomatology. The underlying reasons for these differences are not yet fully understood but involve an interaction of complex genetic and environmental factors that vary among individuals resulting in stress susceptibility or resilience. The chronic social defeat stress (CSDS) model is a well-established mouse paradigm of psychosocial stress, with construct, face, discriminative, and predictive validity for stress-related disorders [ 2 – 4 ]. It involves 10 days of brief daily confrontation of two conspecific male mice, a resident-aggressor and an intruder who reacts with defensive, flight, or submissive behavior [ 5 , 6 ]. Although all defeated mice experience stressful stimuli, only some develop stress-related symptoms, measured as social avoidance, making it an excellent model to investigate mechanisms associated with susceptibility and resilience. We have previously shown that genetic factors strongly affect the behavioral outcome of the CSDS, since different inbred mouse strains vary in the proportion of susceptible and resilient animals as well as in their stress coping behaviors [ 4 ].

To determine if the observed differences in mitochondrial gene and protein expression are associated with changes in mitochondrial morphology in the pre-synaptic or post-synaptic compartments of neurons, we carried out transmission electron microscopy (TEM) in the BNST after CSDS in B6 and D2 mice. We classified mitochondria as “synaptic” if the synaptic density and vesicles within the post-synaptic terminal were clearly visible ( Fig 8A ). We observed stress-associated differences in maximum (length) and minimum (width) diameters of mitochondrial cross-sections in the B6 but not in the D2 strain. The B6 susceptible mice had on average 8.4% shorter mitochondrial cross-sections (maximum diameter) than B6 controls (Padj = 0.003) ( Fig 8B ), but no differences in the width ( Fig 8C ). However, the mean mitochondrial cross section length/width ratio in D2 susceptible, but not B6, mice was 5% larger than in the resilient group (Padj = 0.003) indicative of increased maximum, but decreased minimum, diameter ( Fig 8D ). The mean number of mitochondrial cross-sections was not influenced by stress ( Fig 8E ). However, when we investigated the pre-synaptic and post-synaptic compartments separately ( Fig 8F and 8G ), we detected 39% more pre-synaptic cross-sections in susceptible compared to the control B6 mice and 46% less post-synaptic cross-sections in resilient compared to the susceptible B6 mice. In addition, we observed significant strain differences. In general, B6 mice had wider diameter and smaller number of mitochondrial cross-sections than D2 mice, both pre- and post-synaptically ( Fig 8B, 8F and 8G ). Thus, consistently with our gene and protein differential expression data, we observed significant strain-dependent changes in mitochondrial morphology in the BNST following CSDS.

We next integrated results from the GSEA performed with BNST and blood cell transcriptomic data from stressed mice (data sets A and F) and panic disorder patients exposed to panic attacks (data set H). Due to the large number of data sets, differences in tissue source and organisms, a recommended exploratory false discovery rate (FDR) of 25% [ 43 ] was applied. Again, we found enrichment of several mitochondria and translational control-related gene sets ( Fig 7D ). Intriguingly, the enriched gene sets in the panic disorder patients’ blood cells, both directly and 24 h after exposure-induced panic attack in comparison to the baseline, showed a similar pattern to D2 mice exposed to CSDS, including downregulation of the genes in the enriched pathways. Thus, although we found opposite gene expression patterns in the two mouse strains, the pattern of the highly stress-susceptible mouse strain resembled that of panic disorder patients. These results suggest that the systemic transcriptional regulation of mitochondrial pathways is an evolutionarily-conserved response to anxiety in both humans and mice.

We next examined the number of DE probes (P < 0.05) in panic disorder patients’ blood cells (data set H), collected directly and 24 h after exposure-induced panic attack in comparison to baseline measurement. Of all 47291 analyzed probes, 5791 (corresponding to 4149 genes) were nominally significantly DE in at least one time point. We observed a higher number of DE genes immediately after the exposure-induced panic attack (n = 2753) than 24 h after (n = 2045) ( Fig 7C ). A large number of the identified DE genes (n = 649) was common between the comparisons. Following differential expression analysis, the probes were annotated to HGNC Gene Symbols and DE genes (P < 0.05) and analyzed for overrepresentation of GO terms. The most significant GO term of Biological process shared between both time points was “Translation” (P FDR = 1.72E -7 at 0 h and P FDR = 1.97E -3 at 24 h). The fifth most significant shared term was “Oxidative phosphorylation” (P FDR = 8.24E -5 at 0 h and P FDR = 1.39E -3 at 24 h) ( S10 Table ). Overall, we found several enriched GO terms (P FDR < 0.05) related to translational control and mitochondria in both time points in comparison to the baseline.

To determine which miRNAs are differentially expressed in blood cells in response to CSDS, we performed miRNA-seq followed by DE analysis ( Fig 7B ). Of the 64 DE miRNAs, only three were differentially expressed in both B6 and D2 mice: miR-148a, miR-181b, and miR-592 ( S9 Table ). While miR-592 was expressed at a nominally lower level in both B6 and D2 stress-susceptible mice in comparison to the same-strain controls (FC = -1.71, P = 0.021, Padj = 0.797), miR-148a and miR-181b were expressed at nominally higher levels in B6 susceptible (FC = 1.34, P = 0.012, Padj = 0.731 and FC = 1.66, P = 0.004, Padj = 0.526, respectively), but lower levels in D2 susceptible mice (FC = -1.26, P = 0.042, Padj = 0.328 and FC = -1.41, P = 0.038, Padj = 0.325, respectively), compared to controls. Additionally, miR-181b was expressed at nominally higher levels in B6 susceptible compared to resilient mice (FC = 1.57, P = 0.008, Padj = 0.992). Two of the DE miRNAs, miR-34c and miR-3076, were found in the same comparison (B6 resilient vs. control) in both BNST (D) and blood (G) data sets. Notably, miR-34c was expressed at a nominally lower level in B6 resilient mice in comparison to the controls in both tissues (FC = -2.49, P = 0.040, Padj = 0.955 and FC = -1.77, P = 0.041, Padj = 0.187, in blood and the BNST, respectively). Of these three miRNAs, miR-148a has been associated with panic disorder [ 45 ], and miR-181b proposed as a non-invasive diagnostic marker for schizophrenia [ 46 ].

(A-C) Venn diagrams showing overlap between significantly differentially expressed (A) mouse genes (data set F; P < 0.05 and |FC| ≥ 1.2), (B) mouse miRNAs (data set G; P < 0.05 and |FC| ≥ 1.2), and (C) genes (data set H; P < 0.05) in panic disorder (PD) patients’ blood cells. (D) Merged heatmap from gene set enrichment analysis showing the significantly enriched gene sets (P FDR < 0.25) present in at least 50% of all comparisons. Gene sets which did not pass the cut-off are marked in gray (NA). A positive (or negative) NES for a given gene set implicates its overrepresentation at the top (or bottom, respectively) of the ranked list of upregulated (or downregulated, respectively) genes. Gene sets are ordered alphabetically. B6: C57BL/6NCrl; BNST: bed nucleus of the stria terminalis; D2: DBA/2NCrl.

The blood cell differential gene expression analysis in the B6 and D2 stress-susceptible mice compared to the same-strain controls identified 568 nominally DE genes in the B6 and 771 nominally DE genes in the D2 strain ( S2 Table ). The majority of them (91.1%, n = 102) were expressed in opposite directions between the strains ( Fig 7A ), similarly to the BNST data (data set A). To ask which biological pathways these DE genes belong to, we carried out hypergeometric statistics using C2 curated gene sets [ 43 , 44 ] and GO term enrichment analysis. We found significant enrichment (P FDR < 0.05) of more than 100 gene sets, with the highest number of common DE genes associated with Alzheimer’s disease gene set (n = 17) ( S7 and S8 Tables). The average expression levels of genes that were nominally differentially expressed in the BNST and expressed in blood cells (n = 788), correlated significantly (Pearson’s correlation, r = 0.410, P = 2.93E -33 ), confirming a highly similar transcriptome response in both the BNST and the blood cells after CSDS.

We next asked if the same pathways and gene sets that we found to be stress-responsive in the mouse BNST, could be identified from an accessible tissue, i.e. blood cells, and whether these pathways and gene sets were also DE in blood of humans with anxiety. One week after CSDS, we collected blood samples from stress-susceptible, resilient, and control B6 and D2 mice, and carried out RNA-seq and miRNA-seq (data sets F and G, respectively). Blood samples were collected from the same mice that were used for BNST RNA-seq (data set A), except for D2 resilient mice where we did not obtain a sufficient amount of blood. As a translational human anxiety disorder data set, we collected samples from panic disorder patients who underwent an exposure intervention. We selected panic disorder as our translational target to concentrate on a phenotypically homogeneous sample. We collected blood samples at baseline, 1 h after anxiety peak during exposure and 24 h after exposure-induced panic attack, and carried out microarray-based gene expression profiling.

CYCS is a central element of the electron transport chain in mitochondria [ 41 , 42 ]. Western blot analysis showed that CYCS was more abundant in both the B6 susceptible (P = 0.027) and resilient mice (P = 0.015) in comparison to controls ( Fig 6B ). This was in contrast to the transcriptomic data (data set A), where Cycs was expressed at nominally higher levels in susceptible compared to both the control (P = 0.042, Padj = 0.686) and resilient mice (P = 0.043, Padj = 0.593) ( Fig 6F ). We did not observe differences in the CYCS abundance in the D2 strain by Western blot, although its abundance was lower in the resilient compared to control mice in the LC-MS/MS data (data set B) ( Fig 6D ). Overall, the Western blot analysis confirmed the divergent expression of both PPP1R1B and CYCS observed in the transcriptomic and proteomic analyses following CSDS.

PPP1R1B acts as an integrator of dopaminergic and glutamatergic signaling, and elevated levels of its truncated isoform have been observed in schizophrenia, bipolar disorders, major depression, and poor cognitive functioning [ 40 ]. Western blot analysis confirmed the proteomic results of significantly more abundant PPP1R1B levels in B6 resilient mice in comparison to the control (P = 0.026) and susceptible (P = 0.003) groups, while in the D2 strain, it was significantly less abundant in resilient compared to susceptible mice (P = 0.033) ( Fig 6A ). Furthermore, PPP1R1B protein levels detected by Western blot analysis correlated significantly with SI ratios in the B6 strain (r = 0.689, P = 0.006), but not in the D2 strain (r = -0.244, P = 0.422) ( S4 Fig ).

Merged heatmap showing the expression fold change (FC) and P-value of resilient vs. control, susceptible vs. control and susceptible vs. resilient comparisons from the BNST transcriptomic and proteomic analyses (data sets A and B, respectively). Only DE molecules (P < 0.01) in at least one of the 12 comparisons are shown. “Mitochondria-related gene sets and pathways” include “Mitochondrial Dysfunction,” “Oxidative Phosphorylation” and “Sirtuin Signaling Pathway” shown in Fig 4A and all gene sets shown in S3 Fig . All remaining canonical pathways included in Fig 4A are referred to as “Mitochondria-unrelated pathways”. Genes and proteins significantly DE in the same comparison between at least two BNST data sets are referred to as overlapping between–omic data sets and were selected based on results shown in S6 Table . Molecules with P > 0.05, are shown in diagonal stripe pattern, while those with P < 0.05 are shown with solid colors. Molecules not present in a data set are shown in gray. B6: C57BL/6NCrl; C: control; D2: DBA/2NCrl; R: resilient; S: susceptible.

To identify the key stress-related molecules shared between transcriptomic and proteomic data sets, we focused on the most significantly dysregulated transcriptomic or proteomic pathways from the IPA analysis (P FDR < 0.05; Fig 4A ) and GSEA (P FDR < 0.05; S3 Fig ) ( Fig 5 ). We asked if any molecule specific to those pathways and gene sets shows similar pattern of expression in at least one of the comparisons (susceptible vs. control, resilient vs. control, or susceptible vs. resilient) between at least two BNST data sets ( Fig 2C , S6 Table ). We found three common molecules, Cycs, Glud1, and Atp6v1e1, in mitochondria-related gene sets and pathways, and six molecules, Adcy5, Ppp1r1b, Qdpr, Gad2, Atp2b1 and Ppp3ca, in other canonical mitochondria-unrelated pathways ( Fig 5 ). We selected two of these proteins, PPP1R1B and CYCS (cytochrome c somatic), which have been previously associated with psychiatric disorders, including anxiety disorders [ 36 – 39 ], for validation by Western blot analysis ( Fig 6 ).

To identify active DE miRNAs and their mRNA targets in the BNST of B6 mice following CSDS, we carried out AGO-RIP-seq (data sets D and E; Fig 2C ) and employed the IPA microRNA Target Filter tool [ 27 ]. It uses experimentally validated interactions from TarBase [ 33 ], miRecords [ 34 ] and micro-RNA related findings from the published literature [ 27 ], and microRNA-mRNA interactions predicted by TargetScan [ 35 ]. We observed 72 (resilient vs. control), 59 (susceptible vs. controls), and 36 (susceptible vs. resilient) DE miRNAs, which were either experimentally validated or computationally predicted (high confidence) to repress from one to multiple co-immunoprecipitated AGO2-bound DE mRNAs ( S5 Table ). Only matching pairs of miRNAs and mRNAs that were DE in the same direction were included in the analysis, since they may form functional AGO2-bound miRNA-mRNA pairs. Next, to further increase the confidence of our findings, we selected miRNA-mRNA interactions identified by more than two IPA microRNA Target Filter sources in at least one comparison. We found nominally lower levels of miR-34c and higher levels of miR-15b in both resilient and susceptible mice in comparison to the controls ( Table 1 ). Additionally, we observed nominally lower levels of miR-99b/100 in the stress-susceptible mice compared to both control (FC = -1.50, P = 6.69E -5 , Padj = 0.024) and resilient (FC = -1.22, P = 0.011, Padj = 0.183) mice, as well as in the resilient mice in comparison to the controls (FC = -1.22, P = 0.005, Padj = 0.076) ( S5 Table and Table 1 ). The DE genes Epdr1, Ntrk2, Cldn11, and Ppp3ca were predicted targets of miR-99b/100 ( S5 Table ). Of these genes, Cldn11 was also nominally differentially expressed in the BNST susceptible vs. control transcriptomic analysis in both strains (data set A; FC = 1.61, P = 0.015, Padj = 0.630 and FC = -1.68, P = 0.009, Padj = 0.608), confirming it as a robust stress-affected gene.

The target genes within the beta-estradiol cluster were predicted to be significantly downregulated (P FDR = 0.015) in the D2 resilient mice in comparison to the controls on the transcriptome level. On the proteome level, the same genes within the beta-estradiol cluster were significantly downregulated in the B6 resilient mice in comparison to the B6 controls (P FDR = 0.015) and significantly inhibited in the D2 susceptible mice in comparison to the resilient group (Z-score = -2.40, P FDR = 0.008). Thus, as previously shown [ 32 ], some stress-associated DE genes may be under the regulation of hormones.

To identify transcriptional regulators, which could explain the observed changes in gene and protein expression after CSDS, we performed the IPA Upstream Regulator analysis ( Fig 4B ). It examines the data set for known targets of transcription regulators and compares the directionality of the gene expression differences to the information included in the Ingenuity Knowledge Base [ 27 ]. On the transcriptome level, rapamycin-insensitive companion of mTOR (RICTOR) was a potential upstream regulator with the highest bioinformatic activation score in both B6 and D2 strains among all findings in the analysis (P FDR < 0.001). RICTOR was predicted to be inhibited in the B6 susceptible mice in comparison to controls (Z-score = -4.84, P FDR = 6.37E -12 ) and activated in the same comparison in the D2 strain (Z-score = 6.55, P FDR = 8.96E -14 ). RICTOR’s function in cell growth has been extensively described [ 29 ]. The protein is also known to be involved in dendritic branch formation [ 28 , 29 ]. Furthermore, neuronal knockout mice of Rictor have reduced pre-pulse inhibition, a schizophrenia-like symptom [ 30 ].

To identify the most widely stress-affected pathways, we carried out an integrative analysis of the 12 transcriptomic and proteomic comparisons ( Fig 4A ). We found 15 pathways with at least two of the 12 comparisons significantly DE (P FDR < 0.05) and at least four of the 12 comparisons nominally significantly DE (P < 0.05). Notably, three mitochondria-related pathways ( Fig 4A , written in bold) were significantly dysregulated (P FDR < 0.05) on both transcriptome and proteome levels. Taken together, these results show significant enrichment of pathways related to mitochondrial function and transcriptional control after CSDS, and suggest that their directionality, i.e., upregulation in the B6 susceptible mice and downregulation in the D2 susceptible mice in comparison to the same-strain controls, is affected by the genetic background.

After collecting the individual mRNA, miRNA and protein expression data sets, we carried out several integrative analyses to identify biological pathways significantly affected by chronic psychosocial stress. First, we conducted Ingenuity Pathway Analysis (IPA) and gene set enrichment analysis (GSEA) for DE genes (data set A) and proteins (data set B) of the susceptible vs. control, resilient vs. control and susceptible vs. resilient comparisons in B6 and D2 strains. On the transcriptome level, the most significantly dysregulated pathway was the mitochondria-related oxidative phosphorylation pathway. It was upregulated in susceptible B6 mice compared to controls (Z-score = 4.90, P FDR = 3.38E -10 ), but downregulated in the same comparison in the D2 strain (Z-score = -5.57, P FDR = 5.47E -18 ), a result which was replicated by the GSEA ( S3 Fig ). Opposite molecular signatures for the susceptible vs. control comparison were also observed for the second most significantly dysregulated pathway, the eIF2 signaling pathway responsible for translational control, with upregulation in the B6 (Z-score = 4.47, P FDR = 2.67E -8 ) and downregulation in D2 strain (Z-score = -3.61, P FDR = 1.78E -3 ). On the protein level, the converging cAMP-mediated signaling and dopamine-DARPP32 feedback in cAMP signaling pathways were predicted to be dysregulated (P FDR < 0.01) in both B6 and D2 strains. The dopamine-DARPP32 feedback in cAMP signaling pathway was upregulated in the B6 resilient mice in comparison to the same-strain controls (P FDR = 0.001) and the D2 susceptible vs. resilient mice (P FDR = 0.003).

To identify stress-responsive miRNAs, we carried out parallel RNA-sequencing of active AGO2-associated miRNAs and their bound mRNA targets [ 25 ] in the B6 strain. We detected 233 DE miRNAs and 3565 DE mRNA targets ( Fig 3C and 3F , S2 Table ) between the studied groups (susceptible, resilient and control). Although we observed largely discrete patterns of miRNA expression in stress-susceptible and resilient B6 mice, a significant number of miRNAs were DE in the same direction in both groups indicating a general stress-related response. This finding might reflect the previously described double-role of miRNAs in both restoring homeostasis, upon stress-related environmental changes, and enforcing new phenotype-specific gene expression patterns following CSDS to better adjust to the stressful events [ 26 ].

We next examined protein abundance differences of B6 and D2 mice in the BNST after CSDS. In total, we detected 117 DE proteins in at least one of the comparisons in either or both strains ( Fig 3B and 3E ). A single protein, phosphatase 1 regulatory subunit 1B (PPP1R1B), also known as dopamine- and cAMP-regulated neuronal phosphoprotein (DARPP-32), was nominally DE in the susceptible compared to the resilient mice in both strains (P = 0.011 and Padj = 0.652 in B6 strain; P = 4.95E -9 and Padj = 1.34E -6 in the D2 strain). The protein was less abundant in the stress-susceptible B6 mice and more abundant in the D2 mice compared to the stress-resilient groups. Therefore, also the protein expression patterns after chronic psychosocial stress were significantly different between the B6 and D2 strains.

We first examined the transcriptional stress response in the BNST of stress-susceptible and resilient B6 and D2 mice. The overall number of unique DE genes in both strains was similar (n: B6 = 1638, D2 = 1441; Fig 3A and 3D , S2 Fig ), of which only 194 were common between the strains when stress-susceptible mice were compared to the controls. Notably, all of them were DE in opposite directions between B6 and D2 mice ( Fig 3A and 3D ). Since none of the individual genes were significantly DE after multiple testing correction, we performed Gene set enrichment and Gene Ontology (GO) term enrichment analyses to investigate whether these genes are part of common biological pathways. We found that these nominally DE genes were significantly enriched for mitochondrial and translation-related gene sets ( S3 Table ) and GO terms ( S4 Table ) at P FDR < 0.05. Thus, our results reveal highly divergent gene expression patterns in the two strains after CSDS.

To identify stress-associated transcriptomic and proteomic signatures in the B6 and D2 mouse strains, we profiled the BNST, a key brain region regulating anxiety. Profiling was conducted approximately one week following completion of the CSDS, when we sacrificed all mice and dissected BNSTs for analyses ( Fig 2A and 2B ). We then carried out both gene expression (RNA-seq) and proteomic (liquid chromatography-tandem mass spectrometry) profiling to identify differentially expressed mRNAs (data set A) and proteins (data set B) ( Fig 2C ). Additionally, in the B6 strain, we performed Argonaute 2 RNA immunoprecipitation-sequencing (AGO2 RIP-seq) of active microRNAs (miRNAs) and their mRNA targets (data sets D and E, respectively). AGO2 is the catalytic component of the RNA-induced silencing complex (RISC). Only mature miRNAs can be incorporated into RISC in the presence of AGO2, and serve as a guide molecule for silencing their target mRNAs [ 23 ]. Thus, AGO2 RIP-seq identifies only those miRNAs and their target mRNAs that are bound to the RISC at the time of tissue collection, providing additional specificity compared to the sequencing of all cellular miRNAs and mRNAs. Data sets A and B were collected from the same cohorts of animals, with each cohort being equally divided by the SI ratios between transcriptomic and proteomic experiments. AGO2 RIP-seq (data set D and E) was performed from an additional cohort. In all data sets, we compared the stress-resilient, stress-susceptible, and same-strain control mice. For the resilient and susceptible groups, we selected mice representing the phenotypic extremes, i.e., those with the highest and lowest SI ratios, respectively. Unless stated otherwise, “differentially expressed” (DE) mRNAs, miRNAs and proteins were defined as P < 0.05 and |FC| ≥ 1.2. For all individual findings of DE mRNAs, miRNAs and proteins, we report both nominal P-values and P-values adjusted for multiple testing by the modified Benjamini–Hochberg method (Padj) as defined in [ 24 ]. The total numbers of genes and proteins DE at P < 0.05 and Padj < 0.05 levels for each data set are presented in S2 Table .

To assess if chronic psychosocial stress affects the duration to cease escape-oriented behavior in the face of an acute stressor, we performed the forced swim test (FST) five days after the end of CSDS. The latency to immobility during the FST, used as a measure of active stress coping [ 21 , 22 ], was highly correlated with the SI ratio in the D2 defeated mice (r = 0.920, P = 0.009) ( S1 Table ). In other words, defeated D2 mice with higher resilience to psychosocial stress presented a more active coping strategy than mice with higher social avoidance. We did not observe similar correlations in the D2 control group (r = -0.130, P = 0.759) or in the B6 control or defeated mice (r = -0.079, P = 0.691 and r = 0.026, P = 0.852, respectively) ( S1 Table ). We also found that CSDS had a significant effect on the body weight of the D2, but not the B6 mice ( S1 Fig ). Overall, these results suggest that genetic background has a strong effect on stress susceptibility, with the D2 strain being more susceptible to stress-induced social avoidance than the B6 strain. Furthermore, the two strains used different strategies to cope with stress, as demonstrated by their differences in locomotor and escape-oriented behavior after, and body weight development throughout, the CSDS.

(A) Timeline of behavioral experiments including CSDS, social avoidance (SA) test, forced swim test (FST), and tissue collection. (B) Schematic illustration of 10-day CSDS. The protocol involved a daily confrontation of max. 10 min of two conspecific mice, the experimental mouse in the defeated group (black) and the aggressor mouse (white), followed by sensory contact through a clear perforated plexiglass divider for the remaining of the 24 hours. As indicated by the black arrow, the defeated mice met a new aggressor every day during the 10 days of CSDS. (C) Schematic illustration of the control condition. Two control mice were placed in a single cage separated by a clear perforated plexiglass divider. As indicated by the black arrow, every 24 hours the control mouse was placed into a cage with a new unfamiliar neighbor mouse. Unlike the defeated mice, the control mice were never in physical contact with their neighbors. (D) Schematic illustration of trial I and II of the SA test. (E) SA test results showing social interaction ratios. n = B6: Susceptible = 34, Resilient = 78, Control = 56; D2: Susceptible = 62, Resilient = 8, Control = 56. Outliers criterion: IQR > 3. (F) Percentage of mice resilient and susceptible to stress. n = B6: Susceptible = 34, Resilient = 78 D2: Susceptible = 62, Resilient = 8. B6: C57BL/6NCrl; D2: DBA/2NCrl; FST: Forced swim test; IQR: interquartile range; SA: social avoidance test.

To study how genetic background affects the behavioral response to chronic psychosocial stress, we subjected B6 and D2 mice to a 10-day CSDS, followed by the social avoidance (SA) test 24 hours later ( Fig 1 , S1 Table ). Since these strains differ in their innate social avoidance levels [ 4 ], we evaluated the behavior of the defeated mice in comparison to the same-strain controls. We divided the defeated group of animals and defined stress-susceptible mice as those with social interaction (SI) ratios below one standard deviation from the mean SI ratio of the same-strain controls, as previously described [ 4 ]. We classified all other defeated mice with ratios above those values, i.e. resembling controls, as resilient ( Fig 1E ). The strains showed distinct response to stress since 89% of D2 mice, but only 30% of B6 mice, presented social avoidance behavior, being susceptible to chronic psychosocial stress (Pearson’s chi-square, χ 2 = 60.38, P = 7.76E -14 ) ( Fig 1F ).

Discussion

To identify the key biological pathways mediating resilience and susceptibility to psychosocial stress, a risk factor for onset and recurrence of anxiety disorders, we applied a cross-species multi-omics approach. In a chronic psychosocial stress mouse model, we found differential expression of mitochondrial-related genes and proteins both in the BNST and blood cells. However, the pattern of differential expression was opposite in the B6 and D2 mouse strains. Subsequently, we tested whether the same pathways are involved in acute anxiety provocation in panic disorder patients. Our analyses revealed a consistent convergence of differentially expressed mitochondria-related pathways in the blood samples from panic disorder patients after exposure-induced panic attack. As in the stress-susceptible D2 mice, these genes were downregulated during and after panic attack in patients. Consistently, we observed significant strain-dependent stress-associated differences in mitochondrial morphology in the BNST. Taken together, our results have uncovered an evolutionarily-conserved mitochondrial signature that characterizes anxiety-related behavior in mammals.

Of the mitochondrial genes, especially those related to oxidative phosphorylation were differentially expressed in both BNST and blood cells after chronic stress in mice and during exposure-induced panic attacks in panic disorder patients. These genes, that regulate both ATP production and apoptosis, had lower expression levels in the susceptible D2 mice compared to controls and in the panic disorder patients during and after panic attack. The expression levels of the same genes were higher in the susceptible B6 mice compared to controls. In a bidirectionally bred mouse model of trait anxiety (the HAB/LAB mice), we previously observed increased expression of electron transport chain proteins in the cingulate cortex synaptosomes of the high-anxiety mice [47]. In an outbred strain rat model of social behavior, highly anxious rats that were prone to become subordinate during a social encounter with a rat with low levels of anxiety had lower levels of mitochondrial complex I and II proteins in the nucleus accumbens [48]. In a study that specifically investigated gene expression of mtDNA-encoded genes [49], four of these genes (mt-Nd1, mt-Nd3, mt-Nd6, and mt-Atp6) were downregulated after acute immobilization stress in the hippocampus. However, after chronic immobilization stress mt-Nd6 was upregulated. These effects were mediated by glucocorticoids. Thus, also previous studies have found changes in brain mitochondrial gene expression in anxiety-like behavior and after stress, but the directionality of these changes was depended on the model of anxiety and the duration of the applied stressor (e.g. acute versus chronic stress). Our results extend these previous observations by showing that the directionality of the changes may likely be influenced by the genetic background of the strain, and may be related to the innate anxiety level or stress susceptibility of these strains. Moreover, the greatest advantage of multi-omics studies is their ability to identify the most significantly affected pathways from measurements of thousands of molecules. Although mitochondrial pathways have previously been associated with anxiety, our hypothesis-free approach established their dysregulation as a major brain stress response.

We also observed differences in mitochondrial morphology and/or number of mitochondrial cross-sections in B6 and D2 mice after CSDS in the BNST. Stress-susceptible B6 mice had a larger number but shorter pre- and post-synaptic mitochondrial cross-sections compared to B6 control or resilient mice. In the D2 strain, stress-susceptible mice had slightly increased mitochondrial cross-section length/width ratio compared to the resilient mice. Previously, nocturnal aggression stress has been associated with slightly smaller pineal gland mitochondrial size in gerbils [50]. Chronic, but not acute, immobilization stress in rats leads to increased mitochondrial area in hippocampal mossy fiber terminals [51]. Thus, psychological stress is associated with morphological changes of mitochondria, but how these changes relate to mitochondrial function remains to be revealed. Cellular stress can affect mitochondrial morphology, e.g. in the form of hyperelongation, donut formation, or fragmentation in response to cytochrome c release [52]. In response to ER stress, which leads to downregulation of protein synthesis through the eIF2 pathway involved in mRNA translation, mitochondria reshape and become longer to promote cellular energy levels [53]. It has been proposed that the cumulative effect of stressors over time contribute to mitochondrial allostatic load and overload, and as a consequence lead to recalibration of mitochondrial structure and functional adaptation (e.g. activation of hormonal receptors) as well as dysregulation of gene expression, inflammatory response, and apoptosis [54–56]. In this model, mitochondria interact bidirectionally with stress mediators, but their adaptive capacity and the direction of changes related to mitochondrial function in response to stressors depend on the innate resilience of the organism and the effects of long-term programming during critical developmental periods.

In addition to the opposite directionality of mitochondria-related pathways, we found strain-specific expression patterns of differentially expressed miRNAs. Notably, miR-34c was expressed at a lower level in both blood cells and BNST of B6 resilient mice compared to controls after CSDS. Several members of the miRNA-34 family are differentially expressed in psychiatric conditions in humans: miR-34a is expressed at a higher levels in blood cells of patients with schizophrenia [57], and miR-34b-5p and mir-34c-5p in patients of major depressive disorder [58]. Furthermore, miR-34a expression is higher in the postmortem cerebellum of patients with bipolar disorder [59]. In mice, we have previously shown that both miR-34a and miR-34c are induced by acute and chronic stress in the amygdala, and that injection of miR-34c to amygdala is anxiolytic after a stress challenge [60]. In a rat model of early adolescent traumatic stress, miR-34c was upregulated in the hypothalamus [61]. Many of the miR-34 family effects on stress may be mediated through the CRFR1 [60]. We did not find CRFR1 differential expression in the BNST or blood, suggesting existence of alternative main targets in these tissues. In addition to these two miRNAs differentially expressed only in the B6 strain, we found two out of three miRNAs, miR-148a and miR-181b, to be DE in opposite directions between the strains in blood cells. miR-148a has been previously associated with panic disorder [45] and miR-181b has been identified as a possible marker for schizophrenia [46] and Alzheimer’s disease [62]. Higher expression of miR-181b has been previously correlated with an increase in mitochondrial oxidative stress and oxidative DNA damage [63], an early and systemic process in the pathophysiology of Alzheimer’s disease [64, 65]. Both miRNAs are highly conserved between mice and humans [66].

In addition to differential expression of mitochondria-related pathways and miRNAs, we found distinct molecular responses in the two strains in dopamine- and cAMP-regulated neuronal phosphoprotein (DARPP-32/PPP1R1B) and associated Dopamine-DARPP-32 Feedback in cAMP and Ca2+ Signaling pathway. Increased expression of DARPP-32 protein was previously observed in the prefrontal cortex and the amygdala of socially defeated B6 mice [36]. DARPP-32 has been proposed to act as an integrator of dopaminergic and glutamatergic signaling [67] and elevated levels of its truncated isoform were reported in schizophrenia, bipolar disorders and major depression as well as poor cognitive functioning [40]. These results implicate genetic background-dependent differences in susceptibility to chronic stress between the strains [68, 69]. Overall, our results suggest that transcriptomic response to stress is strongly dependent on the genetic background.

At the behavioral level, innately anxious D2 mice were more susceptible to CSDS than the less anxious B6 mice. While relatively few studies have examined strain differences in response to repeated stress [68, 70–72], and CSDS in particular [4, 73, 74], their findings similarly implicate heightened stress susceptibility in more anxious strains (e.g., BALB/cJ, D2). D2 resilient mice responded to CSDS with higher stress-induced locomotor activity in comparison to the D2 susceptible mice. Additionally, D2 resilient mice with higher social interaction ratios also had an elevated latency to immobility in the FST, suggesting resistance to behavioral despair or lack of adaptation aimed at energy conservation in the face of inescapable situation [75, 76]. In contrast, B6 mice did not differ in either measurement. The two strains also varied in the body weight development during and after CSDS. While the body weight decreased in the defeated D2 mice in comparison to the D2 controls and the baseline weight, the defeated B6 mice gained weight during the CSDS, similarly to the B6 controls. Taken together, we found that genetic background influences adaptation of different stress-coping strategies, as observed previously [4, 77].

For neuropsychiatric diseases, such as panic disorder, obtaining samples from the brain is essentially impossible for a reasonably large and representative sample sets. Therefore, the existence of tissue, which is more accessible and could be used as surrogate for gene expression in the central nervous system, is crucial. In our data set, the overall gene expression pattern between the BNST and blood cells was moderately correlated in both mouse strains. These results are in accordance with previous studies investigating gene expression similarities between the whole blood and multiple brain regions in humans, and suggest that gene expression is not perfectly correlated between brain and blood, but may be useful for studying certain pathways, e.g. related to translational control [78]. On the pathway level, oxidative phosphorylation-related genes were differentially expressed both in the mouse BNST and blood cells, and in the panic disorder patient blood cells. These genes were upregulated in the defeated B6 and downregulated in the defeated D2 mice compared to controls both in the BNST and blood cells. Strikingly, this pathway was downregulated in panic disorder patients directly and 24 h after an exposure-induced panic attack. It is interesting that the gene expression pattern of the patients resembles that of the more stress-susceptible mouse strain, suggesting that stress-susceptibility may involve a general modulation of genes associated with mitochondrial function. There are several caveats in comparing a post-mitotic brain region to a peripheral biospecimen, which are very different in terms of mitochondrial metabolism, but nonetheless, we observed a consistent and converging signature that warrants further investigations. Altogether, our data reinforce the utility of cross-species approaches in the identification of the biological basis of human anxiety disorders but advises careful selection of mouse strains for the best translatability of the findings.

Although we concentrated on mitochondrial pathways and anxiety susceptibility in our follow-up experiments, we also found several pathways in the mouse transcriptomic and proteomic experiments associated with stress resilience. Stress resilience is considered an active, adaptive response to stress, not merely a lack of maladaptive symptoms [79]. Studying of resilience is challenging in humans due to a lack of well-controlled cohorts. However, translational studies on stress resilience would be highly warranted to understand the underlying mechanisms providing putative means to enhance resilience in stress-susceptible individuals.

In conclusion, our cross-species multi–omics approach found a systemic evolutionarily conserved mitochondrial response in anxiety-related behaviors in mice and in panic disorder patients. We have produced a large amount of mRNA, miRNA, and protein expression data and made it publicly available to provide a resource for formulation of additional hypotheses on psychosocial stress-induced anxiety in mice, and panic disorder in humans. Our converging findings on stress-susceptibility in both brain and blood cells indicated dysregulation of translation and mitochondrial-related pathways. Further functional studies underlying the mechanisms behind the observed differences in mitochondrial morphology and the dysregulation of the mitochondria-related genes will provide much needed insight into the molecular basis of panic disorder and other anxiety disorders, a critical step in developing future targeted therapies.