Abstract The molecular mechanisms that link stress and biological rhythms still remain unclear. The habenula (Hb) is a key brain region involved in regulating diverse types of emotion-related behaviours while the suprachiasmatic nucleus (SCN) is the body’s central clock. To investigate the effects of chronic social stress on transcription patterns, we performed gene expression analysis in the Hb and SCN of stress-naïve and stress-exposed mice. Our analysis revealed a large number of differentially expressed genes and enrichment of synaptic and cell signalling pathways between resilient and stress-naïve mice at zeitgeber 16 (ZT16) in both the Hb and SCN. This transcriptomic signature was nighttime-specific and observed only in stress-resilient mice. In contrast, there were relatively few differences between the stress-susceptible and stress-naïve groups across time points. Our results reinforce the functional link between circadian gene expression patterns and differential responses to stress, thereby highlighting the importance of temporal expression patterns in homoeostatic stress responses. Subject terms: Biological sciences, Neuroscience, Molecular neuroscience Introduction Mood disorders are complex heterogeneous disorders that affect cognition, emotional processing and basic physiological functioning [[40]1]. There is growing evidence of a tight connection between circadian rhythms and mood regulation [[41]2]. A majority of depressed patients show more severe symptoms in the morning than in the evening [[42]3]. Diurnal processes including diurnal rhythms in secretion of growth hormones and melatonin, daily core body temperature, and sleep/wake cycles are disrupted in some depressed patients [[43]4, [44]5]. These internal rhythmic processes are maintained on a 24-h cycle by multiple peripheral clocks throughout the body which are synchronised to a single central principal clock in the mammalian hypothalamus—the suprachiasmatic nucleus (SCN) [[45]6]. The SCN itself has been shown to affect depressive and anxious-like responses in mice [[46]7, [47]8]. Furthermore, the SCN projects to several hypothalamic structures that regulate diurnal processes disrupted in depression, including the paraventricular nucleus, subparaventricular zone, and dorsomedial hypothalamus [[48]4, [49]9]. Overall, direct and indirect SCN projections to mood-regulating brain areas are likely responsible for the diurnal rhythmicity of symptoms observed in depressed patients. One such mood-regulating brain area that is thought to receive inputs from the SCN is the habenular complex (Hb), which is composed predominantly of lateral (LHb) and medial habenula (MHb) subregions [[50]10]. The Hb relays information between the forebrain emotion-processing limbic system and the midbrain depression-associated aminergic centres [[51]11]. The Hb is heavily involved in mediating stress responses [[52]12, [53]13] and has also been associated with circadian modulation of depression [[54]14]. Elevated diurnal firing in LHb neurons projecting to the serotonergic dorsal raphe nucleus (DRN) is associated with social avoidance in chronically stressed mice [[55]15], while passive coping, a behaviour mirroring apathy and resignation in humans, has been associated with increased activity in numerous circuits, including: LHb-DRN projections [[56]16], medial prefrontal cortex (mPFC) projections to LHb [[57]17], and LHb projections to the ventral tegmental area (VTA) [[58]18]. The pharmacological inhibition of the LHb rescues depressive-like behaviour [[59]19, [60]20]. In contrast, MHb lesion impairs voluntary exercise and induces despair-like and anhedonic behavioural phenotypes in mice, without affecting the circadian period and sleep architecture [[61]21–[62]23]. Despite the apparent importance of Hb in mood regulation, the effects of stress on molecular pathways in the Hb remain underexplored. Both LHb and MHb exhibit an independent endogenous molecular clock that is synchronised with the SCN clock [[63]10, [64]24–[65]26] which drives the daily variation in their firing [[66]14, [67]24, [68]25]. Hb receives input from the SCN indirectly through the dorsomedial and lateral hypothalamic areas [[69]27–[70]29]. Some evidence suggests that SCN innervates the LHb directly [[71]14, [72]26, [73]30]. Additionally, the Hb complex encompasses the caudal perihabenular region (PHb), which receives direct projections from the retina and has been shown to drive circadian-gated light modulation of a subset of depressive-like behaviours in mice, independently of the SCN [[74]31–[75]34]. Although the PHb does not respond to four-day social defeat stress [[76]31], PHb-dependent light modulation of behaviour occurs through its projections to nucleus accumbens (NAc), mPFC and thalamic reticular nucleus – areas associated with mood, stress-coping and sleep regulation [[77]31, [78]34, [79]35]. Hence, the source of depression-associated diurnal rhythm disruption might be due to perturbations of circadian rhythms at the level of the SCN, the Hb or both. In order to investigate the relationship between mood regulation and diurnal rhythm disruption, we exposed mice to chronic social defeat stress (CSDS). CSDS has two important advantages over other models. Firstly, CSDS is a robust and widely used chronic stress model that induces depressive-like and PTSD-like symptoms in a subset of exposed rodents [[80]36, [81]37]. Experimental mice exhibiting social avoidance are deemed as ‘stress-susceptible’, while those that do not are ‘stress-resilient’ [[82]38]. These two subpopulations show different behavioural, physiological and molecular profiles, allowing for investigations of the underlying mechanisms leading to stress susceptibility or resilience [[83]36]. CSDS in rodents disrupts core body temperature, locomotor activity rhythms, sleep architecture, and other basic physiological functioning, with some biological processes being disrupted in all stress-exposed mice and others only in the stress-susceptible phenotype [[84]39–[85]43]. As a social form of stress, it has higher translational value for modelling human mood disorders [[86]37]. The other important reason for using the CSDS paradigm is that this model utilises a predictable stress exposure with standardised protocol and timing, unlike chronic mild stress (CMS). Hence, it is more suitable for investigation of time-dependent biological processes. Few studies have explored the effect of CMS on core clock gene expression [[87]44, [88]45], but the differences in core clock gene expression patterns in stress-resilient and stress-susceptible mice have not been investigated. CSDS was found to increase the amplitude of Per2 diurnal expression in the SCN [[89]43]. Moreover, LHb cells are responsive to intermittent social defeat [[90]31]. Specifically, LHb cells projecting to the DRN in stress-susceptible mice exhibit blunted diurnal activity compared to stress-resilient and stress-naïve mice [[91]15]. Despite these earlier findings suggesting a significant interaction between diurnal rhythms and chronic stress, the temporal transcriptomic profiles of the SCN and Hb in the stress-resilient and stress-susceptible mice remain unexplored. In this study, we investigated the effects of CSDS on the SCN and Hb transcriptome across the light/dark cycle. Mice exposed to CSDS were sampled at four time points: ZT4 (mid-inactive phase), ZT10 (transition to active phase), ZT16 (mid-active phase), and ZT22 (transition to inactive phase). We hypothesised that the greatest changes would be observed in the stress-susceptible mice in the Hb during the light phase. We expected relatively few differences between the resilient and stress-naïve mice, especially in the SCN. However, we found the largest transcriptomic changes to occur at ZT16 in resilient relative to stress-naïve mice, in both SCN and Hb. Our observations highlight the importance of establishing transcriptomic profiles across several time points to more accurately describe the homoeostatic changes in response to stress. We speculate that homoeostatic plasticity in the SCN may be important for modulating resilience to stress. Material and Methods Ethics approval All experimental procedures were conducted in accordance with the National Institute of Health Guide for Care and Use of Laboratory Animals and approved by the New York University Abu Dhabi Animal Care and Use Committee (IACUC Protocol: 19-0004). Animals and housing Male C57BL/6J mice (8–10 weeks old; Jackson Laboratory) were housed in groups of 5 upon arrival for at least a week before the start of the chronic social defeat stress (CSDS) paradigm. CD1 male mice (retired breeders; Charles River Laboratories) were housed individually under standard housing conditions. All mice were housed with ad libitum access to food and water at 20 ± 1 °C and 50 ± 10% humidity under a 12-h light/dark cycle (7:00 am / ZT0: lights on; 7:00 pm /ZT12: lights off). Chronic social defeat stress CSDS was performed as described previously [[92]38]. Experimental C57 mice were exposed to 10 min of physical stress followed by 24 h of sensory stress in the cage of an unfamiliar CD1 aggressor mouse (Fig. [93]S1A). CSDS was repeated daily for 15 days between ZT8 and ZT9 with a novel CD1. The timing of the daily defeat sessions was selected to complement a long record of published studies that conducted CSDS during the light phase [[94]36, [95]46–[96]48], but failed to report or follow consistent stress-exposure timing across days. The ZT8-ZT9 was specifically selected since mice exhibit low activity levels [[97]15], as well as lower levels of sleep pressure [[98]49, [99]50]. Control mice were housed in pairs in a separate room and moved each day at the time of defeat session. All the mice were single-housed following the end of CSDS until sacrifice. Social avoidance towards an unfamiliar non-aggressive CD1 mouse was measured in a two-stage social interaction (SI) test, as previously reported [[100]38]. The test was conducted between ZT7 and ZT10. TopScan video tracking system (CleverSys. Inc.) was used to automatically track mouse movement. Tissue sampling and region isolation Mice were anaesthetised using isoflurane and sacrificed 48–72 h after the SI test by cervical dislocation. Sampling was done at four time points (ZT4, ZT10, ZT16, ZT22) across the light-dark cycle (Fig. [101]1A, B). Extractions were carried out on a cold sterile surface, and the extracted brains were flash-frozen by immersion into −60 °C isopentane (Sigma Aldrich M3263) for 20 s. The brains were kept on dry ice at -80°C until further processing. The flash-frozen brain samples were then sliced on a Leica CM1950 cryostat (250 µm thickness) (Fig. [102]S1B). The regions were isolated with reference to the Mouse Brain Atlas [[103]51]. For the SCN and Hb, 2–3 bilateral punches per region per biological replicate were combined in a single sample (Fig. [104]S1B). There was no pooling across biological replicates. Fig. 1. Experimental design. [105]Fig. 1 [106]Open in a new tab A Timeline for investigating the effects of stress on temporal changes in gene expression in stress-naïve and stress-exposed animals. B Sampling time points and emerging comparison levels from pooling the sampling time points. RNA isolation, library preparation and RNA-sequencing RNA was isolated from the punches using the Qiagen RNeasy Micro Kit, as per manufacturer’s instructions. Isolated RNA was quantified and checked for quality on an Agilent Bioanalyzer instrument (Agilent Technologies) (Fig. [107]S1C). Samples with RIN values above 7 were used for the library preparation. All the libraries were prepared together using the NEBNext® Ultra™ II RNA Library Prep Kit (New England BioLabs Inc.) as per the manufacturer’s protocol within one sitting. Sequencing was done on NovaSeq S1 flow cell (200 cycles) following Illumina’s recommendations. Differential expression and pathway enrichment analysis Count preprocessing Raw FASTQ sequenced reads were first assessed for quality using FastQC (v0.11.5) [[108]52]. Sequencing depth did not meaningfully differ between the phenotypes and timepoints, averaging 15.79 and 15.56 of TPM-normalised counts for Hb and SCN respectively (Supplementary Table [109]1). The reads were then passed through Trimmomatic (v0.36) for quality trimming [[110]53] and processed with Fastp to remove poly-G tails and Novaseq-specific artefacts [[111]54]. The reads were aligned to the mouse reference genome GRCm38.82 using HISAT2 with the default parameters [[112]55]. The resulting SAM alignments were then converted to BAM format and coordinate-sorted using SAMtools (v1.3.1) [[113]56]. Raw counts were generated by passing the sorted alignment files through HTSeq-count (v0.6.1p1) [[114]57]. Concurrently, the sorted alignments were processed in Stringtie (v1.3.0) for transcriptome quantification [[115]58]. Finally, Qualimap (v2.2.2) was used to generate RNAseq-specific QC metrics for each sample [[116]59]. All of the methods were executed using predefined YAML workflows in the BioSAILs workflow management system [[117]60]. Differential expression Differential expression testing was performed in R package DESeq2 (v1.40.2) [[118]61]. The p-values were adjusted using the Benjamini–Hochberg procedure. To obtain a better estimate of expression changes in low count genes, we performed log[2]-fold-change (LFC) shrinkage using the apeglm method [[119]62]. We used an adjusted p-value of 0.05 (5% FDR) and a log[2]-fold-change (LFC) of 0.5 to identify differentially expressed genes (DEGs). Samples within the same behavioural phenotype from ZT4 and ZT10 time points were computationally pooled together to provide an estimate of the ‘daytime’ expression, while ZT16 and ZT22 time point samples were pooled to provide an estimate of the ‘nighttime’ expression for each phenotype (Fig. [120]1B). Pooling all the samples within the same behavioural phenotype provided an estimate of time-invariant ‘global’ expression corresponding to average expression throughout the day. Differential expression between the behavioural phenotypes was carried out at each of these temporal resolution levels. SCN and Hb samples were assessed for presence of respective canonical markers from published scRNA-seq studies [[121]63–[122]65]. The sequenced SCN samples specifically expressed known SCN markers, while Hb strongly expressed both the LHb and MHb markers, suggesting sampling accuracy (Fig. [123]S2). However, we did not manage to reliably sample the PHb (Fig. [124]S2A). In order to assess whether technical biases were driving the observed differential expression patterns, the raw counts were normalised and stabilised for variance using DESeq2 with a design formula including region, time point and phenotype as covariates. The output was then used to compute sample similarity across all the genes (Fig. [125]S3A), as well as the top 2000 most variable genes (Fig. [126]S3B), which were then clustered based on their Eucledian distance. Beyond the sample region clustering, these heatmaps do not suggest any other specific grouping that would indicate a bias. Furthermore, we performed principle component analysis (PCA) on all the samples (data not shown). These results were replicated using surrogate variable analysis (SVA) (data not shown) [[127]66]. First principal component represented the regional identity of the samples and accounted for 82% of data variability (Fig. [128]S2B). Pathway enrichment Enrichment analysis was conducted in the Ingenuity Pathway Analysis (IPA, v.23.0, Qiagen) on all genes that had passed the 5% FDR significance criterion and was conducted only in those comparisons where the number of DEGs was higher than 20. Pathways were considered significantly enriched if –log(B-H p-value) was greater than 1.3. Based on the known roles of different genes and their relationship in the functionally-predefined pathways, IPA calculates z score which represents the predicted direction of change of the pathway (activation: z > 0; inhibition: z < 0). The analysis included only the pathways associated with the mouse nervous system, neurons, astrocytes, other tissues and unspecified cell types within the IPA framework. Gene ontology We visualised the associated DEGs and the top five IPA pathways from the resilient-to-control comparisons at ZT16 for both regions using gene-concept networks. Gene-concept networks represent the relative coordinates of the enriched pathways based on the number and differential expression of the associated DEGs. Therefore, it is possible to identify overlapping DEGs that have the largest contribution to the enrichment of the top five activated pathways. In silico TF motif enrichment To explore which transcription factors (TF) might be driving the differential expression at specific time points, we used FIMO [[129]67] (fimo --oc {gene.fa} --verbosity 1 --bgfile --nrdb-- --thresh 1.0E-4 H12CORE_meme_format.meme {gene}) with the Mouse TF motif list from HOCOMOCO [[130]68] to generate TF binding motif enrichment within 20 kb upstream of DEG promoters. The FASTA sequences were retrieved using BED2FASTA [[131]67]. FIMO results were then filtered for hits corresponding to q ≤ 0.05 and for TF that were expressed above the median of normalised and transformed expression value. TF hits were then sum-aggregated by TF families. Given the unbalanced number of DEGs by conditions, aggregated TF hits were normalised by the number of hit genes per condition. Finally, the delta (Δ) between upregulated and downregulated hits was computed for each TF family. Positive delta values correspond to enriched binding of the TF family in the promoter of upregulated DEGs, and conversely negative values correspond to binding in the promoter of downregulated DEGs. Detection of rhythmic genes Log-normalised count values were used to determine the rhythmic expression profile of core clock genes across all three behavioural phenotypes using the non-parametric JTK method for identification of expression oscillations in genome-scale datasets [[132]69]. JTK method was run in the R-implemented version of the ECHO (v.4.0) tool [[133]70]. Results CSDS induces social avoidance in a subset of stress-exposed mice To investigate the effects of chronic stress on SCN and Hb transcriptome, we exposed mice to 15 days of CSDS. Following CSDS, mice were classified as stress-resilient or stress-susceptible using SI test (Fig. [134]1A, [135]S4A–C). Stress-naïve and stress-resilient mice exhibited significantly larger social interaction values than the stress-susceptible mice (Fig. [136]S4A–C). Resilient and susceptible mice exhibited decreased movement and velocity relative to the stress-naïve mice (Fig. [137]S4D, E). Stress-resiliency is associated with large gene upregulation in the SCN and Hb at ZT16 In order to investigate differential expression between the behavioural phenotypes at different time points, we sampled SCN and Hb at ZT4, ZT10, ZT16 and ZT22 time points 48–72 h following the SI test (Fig. [138]1A, B). Computationally merging ZT4 and ZT10 provided a daytime expression estimate, while merging ZT16 and ZT22 provided nighttime expression estimate (Fig. [139]1B). Combining all the time points together provided computational estimate of time-independent expression patterns (Fig. [140]1B). Comparing the phenotypes revealed that, in both the SCN and Hb, the highest number of significant genes was observed between resilient and stress-naïve mice at all temporal resolution comparison levels (Supplementary Tables [141]2, [142]3). In the SCN (Fig. [143]2A, B), at global (average gene expression across all 4 time points) and daytime (average gene expression at ZT4 and ZT10) comparisons only a few genes were differentially expressed between any of the groups. Comparison of individual gene expression patterns at ZT4 and ZT10 time points did not reveal substantial differences between the phenotypes. However, at nighttime (average gene expression at ZT16 and ZT22) comparisons, we found 338 DEGs in resilient relative to stress-naïve mice, of which a large majority (~80%) were upregulated (Fig. [144]2A, C). The nighttime differences in the SCN predominantly arose from the differential expression observed at ZT16 where we found 758 DEGs in resilient relative to stress-naïve mice, of which ~60% were upregulated (Fig. [145]2B, E). At ZT22, we observed 34 upregulated DEGs in resilient relative to susceptible mice. However, other phenotypic comparisons at ZT22 did not reveal any major differential expression patterns. These results suggest that the largest changes in stress-induced SCN gene expression occur in resilient mice during nighttime at ZT16. Fig. 2. Differentially expressed genes (DEGs) (FDR < 0.05 and |log[2]FC| > 0.5) between the phenotypes at different time points in the SCN. [146]Fig. 2 [147]Open in a new tab A Number of DEGs at pooled time points. At the global comparison level, there were 16, 5 and 7 DEGs in resilient-to-control, susceptible-to-control and resilient-to-susceptible comparisons, respectively. Daytime comparison revealed 2, 3 and 1 DEGs for resilient-to-control, susceptible-to-control and susceptible-to-resilient comparisons, respectively. Nighttime comparison revealed 338, 11 and 10 DEGs for resilient-to-control, susceptible-to-control and susceptible-to-resilient comparisons, respectively. B Number of DEGs at specific time points. At ZT4, there were 0, 3 and 7 DEGs for resilient-to-control, susceptible-to-control and resilient-to-susceptible comparisons. At ZT10, there were no DEGs for any phenotypic comparison. At ZT16, there were 758, 17 and 23 DEGs for resilient-to-control, susceptible-to-control and resilient-to-susceptible comparisons, respectively. At ZT22, there were 0, 0 and 34 DEGs for resilient-to-control, susceptible-to-control and resilient-to-susceptible, respectively. C Nighttime comparison: Volcano plot showing significant DEGs (FDR < 0.05) in resilient mice relative to controls at nighttime comparison level in SCN. D ZT16 comparison: Volcano plot showing significant DEGs (FDR < 0.05) in resilient mice relative to controls at ZT16 in SCN. Red and blue points represent significantly upregulated and downregulated genes, respectively. In the Hb (Fig. [148]3A, B), differential expression patterns between the phenotypes were somewhat similar to those observed in the SCN. Global comparisons in the Hb revealed 68 upregulated DEGs in the resilient relative to stress-naïve mice (Fig. [149]3A, C). Gene expression patterns did not differ between any of the three behavioural phenotypes during the daytime. Specifically, ZT4 and ZT10 time point comparisons did not yield any discernible differences in gene expression patterns (Fig. [150]3B). However, nighttime comparison revealed 55 upregulated DEGs in resilient relative to stress-naïve mice (Fig. [151]3A, D). We found that the large differences in global and nighttime gene expression profiles are again mainly due to the large numbers of upregulated genes at ZT16 in stress-resilient relative to control and susceptible mice (Fig. [152]3A, B). Specifically, 143 DEGs were identified in the resilient relative to stress-naïve mice at ZT16, with almost all genes being upregulated (Fig. [153]3B, E). We observed 78 upregulated DEGs in resilient mice compared to stress-susceptible mice in the Hb, unlike in the SCN (Fig. [154]3B). Only small differences were observed between the susceptible and control mice at ZT16 in the Hb. No robust differences were observed at ZT22 between any of the behavioural phenotypes in the Hb (Fig. [155]3B). Hence, stress-induced differences in Hb gene expression patterns in the resilient mice, both at global and nighttime analysis, are driven by changes at the ZT16 time point, similar to the SCN. Fig. 3. Differentially expressed genes (DEGs) (FDR < 0.05 and |log[2]FC| > 0.5) between the phenotypes at different time points in the Hb. [156]Fig. 3 [157]Open in a new tab A Number of DEGs at pooled time points. At the global comparison level, there were 69, 0 and 1 DEGs for resilient-to-control, susceptible-to-control and resilient-to-susceptible comparisons, respectively. At the daytime comparison level, there were 0, 0 and 2 DEGs for resilient-to-control, susceptible-to-control and susceptible-to-resilient comparisons, respectively. At the nighttime comparison level, there were 55, 0 and 0 DEGs for resilient-to-control, susceptible-to-control and susceptible-to-resilient comparisons, respectively. B Number of DEGs at specific sampled time points. At ZT4, there were 1, 0 and 0 DEGs for resilient-to-control, susceptible-to-control and resilient-to-susceptible comparisons. At ZT10, there were no DEGs for resilient-to-control and resilient-to-susceptible comparisons, and only 1 DEG in susceptible-to-control comparison. At ZT16, there were 143, 2 and 80 DEGs for resilient-to-control, susceptible-to-control and resilient-to-susceptible comparisons, respectively. At ZT22, there were no DEGs in any of the phenotypic comparisons. C Global comparison: Volcano plot showing significant DEGs (FDR < 0.05) in resilient mice relative to controls at global comparison level in Hb. DEGs were labelled based on |log2FC| > 1.0 and biological relevance. D Nighttime comparison: Volcano plot showing significant DEGs (FDR < 0.05) in resilient mice relative to controls at nighttime comparison level in Hb. E ZT16 comparison: Volcano plot showing significant DEGs (FDR < 0.05) in resilient mice relative to controls at ZT16 time point in Hb. Red and blue points represent significantly upregulated and downregulated genes, respectively. Given the strong differential expression at ZT16 between resilient and stress-naïve mice, we examined in sillico possible TF candidates among DEGs that could drive these differences (Fig. [158]S5). We found a greater number of possible TFs associated with downregulated than upregulated DEGs in both SCN and Hb. The two brain regions showed enrichment of 35 common TF families. Differential expression in the SCN and Hb was associated with 8 and 39 region-specific TF families, respectively. Overall, these results show robust transcriptomic differences in resilient relative to stress-naïve mice in both the SCN and Hb. The largest observed differences occur at ZT16 in both regions, with SCN showing larger transcriptional response than the Hb in terms of DEG numbers and associated LFC. Overall, susceptible mice did not show significant differential expression patterns in most comparisons relative to either the resilient or control mice. Taken together, the results suggest resiliency-associated transcriptional signature to be temporally-restricted to the middle of the active phase. Stress resilience is associated with enrichment of signalling and plasticity pathways in SCN at ZT16 To better understand the molecular pathways and functions of DEGs between the resilient and stress-naïve mice, we performed pathway enrichment analysis in the SCN at nighttime and at ZT16 and found a number of enriched pathways (Supplementary Table [159]4, Fig. [160]4). At ZT16, the enriched pathways included an expanded list of those observed from the nighttime comparison (Supplementary Table [161]4, Fig. [162]4A, B). All the enriched pathways in the resilient mice can be divided into two broad functional groups: signalling and plasticity-related pathways (Fig. [163]4A). Furthermore, at ZT16, ‘circadian rhythm signalling’ pathway was significantly enriched in the resilient mice. Both signalling and plasticity-related pathways in the resilient mice were enriched by common DEGs (subset shown in Supplementary Table [164]5). Many upregulated DEGs were associated with most enriched pathways (Fig. [165]4C). These results suggest that stress resilience is strongly associated with upregulation of relatively few genes in SCN at ZT16 which subsequently have broad impact on signalling and plasticity regulation, possibly through their role as TFs. Fig. 4. Enriched pathways in the SCN of resilient mice relative to controls. [166]Fig. 4 [167]Open in a new tab The activation pattern of enriched biological pathways in the SCN. All lists are sorted by adjusted p-value. Red denotes activation (z > 0), blue denotes inhibition (z < 0), white denotes no net effect (z = 0), and grey denotes uncertain activation pattern. A Enriched signalling pathways. B Enriched plasticity-associated pathways. C Gene-concept networks of the 5 most enriched IPA pathways at ZT16 between resilient and control phenotypes in the SCN. Stress resilience is associated with enrichment of synaptogenesis and calcium signalling pathways in the Hb at ZT16 In the Hb, we explored pathway enrichment at global, nighttime and ZT16 comparison levels between the resilient and stress-naïve mice (Fig. [168]5). Fewer pathways were as highly enriched and activated in the Hb of the resilient mice at ZT16 compared to the SCN (Supplementary Table [169]6, Fig. [170]5A, B). The most enriched pathways in this comparison were associated with mechanisms of synaptic plasticity (Supplementary Table [171]6, Fig. [172]5B). Similar to the SCN, ‘circadian rhythm signalling’ pathway was enriched in the resilient mice at global, nighttime and ZT16 comparisons. Relatively few DEGs were associated with a high number of enriched pathways in the Hb of resilient mice (Supplementary Table [173]5). Furthermore, the most enriched pathways were associated with a handful of genes (Fig. [174]5C). These results suggest that stress resilience is associated with upregulation of plasticity-associated genes in the Hb at ZT16. Although both the SCN and Hb exhibited robust changes specifically at ZT16, the SCN shows larger and functionally more diverse transcriptomic response to stress than the Hb. Fig. 5. Enriched pathways in the Hb of resilient mice relative to controls. [175]Fig. 5 [176]Open in a new tab The activation pattern of enriched biological pathways in the Hb. All lists are sorted by adjusted p-value. Red denotes activation (z > 0), blue denotes inhibition (z < 0), white denotes no net effect (z = 0), and grey denotes uncertain activation pattern. A Enriched signalling pathways. B Enriched plasticity-associated pathways. C Gene-concept networks of 5 most enriched IPA pathways at ZT16 between resilient and control phenotypes in the Hb. Stress induces arrhythmicity of Per1 and Arntl in the SCN of stress exposed mice Given the enrichment of ‘circadian rhythms signalling’ pathway, we explored whether chronic stress affects the rhythmicity of core clock genes in the SCN and Hb. To this end, we performed JTK analysis and found rhythmic expression of the core clock genes (Per1, Per2, Per3, Cry1, and Arntl) in the SCN (p[adj] <0.05) of stress-naïve mice (Fig. [177]6A, B, Supplementary Table [178]7). Clock was arrhythmic in both brain regions across all the phenotypes. In the SCN, Per1 was arrhythmic specifically in the stress-susceptible mice. Arntl was arrhythmic in both stress-resilient and stress-susceptible mice. In contrast Cry2 was rhythmic only in the SCN of resilient mice. In the Hb, the expression of core clock genes was more uniform across the phenotypes. With the exception of Cry1 and Clock, all genes were rhythmic in stress-resilient and stress-susceptible mice. Moreover, Per1, Per3 and Cry1 were arrhythmic in stress-naïve mice, but rhythmic in resilient and susceptible mice (Fig. [179]6A, B). Beyond core clock genes, large differences in numbers of phenotype-specific rhythmic genes were observed in both the SCN and Hb (Fig. [180]S6). Fig. 6. Clock gene rhythmicity. [181]Fig. 6 [182]Open in a new tab A Temporal expression of core clock genes in the SCN and Hb of control, resilient and susceptible mice. All the clock genes show previously reported temporal expression in the SCN of control mice. B Dotplot showing rhythmic expression of clock genes in SCN and Hb. SCN: Per1 is arhythmic stress-susceptible mice. Arntl is arrhythmic in both resilient and susceptible mice. Cry2 is rhythmic in resilient mice. Hb: Per1, Per3 and Cry1 are arrhythmic in the stress-naïve mice, while being rhythmic in both the resilient and susceptible mice. Discussion Stress resilience and homoeostatic response The close relationship between diurnal rhythms and mood disorders warranted a systemic investigation of transcriptomic changes in the principal clock (SCN) and mood-regulating brain regions. We focused on the Hb because it is a semi-autonomous oscillator with a strong role in stress regulation and mood disorders [[183]14]. We found that chronic stress elicited a robust change in the gene expression in the SCN and to a lesser degree in the Hb (Figs. [184]2–[185]3). This transcriptomic signature was specific to resilient mice and was strongest at ZT16. In contrast, comparing susceptible and stress-naïve mice did not reveal strong differential expression patterns. Large differential upregulation in resilient mice was associated with enrichment and predicted activation of signalling and plasticity-related pathways in both regions (Figs. [186]4–[187]5). By computationally pooling the samples and decreasing the temporal resolution, we demonstrated that such strong transcriptomic signatures are lost. Resilience is an active homoeostatic process and involves transcriptional changes in a diffuse network of interconnected brain regions [[188]71]. The current understanding is that chronic stress imposes an adaptive challenge that requires behavioural flexibility and physiological homoeostatic buffering [[189]72]. The inability to elicit such a response leads to failure in stress-coping and strongly contributes to associated stress-induced pathologies [[190]72, [191]73]. Resilient mice do not show behavioural deficits likely because of the large transcriptomic changes such as those observed in our study. Pioneering studies found that resilient mice exhibit large and unique upregulation of neurotrophic factors and stress-responsive genes in the prefrontal cortex [[192]36]. Others found strong gene upregulation in resilient mice relative to stress-naïve and susceptible mice in midbrain regions such as the VTA and NAc [[193]74, [194]75]. The physiological consequences of large transcriptomic changes in resilient mice are also likely responsible for observed differences in circuit-specific neural activity between resilient and susceptible mice. For example, VTA dopaminergic (DA) cells in susceptible, but not resilient, mice exhibit elevated firing rate [[195]76]. The pathophysiological elevated firing is counteracted in resilient mice by compensatory mechanisms that are absent in susceptible mice [[196]77]. We speculate that the relatively negligable differences between resilient and susceptible mice were observed because chronic stress induces similar transcriptional changes in both phenotypes. However, larger differential expression in resilient mice would be sufficient in driving necessary homoeostatic response for acquisition of the resilient phenotype. Susceptible mice do not adapt to the same degree as resilient mice and are transcriptomically found in an ‘intermediate’ state [[197]36, [198]77], explaining negligible differences in relation to both resilient and stress-naïve mice. Our findings that large transcriptomic changes occur in resilient mice at a specific night time point (ZT16) add a significant temporal dimension to the current understanding of the molecular mechanisms that drive stress-resilience. Our data, along with previously published physiological and circuit studies, suggest that temporal patterns of regional activity may play a crucial role in driving a particular behavioural phenotype [[199]15]. It provides further evidence of a tight association between large transcriptional changes and successful stress-coping. Future studies with sampling performed over several day/night cycles following chronic stress would allow for temporal stability estimation of this transcriptomic signature. SCN transcriptomic profile We observed significant upregulation of Adora2a in the SCN of resilient mice at ZT16. Adenosine is known to regulate Per1 and Per2 expression via A1/A2A receptor activation which is hypothesised to integrate light and sleep signalling to track homoeostatic sleep pressure and regulate sleep architecture - both of which were previously found to be differentially impacted in resilient and susceptible mice [[200]39, [201]41, [202]78]. Upregulation of Adora2a in the SCN of resilient mice at ZT16 may be a subcomponent of adaptive homoeostatic response that maintains consolidated sleep [[203]39, [204]41]. Interestingly, the upregulation of A2A receptors in the forebrain is associated with depressive-like behaviour and decreased synaptic plasticity following unpredictable CMS [[205]79]. Resilient mice also displayed strong upregulation of Ntrk1 and Ngfr in the SCN at ZT16. NGF binds to Ntrk1 and Ngfr [[206]80, [207]81]. NGF levels are significantly lower in depressed patients [[208]82] while increasing its levels has antidepressant properties [[209]83, [210]84]. Furthermore, BDNF, which has functional roles in synaptic plasticity, stress response, and antidepressant action in several regions, also binds to Ngfr [[211]85, [212]86]. Ngfr is specifically synthesised in cholinergic neurons which strongly modulate SCN neuronal excitability and circadian rhythmicity at night [[213]87–[214]89]. We also observed high upregulation of cholinergic muscarinic receptors Chrm 1-5. These receptors have a functional role in fear and emotional responsiveness [[215]90, [216]91]. Cholinergic signalling in the SCN can induce phase shifts in locomotor rhythms [[217]92]. However, at present it is unclear how increased muscarinic receptor expression at ZT16 in the SCN drives resilience. A recent study demonstrated that the stress-associated neuropeptide CRH increases cholinergic interneuron firing that causes potentiation of dopamine transmission [[218]93, [219]94]. Besides cholinergic receptors, Drd2 and Crh were also upregulated in the SCN of resilient mice at ZT16. Dopaminergic input from VTA to SCN modulates circadian entrainment through activation of excitatory Drd1 receptors [[220]95]. The putative homoeostatic role of inhibitory Drd2 in circadian regulation in resilient mice remains unknown. Our observation that Crh is upregulated in resilient mice adds to the recent findings that CRH-expressing neurons are also present in the SCN [[221]96]. CRH increases gating of BDNF signalling [[222]97]. Thus, increased CRH expression at ZT16 may be responsible for putative BDNF-induced antidepressant-like effects in stress-resilient mice. The resilient phenotype has also been associated with increased cAMP and GPCR-mediated signalling [[223]74, [224]98]. These plasticity-related pathways were also enriched in our dataset and predicted to be activated in the SCN of resilient mice. However, the functional consequences of these changes are unclear since plasticity-related transcription factor Creb3l1 was downregulated, while many others, such as Gbx1, were strongly upregulated in the resilient mice. Previously, however, Creb downregulation has been associated with stress resilience [[225]99, [226]100]. Overall, few studies have examined stress effects on SCN physiology and those predominantly focused on the expression of clock genes [[227]7, [228]101–[229]103]. Our observations highlight potentially novel molecular mechanisms in the SCN that promote resilience to stress, warranting further investigation. Hb transcriptomic profile Adcy1 and Adcy9 were upregulated in the Hb of resilient relative to stress-naïve mice at ZT16. Adenylyl cyclases play a pivotal role in both cAMP-mediated and the ERK½-CREB signalling cascade which are known to modulate neuronal physiology [[230]104–[231]106]. Overexpression of Adcy1 in the forebrain prevents BDNF downregulation in stressed mice, thereby promoting stress resiliency [[232]107]. Adcy1 forebrain overexpression also induces stronger long-term potentiation (LTP) [[233]106]. We also observed ZT16-specific upregulation of other plasticity-modulating genes in resilient mice, such as Grin2a, Grin2b, Cacna1h, Camk2b, and Ppp3ca which correlate with large predicted activation of ‘synaptogenesis’ and ‘LTP signalling’ pathways. Resilient mice were found to exhibit larger hippocampal LTP relative to susceptible mice because of higher expression of Grin2b and Camk2 [[234]108]. Upregulation of plasticity-related genes may be a part of a larger homoeostatic response associated with stress resiliency [[235]109, [236]110]. We also observed upregulation of Gad1, Gabra5, Kcnq3 and Npy1r genes in stress-resilient mice. Sub-chronic mild stress decreases long-term depression (LTD) of GABAergic neurons in the LHb [[237]111]. Furthermore, male mice with Kcnq3 mutation that inhibits GABA binding exhibit reduced self-grooming, behaviour indicative of apathy [[238]37, [239]112]. Neuropeptide Y action through its Y1 receptor has been associated with stress resilience and anxiolysis in several brain regions, including the Hb [[240]113, [241]114]. In the LHb, Npy1r-expressing neurons play a critical role in chronic stress adaptation by driving palatable food intake and reducing anxious-like behaviour [[242]115]. Thus, it is possible that elevated Npy1r expression at ZT16 in the Hb drives increased food intake in resilient mice to compensate for lower levels of adipose tissue and serum leptin [[243]116]. Our data suggests that Hb might play an important role in promoting stress-resilience, beyond its known role in promoting depressive-like phenotype. Clock genes Since our analysis revealed a strong effect of stress on transcriptional profiles in the SCN and Hb at a specific time point, we examined the rhythmic expression of clock genes in both of these regions. We found that the resilient and susceptible mice exhibit similar rhythmic expression profile of the core clock genes in the SCN with the exception of Per1 and Cry2 which were differentially modulated by CSDS. Chronic stress induced arrhythmicity of Arntl in the SCN of both resilient and susceptible mice which might imply that putative homoeostatic adaptations in resilient mice do not protect against stress-induced disruption in rhythmic expression of this specific clock gene. Despite previous observation that knockdown of Arntl leads to despairful-, anxious- and helpless-like behaviour following CMS, our results suggest that Arntl arrhythmicity in the SCN does not directly impact mood regulation [[244]7]. Interestingly, Cry2 was rhythmically expressed in the stress-resilient mice, unlike stress-naïve and stress-susceptible mice, potentially pointing to a compensatory mechanism involving induced rhythmic expression of Cry2. Using in vitro bioluminescence methods, one study found that Per2 expression amplitude was increased in the SCN of mice that underwent nighttime CSDS, but not daytime CSDS [[245]43]. Although inconclusive due to low temporal resolution, our data also suggests amplitude blunting. Others have shown that exposure to CSDS during night induces differences in locomotor activity and tissue-specific Per2 expression profile [[246]117, [247]118], but also that resting-phase CSDS more strongly supresses locomotor activity [[248]101]. In contrast no difference in Per2 expression was observed in mice kept in constant dim red light where clock rhythms are free running [[249]119]. These observations imply that phase-specific CSDS exposure modulates rhythmic behavioural and gene expression profiles. Future studies will investigate whether the transcriptomic profile of mice exposed to chronic stress at night induces similar transcriptomic profile to the one observed in resilient mice in our study. Finally, rhythmic expression profile of SCN of stress-naïve mice in our study match previously published, more spatially-resolved findings [[250]64]. In our study, Per1, Per3, and Cry1 were arrhythmic in the Hb of stress-naïve mice, but not resilient and susceptible mice. Although unexpected, these results can be attributed to the subdivisional differences at the level of MHb and LHb [[251]10, [252]24–[253]26]. Clock gene expression in the MHb has yet to be systematically investigated in mice, though it has been reported previously in hamsters and rats [[254]120, [255]121]. Despite both the LHb and MHb showing diurnal variation in firing [[256]24, [257]25], the LHb displays more robust rhythms with higher amplitude than the MHb [[258]10]. In the present study we did not differentiate between the two subregions and the reduced spatial resolution might have contributed to the observed arrhythmicity of clock genes in stress-naïve mice. Moreover, different stressors have been previously shown to differentially impact neuronal dynamics [[259]37]. Thus, the differential effect of CSDS and CMS on the transcriptomic landscape may be explained by the utilised stressor type. Future studies using higher sampling resolution would be required to fully describe effects on the phase of core clock gene expression. Limitations There are several limitations to our study. An increasing number of studies report stronger perturbations following nighttime, rather than daytime, defeat [[260]42, [261]43, [262]117, [263]118, [264]122]. Our daytime defeat may miss potential molecular differences between resilient and susceptible mice that might be evident in mice exposed to nighttime stress. Furthermore, bulk RNA-seq approach offers limited spatial resolution where neither the subregions nor neuronal subpopulations can be clearly differentiated. It has been previously shown that cell types and subpopulations can widely differ in rhythmic profiles of core clock genes [[265]64]. Even though PHb does not respond to short-duration CSDS [[266]31], it was not reliably sampled in our study. Lastly, sampling at additional time points to increase the temporal resolution would have allowed for more robust temporal gene expression profiling. Future studies will elucidate subregion-specific effects using techniques of greater spatial resolution, such as scRNA-seq, and with larger number of time points for increased temporal resolution. Concluding remarks The SCN and Hb are heterogeneous regions with many functionally and molecularly distinct subpopulations. Previous transcriptomic studies have only examined specific subpopulations of the Hb for purposes of circuit manipulations [[267]18, [268]35, [269]123]. Likewise, studies linking SCN and stress were limited to investigating a limited predetermined list of genes [[270]7, [271]43]. To the best of our knowledge, no study has systematically examined the SCN or Hb bulk transcriptome in the context of murine chronic stress models across multiple time points throughout the diurnal cycle. Our study is the first to show that the largest transcriptomic changes in the SCN and Hb of resilient mice occur at night with a peak at specific time point. This highlights the importance of sampling over a 24-hour cycle to accurately determine dynamic changes in gene expression following stress exposure, since failing to do so conceals strong transcriptomic signatures as demonstrated with our computational pooling of day and night time points. Most studies fail to consider diurnal changes in expression profiles. Our time-dependent transcriptomic analysis suggests novel potential molecular mechanisms for driving resilience to stress that would not be evident with daytime-limited sampling. Such findings further support time-dependent medication treatment administration, which might enhance its efficacy. Crucially, our results highlight the link between temporal expression patterns and homoeostatic stress responses that remains to be investigated further. Supplementary information [272]Supplemental Figures Titles^ (20.8KB, docx) [273]Supplemental Figure 1^ (439.9KB, pdf) [274]Supplemental Figure 2^ (1.3MB, pdf) [275]Supplemental Figure 3^ (338KB, pdf) [276]Supplemental Figure 4^ (78.8KB, pdf) [277]Supplemental Figure 5^ (183.8KB, pdf) [278]Supplemental Figure 6^ (92.4KB, pdf) [279]Supplemental Table 1^ (22KB, docx) [280]Supplemental Table 2^ (24.1KB, docx) [281]Supplemental Table 3^ (24.1KB, docx) [282]Supplemental Table 4^ (21.7KB, docx) [283]Supplemental Table 5^ (20.4KB, docx) [284]Supplemental Table 6^ (20.8KB, docx) [285]Supplemental Table 7^ (24.8KB, docx) Acknowledgements