Abstract
Anxiety disorders often manifest in genetically susceptible individuals after psychosocial stress, but the mechanisms underlying these gene-environment interactions are largely unknown. We used the chronic social defeat stress (CSDS) mouse model to study resilience and susceptibility to chronic psychosocial stress. We identified a strong genetic background effect in CSDS-induced social avoidance (SA) using four inbred mouse strains: 69% of C57BL/6NCrl (B6), 23% of BALB/cAnNCrl, 19% of 129S2/SvPasCrl, and 5% of DBA/2NCrl (D2) mice were stress resilient. Furthermore, different inbred mouse strains responded differently to stress, suggesting they use distinct coping strategies. To identify biological pathways affected by CSDS, we used RNA-sequencing (RNA-seq) of three brain regions of two strains, B6 and D2: medial prefrontal cortex (mPFC), ventral hippocampus (vHPC), and bed nucleus of the stria terminalis (BNST). We discovered overrepresentation of oligodendrocyte (OLG)-related genes in the differentially expressed gene population. Because OLGs myelinate axons, we measured myelin thickness and found significant region and strain-specific differences. For example, in resilient D2 mice, mPFC axons had thinner myelin than controls, whereas susceptible B6 mice had thinner myelin than controls in the vHPC. Neither myelin-related gene expression in several other regions nor corpus callosum thickness differed between stressed and control animals. Our unbiased gene expression experiment suggests that myelin plasticity is a substantial response to chronic psychosocial stress, varies across brain regions, and is genetically controlled. Identification of genetic regulators of the myelin response will provide mechanistic insight into the molecular basis of stress-related diseases, such as anxiety disorders, a critical step in developing targeted therapy.
- anxiety
- chronic social defeat stress
- inbred mouse strain
- myelin
- RNA-sequencing
- transmission electron microscopy
Significance Statement
Chronic psychosocial stress is a well-established risk factor for anxiety disorders, but the development of targets for therapeutic intervention is limited by ignorance of the underlying molecular and cellular mechanisms. We used inbred genetically defined mice to identify neurobiological pathways that underlie stress-induced social avoidance (SA), a type of anxiety. We found genetically controlled differences in myelin-related gene expression in stress-exposed mice, with concurrent differences in myelin thickness, suggesting that myelin plasticity is a major stress response of the brain. The adaptive response to stress may increase or decrease myelin thickness, depending on the demands of the specific circuit. Our findings provide a foundation for the identification of specific genetic regulators of chronic stress-induced myelin plasticity.
Introduction
Anxiety disorders, including panic disorder, social anxiety disorder, specific phobias, and generalized anxiety disorder, are the most common mental disorders with a prevalence of 14% (Wittchen et al., 2011). Human genetic studies of anxiety disorders have confirmed a modest heritability and considerable environmental component (Hettema et al., 2001). However, identification of replicated risk variants is challenging due to genetic heterogeneity, environmental factors that cannot be well-controlled in human settings, and poor availability of large patient cohorts with accurate phenotypes (Hovatta and Barlow, 2008; Smoller, 2016). Consequently, many investigators have resorted to animal models, which allow controlled experiments on both genetic and environmental risk factors, to reveal the biological mechanisms underlying these diseases.
Chronic psychosocial stress is a well-established risk factor for anxiety disorders (Kemeny and Schedlowski, 2007; Moffitt et al., 2007; Donner et al., 2012). It can be modeled in mice by the chronic social defeat stress (CSDS) paradigm, that has etiological, predictive, and face validity for affective and anxiety disorders (Hammels et al., 2015). It leads to social avoidance (SA) and long-term plastic changes in the brain (Avgustinovich et al., 2005; Krishnan et al., 2007). However, only a portion of mice exhibit SA after CSDS, and thus the defeated animals can be divided into stress-susceptible and resilient. Comparison of these groups allows investigation of the mechanistic basis of stress-induced SA and anxiety-like behavior and resilience to it. Understanding the underlying risk factors and promotion of the resilience factors in susceptible individuals should facilitate development of secondary prevention methods of anxiety disorders after traumatic events, and selective pharmacological treatment of anxiety (Howlett and Stein, 2016).
Most insight into the molecular mechanisms of CSDS comes from the C57BL/6 mouse strain, in which specific genetic, epigenetic, and neurophysiological mechanisms underlie stress susceptibility and resilience (Han and Nestler, 2017). Because genetic background strongly modulates mouse behavior (Threadgill et al., 1995; Hovatta et al., 2005; Sittig et al., 2016), hypotheses of any mechanisms should be tested across different mouse strains. Such studies are also critical for translating the results to genetically heterogeneous humans. CSDS has different behavioral and physiologic consequences in C57BL/6J and BALB/c strains (Razzoli et al., 2011a,b; Savignac et al., 2011), reflecting underlying genetic, and by extension, gene expression differences between the strains. Therefore, comprehensive unbiased transcriptomic approaches applied to mouse strains with different responses to CSDS should reveal gene-environment interactions that explain why some individuals are more susceptible to stress-induced maladaptive behaviors than others.
To examine how genetic background affects the behavioral response to chronic stress, we performed behavioral testing in four inbred mouse strains after CSDS. We selected C57BL/6NCrl (B6) and BALB/cAnNCrl (BALB) strains, which differ in behavioral and metabolic responses to CSDS (Razzoli et al., 2011a,b; Savignac et al., 2011), and DBA/2NCrl (D2) and 129S2/SvPasCrl (129) as they are sensitive to stress-induced anxiety- and depression-like behavior in general (Ducottet and Belzung, 2005; Millstein and Holmes, 2007). To identify the major underlying biological pathways, we performed unbiased transcriptomic analysis in the B6 and D2 strains that showed the largest differences in susceptibility to stress. We studied three brain regions, the medial prefrontal cortex (mPFC), ventral hippocampus (vHPC), and bed nucleus of the stria terminalis (BNST), which critically regulate anxiety and the stress response (Garakani et al., 2009; Calhoon and Tye, 2015; Tovote et al., 2015), and are activated by CSDS (Vialou et al., 2014; Laine et al., 2017). In the follow-up gene expression and histologic analyses, we found that myelination was significantly altered after stress in a strain and brain region-dependent manner. Our results illustrate that genetic background has a large effect on both the behavioral and brain transcriptomic response to chronic psychosocial stress. Moreover, we demonstrate that the pattern of stress-induced myelination changes is dependent on the genetic background and varies across brain regions.
Materials and Methods
Animals
We ordered five-week-old male mice from four inbred strains [DBA/2NCrl (D2), 129S2/SvPasCrl (129), BALB/cAnNCrl (BALB), and C57BL/6NCrl (B6); Charles River Laboratories] for all CSDS experiments and let them acclimatize for 10 d before CSDS, housed in groups in a temperature (22 ± 2°C) and humidity (50 ± 15%) controlled facility on a 12/12 h light/dark cycle (lights on 6 A.M. to 6 P.M.). From the end of CSDS to the time of dissection, all mice were single-housed. As aggressors for CSDS, we used male Clr-CD1 mice (CD1, Charles River Laboratories), aged 13–26 weeks. All mice had ad libitum access to food and water throughout the experiment, except for the durations of behavioral tests. Aspen chip bedding (Tapvei Oy) in the cages was changed weekly (except for the duration of the CSDS) and standard environmental enrichment [aspen strips as nesting material (Tapvei Oy) and an aspen brick (Tapvei Oy)] was provided throughout the experiment. Animal procedures were approved by the Regional State Administration Agency for Southern Finland (ESAVI-3801-041003-2011 and ESAVI/2766/04.10.07/2014) and conducted in accordance to directive 2010/63/EU of the European Parliament and of the Council, and the Finnish Act on the Protection of Animals Used for Science or Educational Purposes (497/2013).
Behavioral experiments
CSDS
We conducted CSDS as previously described (Golden et al., 2011; Laine et al., 2017). Briefly, aggressor CD1 mice were first checked for appropriate aggression levels during a 3-d screening before all social defeat experiments. The selected aggressors had to meet the following criteria: attack in at least two consecutive sessions, had latency to attack of <90 s and do not attack within 1–5 s in any of the sessions. For CSDS, each defeated mouse was placed into the cage of an aggressor mouse for max. 10 min. The mouse was then transferred to another compartment of the cage, separated from the aggressor by a perforated Plexiglas wall, for 24 h. This procedure was repeated for 10 d, and each day the experimental mouse encountered a novel aggressor. Physical contact time was shortened in the case of severe physical aggression. Control mice were housed in similar cages but with another control mouse as a cage-mate, and without physical contact, switching cage-mates daily. The day after the last defeat session, all mice were separated into single-housed cages and maintained like this until dissection. The order of testing for all consecutive behavioural tests was randomized.
SA test
Twenty-four hours after the last social defeat session, at the start of the dark phase of the light cycle, we tested both defeated and control mice in the SA test. All animals were brought to the experiment room at least 30 min before the start of the test, and animals performing the test were separated from experimenters and the remaining animals by a screen. For the first trial (no-target), the mouse was placed in the center of an open arena (42 × 42 cm) with an empty perforated Plexiglas cylinder located next to one of the walls. The movements of the mouse were tracked using a camera and EthoVision XT10 software (Noldus Information Technology) for 150 s, after which the mouse was placed back into the home cage. The arena was cleaned and the Plexiglas cylinder was replaced with another one containing an unfamiliar CD1 social target mouse. The test mouse was then placed back in the middle of the arena and tracked for 150 s (target trial). The amount of time the mice spent in the interaction zone (IZ), defined as a semicircle (370 cm2) around the perforated Plexiglas cylinder, was measured and a social interaction (SI) ratio calculated by dividing the IZ time of the social target trial with the IZ time of the no-target trial, multiplied by 100. Thus, a low SI ratio indicates high SA.
To account for strain differences in baseline social behavior, we assessed the response of defeated mice in relation to same-strain controls, following the statistical approach previously implemented by Nasca et al. (2015). We calculated mean SI ratios of control mice from each strain based on several cohorts (n: 129 = 8, BALB = 40, B6 = 126, D2 = 114). We used log-transformation to normalize the distribution, removed outliers (>3 IQRs from the median), and divided the defeated mice to stress-resilient (resembling controls) and susceptible (showing SA) based on SI ratios, with the border determined as the controls’ mean score minus 1 SD. SI ratio border values for each strain were: 129 = 62.68, BALB = 81.76, B6 = 76.49, D2 = 105.99.
Body weight
Body weight was recorded at the beginning of CSDS (day 1) and on every second day throughout the CSDS (days 2, 4, 6, 8, and 10). The amount of body weight gain was calculated as the difference between the first and the last measurement.
Open field (OF) test
We assessed spontaneous locomotor activity and anxiety-like behavior with an automatic MedAssociates system at the start of the light phase of the light cycle. Individual mouse cages were brought to the experimental room in two groups. After 30 min adaptation to the experimental room (lit at 100 lux), the mouse was released to the corner of the experimental chamber (27 × 27 × 20 cm, transparent walls and white floor virtually divided into a 19 × 19 squares grid) and allowed to explore freely for 5 min.
Elevated zero maze (EZM)
EZM was performed at the start of the light phase of the light cycle. All animals were brought to the experiment room at least 30 min before the test, and animals performing the test were separated from experimenters and the remaining animals by a screen. The apparatus consisted of plastic annular runway (diameter = 50 cm, width = 5 cm) elevated 40 cm above the floor. The runway was divided into four sectors: two open sectors opposing each other and two opposing closed sectors protected by inner and outer non-transparent walls (height = 15 cm). After 30 min of adaptation to the experimental room (dimly lit at 15–20 lux), the mouse was placed in the middle of one of the closed sectors and allowed to explore the maze freely for 5 min. The mouse was video-tracked with EthoVision XT10.
Forced swim test (FST)
All animals were brought to the experiment room at least 30 min before the test, and animals performing the test were separated from experimenters and the remaining animals by a screen. After adaptation to the experimental room (150 lux), the mouse was placed in a glass cylinder (diameter = 18 cm, height = 25 cm) filled with water (room temperature) to the height of 15 cm. The immobility time (passive floating) was detected with EthoVision XT10 system for 6 min with 2-min time bins. Data from the last 4 min were used.
Gene expression profiling
Dissections
We dissected the mPFC (B6: Con n = 6, Res n = 6, Sus n = 6; D2: Con n = 6, Sus n = 8), BNST (B6: Con n = 5, Res n = 5, Sus n = 5; D2: Con n = 5, Res n = 3, Sus n = 5), and vHPC (B6: Con n = 6, Res n = 8, Sus n = 3; D2: Con n = 6, Sus n = 5) 6–8 d after the last CSDS (mice aged eight weeks). Mice were killed by cervical dislocation between 8 and 11 A.M. to avoid circadian differences in gene expression, and the order was counterbalanced across groups (resilient, susceptible, and control mice). Dissections were performed on a sterile chilled Petri dish within 7 min, and tissue was flash frozen in liquid N2.
RNA-sequencing (RNA-seq)
Total RNA was extracted with TriReagent (Molecular Research Center Inc.) and RNA quality was assessed with a 2100 Bioanalyzer (Agilent Technologies) using Agilent RNA 600 Nano Chip kit (Agilent Technologies). rRNA was depleted with Ribo-Zero Gold rRNA Removal kit (Illumina Inc; mPFC and vHPC) or custom Insert Dependent Adaptor Cleavage (InDA-C) primers (BNST). RNA was fragmented using the S2 ultrasonicator (Covaris Inc.) and sequencing libraries were prepared with Nextera (Illumina; vHPC), ScriptSeq v2 (Epicentre; mPFC), or Ovation Universal RNA-Seq System (NuGEN; BNST) RNA-seq library preparation kits. Libraries were size-selected with Pippin Prep (Sage Science) and sequencing was conducted on HighSeq 2000 (vHPC, paired-end 91 bp, Illumina) or NextSeq 500 platforms (mPFC and BNST, single-end 96 bp; Illumina).
The RNA-seq reads were trimmed for adapters with Cutadapt v1.8.3 (vHPC) and FastX toolkit (mPFC, BNST) and PCR duplicates were removed with PRINSEQ v0.20.4. Reads were aligned using STARv 2.5.0c (Dobin et al., 2013) with default settings to mouse genome GRCm38, and annotated to gene exons with HTSeq v0.6.1 (Anders et al., 2015) using GTF release 86 (update 2016-10).
Differential expression (DE) analysis
DE analysis was conducted using limma eBayes (Ritchie et al., 2015; Phipson et al., 2016) comparing resilient and susceptible mice to same-strain controls within brain regions. RNA-seq data were filtered to remove low-abundance genes, keeping genes with at least 1 count per million (CPM) in at least six samples (Anders et al., 2015). Subsequently, the data were normalized with voom (Law et al., 2014), and adjusted for sequencing (vHPC) and library preparation batches (vHPC, mPFC, and BNST) with ComBat (Johnson et al., 2007). To identify top DE genes, we calculated -log10(p)*logFC for each gene and ranked them accordingly (Xiao et al., 2014).
Gene expression data are available in Gene Expression Omnibus (GEO; GSE109315).
Rank-rank hypergeometric overlap (RRHO) test (Plaisier et al., 2010)
RRHO infers pair-wise similarity between two DE result lists, where genes are ranked by DE (-log10(p)*logFC) between resilient and control, or between susceptible and control mice, by calculating significance of overlapping genes at different rank bins (Plaisier et al., 2010). We used step size of 100 genes to bin the genes and applied the same -log(p) scale to comparisons within brain regions. Significant overlap of genes in rank groups containing up- and downregulated genes in both lists (lower left corner and upper right corner of square matrix, respectively, see Fig. 2C,D) shows a shared transcriptome-wide gene expression pattern in response to CSDS.
Gene set enrichment analysis (GSEA)
We conducted GSEA using the GSEA Preranked module implemented in GSEA Desktop v3.0 (Mootha et al., 2003; Subramanian et al., 2005) and the curated gene sets (C2) of the Molecular Signature Database (MSigDB) v6.0 (http://www.broad.mit.edu/gsea/). The pre-ranked GSEA was performed with 1000 permutations. The top five gene sets with the highest positive and negative normalized enrichment scores (NESs; pFDR < 0.05) within each comparison were selected for further analysis. From the selected top up- and downregulated gene sets, all present in at least two comparisons and two brain regions were visualized using Circos software (Krzywinski et al., 2009). The overlap between the top enriched gene sets, presented on the Circos plot, was further investigated with the hypergeometric test implemented in the MSigDB v6.0.
Gene ontology (GO) term enrichment
We analyzed GO term enrichment separately for top 300 upregulated and 300 downregulated differentially expressed genes from each comparison with topGO (Alexa et al., 2006), using the weight01 model to account for GO term dependencies. In RNA-seq, long transcripts yield more read counts and are more easily passed through low-abundance filtering than short genes with the same expression level, and genes with high number of counts have greater statistical power being detected as DE than genes with low number of counts. To minimize these selection biases (Young et al., 2010), we matched the background genes, i.e., the gene universe used to compare the top DE genes with, with the top DE genes using R package genefilter (Gentleman et al., 2017), resulting in 6882 (±290 SEM) background genes.
Visualizing oligodendrocyte (OLG) progenitor cell (OPC) and OLG marker genes
We manually curated a list of OPC and OLG-specific marker genes, based on prior publications (Zhang et al., 2014; Marques et al., 2016). Figures were constructed based on previously published scripts (Haarman et al., 2014).
q-RT-PCR
q-RT-PCR was applied to validate five myelin-related genes from RNA-seq. We used published primers to amplify Mobp and Plp1 (Liu et al., 2012) and designed primer pairs (5’-3’ forward, reverse) using NCBI primer designing tool (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) for Opalin (ACTGCCATCGAATACGACATC, CCTCTACGGGCTCATCATCG), Ermn (AACCAGGCAGGAGACAACTG, GATGGCCTGGTGAACAACGA), and Mbp (ACACACGAGAACTACCCATTATGG, AGAAATGGACTACTGGGTTTTCATCT). RNA samples of mPFC and BNST were the same as used in RNA-seq, and for vHPC samples overlapped by 37.5% (12/32 samples); 250 ng of DNase I (Thermo Scientific)-treated total RNA was converted to cDNA with iScript select cDNA synthesis kit (Bio-Rad Laboratories) and amplified with 250 nM primers in CFX384 Real-Time PCR cycler using IQ SYBR Green supermix (Bio-Rad Laboratories). Expression levels were normalized to Ppib (GGAGATGGCACAGGAGGAAA, CCCGTAGTGCTTCAGCTTGAA). Each reaction was run in triplicate and relative expression level was calculated using a standard curve (7.15, 10.0, 5.0, 2.0, 1.0, 0.5, and 0.25 ng of cDNA) present on each assay plate with CFX Manager (Bio-Rad Laboratories). Statistical analysis (Pearson’s r) was conducted with GraphPad Prism v7.02 (GraphPad Software Inc.).
Immunohistochemistry
Mice were anaesthetized 6–8 d after CSDS with a lethal dose of pentobarbital (Mebunat Vet 60 mg/ml, Orion Pharma) and transcardially perfused with 4% paraformaldehyde (PFA) in PBS. After postfixation in 4% PFA overnight (+4°C), the brains were cut into 20-µm coronal sections using a Leica VT-1200S vibratome (Leica Biosystems) and stored at –20°C free-floating in cryoprotectant. Sections were washed 3× in PBS and mounted on Superfrost Ultra Plus (ThermoFisher Scientific) slides. We performed antigen retrieval by submerging the slides in 0.01 M citrate buffer, heated to a boil for 20 min. Slides were blocked in 2.5% BSA in 0.5% PBST + 7.5% normal goat serum, followed by primary antibody incubation with mouse anti-CNPase (1:250, Merck Life Science, #MAB326R) overnight at +4°C. Slides were washed in PBS before goat anti-mouse Alexa Fluor 488 secondary antibody incubation (1:400, ThermoFisher Scientific, #A-11029) for 2 h at room temerture. After washing in PBS, we coverslipped the slides with Vectashield + DAPI mounting medium (Vector Laboratories, #H-1200). We acquired images with ZEISS Apotome.2 system (Zeiss) and analyzed them with ImageJ software (National Institutes of Health). We measured corpus callosum thickness on both sides of the midline from each section (two to six sections per animal, distance from bregma between 0.22 and -0.10) using ImageJ, and calculated their mean.
Transmission electron microscopy (TEM)
Mice were anaesthetized 6–8 d after CSDS with a lethal dose of pentobarbital (Mebunat Vet). We transcardially perfused the mice with PBS followed by fixation with 100 ml 2% glutaraldehyde (GA)/2% PFA in 0.1 M sodium cacodylate (NaCac) buffer (GA, PFA, and NaCac: Sigma Aldrich), heated to +37°C. The brains were postfixed in the same fixative for 2–4 h and immersed in 0.1 M NaCac buffer for 2–24 h, both at +4°C. We cut them into 200-μm sagittal slices with a Leica VT-1200S vibratome (Leica Biosystems) in 0.1 M phosphate buffer. Regions of the mPFC, BNST, and vHPC were cut manually by using anatomic landmarks (Extended Data Fig. 5-1).
We postfixed the sections in osmium tetroxide [1%, +1.5% K4[Fe(CN)6] in 0.1 M NaCac] for 2 h at +4°C, stained en bloc with uranyl acetate for 1 h at +4°C, dehydrated with EtOH and acetone, and embedded them into hard Epon; 60- to 70-nm sections were cut with an ultramicrotome and placed on copper grids for microscopy and stained with lead citrate. We imaged the ultrathin sections using a Jem-1400 transmission electron microscope (Jeol) and randomly selected and imaged myelinated axons at 5000× magnification.
We measured myelin thickness, axon diameter, and g ratio using ImageJ. We calculated the diameter by measuring the area of the whole fiber and the area of the axon (inside the compacted myelin sheath). The diameters of geometric circles with the same areas were calculated for both parameters. We calculated the g ratio by dividing the diameter of the axon with the diameter of the whole fiber. Myelin thickness was measured at three fully compacted positions and their average calculated.
Statistical analysis
Statistical analyses were conducted with SPSS Statistics 24 (IBM) or GraphPad Prism 7.02 (GraphPad Software Inc.). Planned post hoc comparisons (control vs resilient, control vs susceptible, and resilient vs susceptible) were conducted by Fisher’s LSD. For these tests we report nominal p values, evaluated for significance against an α-level adjusted for multiple corrections with test-wise Bonferroni correction (Table 1). Only p values which survive this correction are shown.
Group differences in TEM data were assessed using generalized estimating equations (GEEs) to control for within-subject dependencies of individual axons measured from the same animal (ranges for each region: mPFC = 93–104, BNST = 31–69, and vHPC = 54–62 axons per animal). We selected this approach due to the low number of animals per group, which could not be reliably analyzed by ANOVA. GEE has been proposed as a suitable approach for analyzing data with non-independent features (Hanley et al., 2003), such as axons measured from the same individual. Pair-wise contrasts were computed for comparing groups with Fisher’s LSD and significance determined against the Bonferroni corrected α-level as above.
Data analysis of RNA-seq data were performed as described above. Multiple testing correction was done by the Benjamini–Hochberg method (Benjamini and Hochberg, 1995).
Results
Genetic background influences behavioral response to CSDS
To determine whether genetic background affects the behavioral response to psychosocial stress, we conducted 10-d CSDS in mice from four inbred strains, D2, 129, BALB, and B6 (Fig. 1A). Twenty-four hours after the last defeat session, we conducted the SA test to assess their SA phenotype. To account for the strain differences in baseline social behavior during SA test, we evaluated the response of defeated mice in relation to same-strain controls. We divided the defeated mice into stress-susceptible and resilient groups, considering mice with SI ratios within 1 SD, or above, of the same-strain control mean as resilient (i.e., resembling the control mice), and those with the SI ratio below 1 SD from the mean as susceptible. We observed a significant difference in the distribution of susceptible and resilient mice between the strains (χ2 = 63.401, p = 1.102 × 10−13), with the B6 defeated mice being mostly resilient and the 129, BALB, and D2 mice mostly susceptible to stress (Fig. 1B). During the social target trial of the SA test, susceptible mice from all strains spent significantly less time in the IZ (Fig. 1C) and more time in the corners of the arena (Fig. 1D) than during the no-target trial.
Extended Data Figure 1-1
Exact p values for all behavioral tests. The test used for each parameter depicted in Figure 1C–I are outlined along with results of post hoc contrasts. D2: DBA/2NCrl strain; BALB: Balb/cAnNCrl; 129: 129S2/SvPasCrl; B6: C57BL/6NCrl. Download Figure 1-1, XLSX file.
We next determined how CSDS influences other behaviors. To assess locomotor behavior, we measured the distance moved during the no-target trial of the SA test (Fig. 1E) and the OF test (Fig. 1F). B6 and D2 defeated mice moved significantly less than control mice during both tests. We did not observe differences in distance traveled between defeated and control 129 or BALB mice. For the B6 strain, we also performed the EZM and FST to assess anxiety-like and despair behavior, respectively. Susceptible mice showed increased anxiety-like, but not despair behavior, compared to controls (Fig. 1G,H). To study metabolic effects of CSDS, we measured body weight before and after CSDS (Fig. 1I). All B6 mice gained weight (mixed ANOVA post hoc comparison before versus after CSDS: controls p = 9.98 × 10−28, resilient p = 2.21 × 10−26 and susceptible p = 6.23 × 10−20). BALB and 129 controls gained weight (p = 4.90 × 10−5 and p = 0.012, respectively), while the weight of the defeated mice of these strains did not change during CSDS. Both D2 resilient and susceptible mice lost weight (p = 0.004 and p = 2.53 × 10−6, respectively).
Oligodendrocyte (OLG)-related genes are differentially expressed in response to stress
To establish which biological pathways were affected by chronic stress, we conducted RNA-seq one week after CSDS in mPFC, vHPC, and BNST (Fig. 2A Extended Data Fig. 2-1). We selected the B6 and D2 strains for this analysis as they represented the phenotypic extremes in the proportions of susceptible and resilient mice. We were not able to analyze gene expression levels of resilient D2 mice for mPFC and vHPC due to low number of animals in this group. We always compared the stress-susceptible or resilient mice to the same-strain controls. We first determined the overlap of the top 300 up- and downregulated genes between the strains in each brain region (Fig. 2B; Extended Data Fig. 2-2), followed up by the RRHO analysis, which shows the overall similarity and direction of DE of all genes (Fig. 2C,D; Extended Data Fig. 2-3). In the mPFC, only 26 (2.2%) of the top differentially expressed genes between the susceptible versus control mice were shared between the B6 and D2 strains, and the RRHO analysis confirmed the highly divergent stress response of the two strains. In the BNST, the transcriptomic response of the two strains was marginally more similar as resilient B6 and D2 mice shared 67 (5.9%), and susceptible B6 and D2 mice shared 101 (9.2%) of the top genes. We detected the greatest overlap of the gene expression response between the strains in the vHPC, where B6- and D2 susceptible mice shared 390 (48.1%) of the top genes. Unlike in the mPFC or BNST, the stress effect was stronger than the strain effect in the vHPC. However, in the vHPC, several genes were downregulated both in the B6 susceptible and resilient mice.
Extended Data Figure 2-1
SI ratios, division of mice to phenotypic groups, library preparation and sequencing batches, and the number of analyzed RNA-seq reads. mPFC: medial prefrontal cortex; vHPC: ventral hippocampus; BNST: bed nucleus of the stria terminalis; B6: C57BL/6NCrl; D2: DBA/2NCrl. Download Figure 2-1, XLSX file.
Extended Data Figure 2-2
Top 300 upregulated and top 300 downregulated genes as ranked by -log10(p)*|logFC| for each comparison shown in Figure 2B. mPFC: medial prefrontal cortex; vHPC: ventral hippocampus; BNST: bed nucleus of stria terminalis; B6: C57BL/6NCrl; D2: DBA/2NCrl; Res: resilient; Sus: susceptible; Con: control. Download Figure 2-2, XLSX file.
Extended Data Figure 2-3
Overlap of the top differentially expressed and OLG-related genes after CSDS. A, B, Overlap of the top 300 downregulated (A) and top 300 upregulated (B) genes between resilient versus control and susceptible versus control mice, separately for each brain region (Fig. 2B). C, Merged heat map showing the expression fold change (FC) of resilient versus control and susceptible versus control groups. Genes belonging to the LEIN_OLIGODENDROCYTE_MARKERS gene set in the GSEA (Fig. 2E) are shown. FCs are shown only for genes with nominal p < 0.05. Genes which did not pass the cut-off are marked in grey (NA). mPFC: medial prefrontal cortex; vHPC: ventral hippocampus; BNST: bed nucleus of the stria terminalis; B6: C57BL/6NCrl; D2: DBA/2NCrl; Con: control; Res: resilient; Sus: susceptible. Download Figure 2-3, PDF file.
Extended Data Figure 2-4
Pathway analysis results from GSEA for Figure 2E. Only pFDR < 0.05 are shown in the table. pFDR: a false discovery rate corrected p value; pFWER: familywise-error rate. mPFC: medial prefrontal cortex; vHPC: ventral hippocampus; BNST: bed nucleus of the stria terminalis; B6: C57BL/6NCrl; D2: DBA/2NCrl; Res: resilient; Sus: susceptible; Con: control. Download Figure 2-4, XLSX file.
Extended Data Figure 2-5
Significantly (p < 0.05) enriched GO terms within the top 300 upregulated and top 300 downregulated genes for each comparison, converging with results from GSEA (Fig. 2E). mPFC: medial prefrontal cortex; vHPC: ventral hippocampus; BNST: bed nucleus of the stria terminalis; B6: C57BL/6NCrl; D2: DBA/2NCrl; Res: resilient; Sus: susceptible; Con: control; BP: biological process; CC: cellular component; MF: molecular function. Download Figure 2-5, XLSX file.
To ask which biological pathways were affected by CSDS, we conducted GSEA and GO term enrichment analysis. GSEA showed significant enrichment of several gene sets (Extended Data Fig. 2-4), of which aging- and OLG-related sets were enriched in nearly all comparisons in all brain regions and both strains (Fig. 2E; Extended Data Fig. 2-3C ). Although functionally diverse, the genes included in the aging-related gene sets were significantly overrepresented in the “Lein OLG Markers” gene set (pFDR = 1.2 × 10−16). In accordance with GSEA, the most enriched GO terms were related to myelination and OLG development, in particular among the downregulated genes of B6 susceptible mice in the mPFC (Extended Data Fig. 2-5). These combined results from both enrichment analyses prompted us to further investigate OLG-related genes.
Mature OLGs develop from OPCs, which persist even in the adult brain as committed precursors. We asked whether either OPC or OLG cell populations dominantly contributed to the observed DE. The transcriptomic signature associated with CSDS was stronger in the mature OLG markers (Fig. 2G) than in the OPC markers (Fig. 2F). To validate the RNA-seq findings, we analyzed Opalin, Ermn, Mobp, Plp1, and Mbp gene expression levels with q-RT-PCR in mPFC, BNST, and vHPC and found high correlation with RNA-seq and q-RT-PCR measurements of these myelination-related genes (mean r = 0.82; Fig. 3).
Myelin-related variation after stress is not observed globally in the brain
To test whether differences in myelin-related gene expression were observed across the brain, we measured expression levels of Opalin, Ermn, Mobp, Plp1, and Mbp in the whole cortex (lacking the mPFC), hypothalamus, and dorsal hippocampus (dHPC) by q-RT-PCR. Using mixed ANOVA we determined that there was no significant main effect by the group (control, resilient or susceptible) on myelin-related gene expression in any of the brain regions of either strain. While Opalin expression was lower in D2 susceptible mice compared to controls in the hypothalamus as determined by post hoc comparison (p = 0.006), the expression levels of the other genes did not differ between the stressed and control mice in any region (Fig. 4A–F). We also measured the thickness of the corpus callosum in brain sections stained with a myelin-binding antibody (anti-CNPase) (Fig. 4G). We observed no differences in stress-susceptible or resilient mice compared to controls (Fig. 4H,I).
Myelin thickness and g ratio differ in brain region and genetic background-dependent manner
To determine whether CSDS affects myelination, we conducted TEM of myelinated axons in mPFC, BNST, and vHPC. In addition to analyzing axons of different diameters within each brain region and stress group together (Extended Data Fig. 5-2), we divided the axons to three size groups given that different types of neuronal projections may differ in axon diameter (Fig. 5; Extended Data Fig. 5-3). We discovered several brain region- and strain-specific differences between stress groups in g ratio, i.e., the ratio of the inner axonal diameter to the total outer diameter, myelin thickness, and axon diameter. Overall, D2 resilient mice had higher g ratio and thinner myelin in the mPFC compared to susceptible mice. B6 susceptible mice had thinner myelin in the vHPC compared to controls and thicker myelin in the BNST compared to resilient mice (Extended Data Fig. 5-2). In addition to these general findings, we found significant differences between stress susceptible or resilient mice compared to controls in axons of certain diameter (Fig. 5). We also observed modestly but significantly smaller axon diameter (without the myelin sheath) in the mPFC of D2 susceptible mice compared to controls (Extended Data Fig. 5-2).
Extended Data Figure 5-1
Regions dissected for TEM. Purple shaded squares outline the dissected samples from 200-µm sections. Atlas outlines are based on Franklin and Paxinos (2008), and the distance from the midline (sagittal) in millimeters is shown below each image. Download Figure 5-1, PDF file.
Extended Data Figure 5-2
Effect of CSDS on g ratio (A) and myelin thickness (B) when comparing all axons without division into axon size categories (for results with division into small, medium, and large axon size categories, see Fig. 5) and on axon diameter (C). D, Representative TEM images for each group, scale bar = 0.5 µm. Error bars ± 1 SEM. All nominal p values surviving Bonferroni correction against α = 0.0167 are shown. Download Figure 5-2, PDF file.
Extended Data Figure 5-3
Details of axons assessed by TEM (depicted in Fig. 5). Calculation of axon diameter category boundaries (small, medium, large) in each brain region, and the number of animals and axons per category and group. mPFC: medial prefrontal cortex; vHPC: ventral hippocampus; BNST: bed nucleus of the stria terminalis; B6: C57BL/6NCrl; D2: DBA/2NCrl; Res: resilient; Sus: susceptible; Con: control. Download Figure 5-3, XLSX file.
Discussion
We established that the behavioral and brain transcriptomic responses to chronic psychosocial stress are genetically controlled. We discovered that BALB, 129, and D2 mice were more susceptible to chronic stress than B6 mice. We demonstrated, by unbiased RNA-seq, that after CSDS the most significantly affected gene sets and biological pathways in the mPFC, vHPC, and BNST were related to myelination. Consistently, we observed significant brain region and strain-dependent differences in myelin thickness after stress. Neither myelin gene expression in other cortical regions or the dHPC, nor corpus callosum thickness, differed between stressed and control mice, suggesting no overall white matter changes due to CSDS.
We first established the CSDS paradigm in four inbred mouse strains and demonstrated that genetic factors control their adaptive behavior to stress. Consistently with prior findings (Razzoli et al., 2011a; Savignac et al., 2011), BALB mice were more sensitive to CSDS-induced SA than B6 mice. Additionally, innately anxious D2 and 129 mice were also highly susceptible to CSDS, suggesting that anxious strains may in general be stress-sensitive. CSDS did not affect locomotor activity of 129 or BALB mice, but both resilient and susceptible mice failed to gain weight during the defeat period. By contrast, the defeated D2 mice had lower locomotor activity than controls, and they lost weight during defeat. Differently to the other strains, the defeated B6 mice had lower locomotor activity compared to controls but they gained weight during defeat, similarly to controls. We also investigated other behaviors in this strain, and found that the stress-susceptible mice had increased anxiety-like behavior, but no difference in despair behavior, as in previous studies (Krishnan et al., 2007; Razzoli et al., 2011a). We acknowledge that due to limitations of our animal facility, all mice were brought into the experimental room together before behavioral testing, where they stayed behind a screen during the testing of other animals, possibly introducing confounding factors. Overall, our results suggest that the studied strains use different coping strategies to stress, and that such behavior has a strong genetic basis.
In addition to different behavioral responses to CSDS, we observed large differences in the brain transcriptomic response to stress between B6 and D2 strains. mPFC and BNST expression patterns were more similar within strains, between the resilient and susceptible mice, than between susceptible or resilient mice of different strains. However, in vHPC, the transcriptomic response of B6 and D2 susceptible mice was highly similar. This result may reflect the different organization and roles of these three brain regions in processing stress-related signals. It may be that the vHPC relays primary information regarding stress, while the BNST and mPFC may have more interpreting and processing roles, and therefore little genetic variation is tolerated within the hippocampal stress response. vHPC gene expression is influenced by glucocorticoids through glucocorticoid receptors, which are strongly expressed in the hippocampus (Ahima and Harlan, 1990; Ahima et al., 1991; Gray et al., 2017). BNST receives projections from several amygdalar nuclei and mediates anxiety-related information to hypothalamic and brainstem targets (Davis et al., 2010). The mPFC influences stress-associated social behavior by evaluation of perceived threats, top-down control of goal-directed behavior (Carlén, 2017), and emotion regulation (Etkin et al., 2015). These differences reflected on the transcriptomic response to stress may contribute to the distinct coping strategies of the strains.
We found significant enrichment of OLG- and aging-related genes within the differentially expressed genes in both strains and all brain regions. The highly significant overlap between both of these gene sets suggests similar changes in OLG-related gene expression may be associated with both aging and chronic psychosocial stress. This OLG-related signal was derived mostly from genes expressed in mature OLGs, which are the myelin-producing cells in the central nervous system. Myelin plasticity, a response of OLGs to neuronal activity (Gautier et al., 2015; Koudelka et al., 2016; Purger et al., 2016; Mitew et al., 2018), is retained in adulthood (Bengtsson et al., 2005; McKenzie et al., 2014). Downregulation of OLG-related genes was especially pronounced in the vHPC, and also in the mPFC of B6 susceptible mice. On the structural level, vHPC myelin sheaths of B6 susceptible mice were thinner than in control mice, concurring with the downregulation of myelin-related genes. Interestingly, in the BNST, B6 susceptible mice had thicker myelin sheaths and smaller g ratio compared to both resilient and control mice, also concurring with the gene expression patterns. D2 susceptible mice had an opposite gene expression pattern in the BNST, which, however, was not supported by changes in myelin sheath thickness or g ratio. Opposite gene expression changes in different inbred mouse strains in response to stress have previously been observed in other models, but the underlying mechanisms are not fully understood (Mozhui et al., 2010; Malki et al., 2015). In mice, various stressors have diverse effects on myelination, likely reflecting involvement of distinct neural processes, developmental stage at the time of stress exposure, and duration of stress. Early life stress affects myelin-related gene expression in the mPFC (Bordner et al., 2011; Makinodan et al., 2012), and myelination of the amygdala (Ono et al., 2008). In adult mice, social isolation (Liu et al., 2012) and 14-d CSDS in susceptible mice (Lehmann et al., 2017) reduce myelin-related gene expression and myelination within the frontal cortex. Intermittent social defeat stress leads to reduced MBP-stained myelin area in the mPFC (Zhang et al., 2016). Corpus callosum myelination is decreased after chronic restraint stress (Choi et al., 2017). Also chronic variable stress induces temporally variable myelin-related gene expression changes, and these changes were regionally selective in the mPFC, nucleus accumbens, and corpus callosum (Liu et al., 2018). Therefore, it is likely that alterations in gene expression and physical parameters of the myelin sheaths occur dynamically over days and weeks (Montesinos et al., 2015; Almeida and Lyons, 2017), possibly explaining thicker myelin in B6 resilient mice in the mPFC, without differences in myelin-related gene expression compared to controls. Overall, our results concur with earlier studies and suggest that stress does not simply cause widespread downregulation of myelin-related genes or myelin loss, but that stress-influenced myelin-plasticity involves specific stress-associated brain circuits.
Our strategy to divide the stressed animals to susceptible and resilient groups allowed identification of specific patterns of myelination-related differences between these groups. We observed that chronic stress associates with thicker myelin in the susceptible (in the BNST of B6 mice) and thinner myelin in resilient mice (D2 mPFC). No prior studies exist, to our knowledge, on resilience-related myelin thinning. Although DTI measures are only considered proxies for myelination, in humans, baseline fractional anisotropy has been demonstrated to correlate positively with state anxiety, but following exposure to a traumatic event, correlation is negative (Sekiguchi et al., 2014). Furthermore, increased fractional anisotropy, suggestive of enhanced integrity, has been reported in the cingulate gray matter of panic disorder patients (Han et al., 2008). Also, the number of mature OLGs in the ventromedial PFC of depressed subjects with childhood abuse is increased (Tanti et al., 2017). Our results suggest that these complex patterns are strongly modulated by genetic factors.
Although myelination studies on human anxiety disorders and chronic stress are scarce, major depression has been consistently associated with white matter disruptions (Wang et al., 2014). Depressed suicide completers with childhood abuse have impaired myelin-related gene expression and reduced myelin thickness in the anterior cingulate cortex (Lutz et al., 2017). Interestingly, this effect was only seen in small caliber axons, similarly to our finding of mPFC myelin thickness being larger in B6 resilient mice compared to controls. Axons of different diameter may represent specific types of neurons (Perge et al., 2012) and therefore be differentially vulnerable to the effects of stress. For example, in the macaque cortex, axons which project longer distance from the point of origin have larger diameters than those projecting to proximal targets (Innocenti et al., 2014). Based on the axon diameters measured by TEM, it is not possible to infer the type of neurons affected in our experiment but this question could be addressed with immuno-electron microscopy.
How may structural changes in myelin affect behavior? Myelin thickness is a key component in determining axonal conduction speed, thereby influencing circuit function. Myelin deficient rats have altered conduction time and lower average synchrony levels in the cerebellum (Lang and Rosenbluth, 2003). Mice have increased synchronous activity between the vHPC and mPFC in an anxiogenic environment (Adhikari et al., 2010). Furthermore, in primates, increased power and phase synchrony in the theta range has been detected in the amygdala-prefrontal circuit during aversive conditioning (Taub et al., 2018). Altered synchrony has also been recorded in human psychiatric disorders (Schulman et al., 2011; Leuchter et al., 2015). Genetic factors affect synchronous oscillations during sleep in the inbred mouse strains (Franken et al., 1998; Tafti et al., 2003). The role of brain oscillations in myelination and psychiatric phenotypes are still poorly understood, but myelin plasticity may be a mechanism that allows regulation of axonal conduction and spatial connectivity.
Overall, our unbiased brain gene expression analysis suggests that myelin plasticity is among the most significant responses to chronic psychosocial stress. Our strategy to divide the mice into stress-susceptible and resilient groups allowed us to demonstrate significant differences in myelination-related gene expression and myelin thickness between these groups. Furthermore, variable myelin plasticity across brain regions suggests that chronic stress has localized effects on myelination. Importantly, by using two different inbred mouse strains we demonstrated that stress-induced myelin plasticity is genetically controlled. Identification of the genetic regulators of the myelin response will provide mechanistic insight into the molecular basis of stress-induced anxiety, a critical step in developing targeted therapy for anxiety disorders.
Acknowledgments
Acknowledgments: We thank Saija-Anita Callan, Sanna Kängsepp, and Laura Salminen for help with mouse work; Maria Razzoli for help with setting up the CSDS protocol; Vootele Voikar from the Mouse Behavioral Phenotyping Facility [supported by University of Helsinki (HiLIFE) and Biocenter Finland] for support with behavioral testing; Jenni Lahtinen for RNA-seq library preparation; Juho Väänänen and Birgitta Paranko for help with bioinformatic analyses; Mervi Lindman, Arja Strandell and Helena Vihinen for help with electron microscopy; Päivi Laamanen, Harri Kangas and Matias Rantanen from the Institute of Biotechnology (University of Helsinki) for RNA-seq; the IT Center for Science (CSC) for computing facilities; and Petri Hyytiä, Anders Paetau, Eva Rikandi, Outi Mantere, Tuula Kieseppä, Jaana Suvisaari, Tuukka Raij, and Hovatta lab members for helpful discussions. We also thank the Sequencing Unit at FIMM Technology Center supported by University of Helsinki and Biocenter Finland.
Footnotes
The authors declare no competing financial interests.
This work was supported by the European Research Council Starting Grant GenAnx 281559 (to I.H.), ERA-NET NEURON (AnxBio; I.H.), University of Helsinki (I.H.), Sigrid Jusélius Foundation (I.H.), Doctoral Program Brain and Mind, University of Helsinki (M.L. and Z.M.), and the A*MIDEX Grant ANR-11-IDEX-0001-02 funded by the French Government “Investissements d’Avenir” Program (to P.A.). The EM-unit is supported by the Institute of Biotechnology (University of Helsinki), Helsinki Institute of Life Science and Biocenter Finland (E.J.).
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.