Distinct mechanisms of transcriptomic habituation to repeated stress in the mouse hippocampus

Domain: biology
Published: 2025-12-22
Editorial outcome: minor_revision

Abstract
Chronic stress is a risk factor for neuropsychiatric disorders, making the ability to adapt to repeated stress a crucial determinant of mental health. It remains unclear whether repeated stress exposure leads to habituation or new adaptive responses. Here, we investigated how the transcriptomic response to acute restraint stress changes – in the mouse ventral hippocampus – with repeated exposure. Using bulk RNA-seq across multiple time points, we observed widespread transcriptional habituation: stress-induced gene expression was blunted without the emergence of new response profiles or changes in baseline gene expression. Temporally resolved single-nucleus multi-omic profiling (RNA and chromatin accessibility) confirmed that the transcriptional habituation occurs across cell types, revealing cell type-specific patterns. We identify gene clusters whose expression habituates with repeated exposure, which appear primarily related to two distinct mechanisms: an early blunting of cAMP-associated genes chiefly linked to a reduced number of activated cells, and a shortened transcriptional response of corticosterone-associated genes, independent of the number of activated cells. Our study provides a comprehensive resource, accessible as an interactive app, that characterizes the dynamic transcriptomic response to stress in the mouse hippocampus and how it habituates with repeated exposure to stress.

Introduction
In moments of acute stress, our bodies rapidly mobilize energy resources to ensure survival and initiate a complex molecular response to meet increased metabolic demands1. These widespread molecular changes resemble a symphony orchestrated by numerous neurotransmitters, peptides and hormones1,2. Temporally, the response can be broadly divided into fast and slow components. The fast component, occurring within milliseconds, is mediated by the release of neurotransmitters and peptides like glutamate, GABA, corticotropin releasing hormone (CRH), noradrenaline or dopamine. The slow component is mediated primarily by the hypothalamus-pituitary-adrenal (HPA) axis, with corticosterone (CORT) reaching the brain only several minutes after stress initiation3–7. At the transcriptional level, where our molecular understanding of the acute stress response is most detailed, gene expression unfolds in waves and concludes within several hours after stress initiation8,9. By manipulating individual signaling pathways, studies begin to reveal which stress-mediators contribute to these stress-induced changes4,10–12.

To prevent wear-and-tear, exhaustion, and the onset of stress-related disorders, the stress response must adapt when stress occurs repeatedly and becomes chronic13,14. However, adaptation to stress can come in many flavors. First, repeated stress exposure can lead to habituation, the reduced responsiveness to a repeated stimulus15,16. Habituation has been described, particularly for the HPA axis, across a variety of different stressors17–20. On a molecular level, stress-induced expression of immediate early genes (IEGs, like Fos or Arc) also typically habituates after chronic stress exposure21–25. Second, stress sensitization can occur when repeated presentation of the same stressor, or different stressors, triggers a more intense response than the first presentation. This process is often associated with heightened arousal in post-traumatic stress disorder, or with increased cortisol release from hypertrophied adrenal glands26–30. Third, chronic stress exposure can lead to the emergence of novel molecular and circuit adaptations that differ qualitatively from acute stress exposures, effects that have been studied particularly in the context of stress resilience14.

On a genome-wide level, it remains unknown how the acute stress response itself adapts after repeated stress exposure, and which mechanisms might be driving it. Most of the work on the transcriptomic consequences of chronic stress has focused on long-term alterations that manifest in the brain after extended periods of chronic stress31. Large gene expression screens have revealed widespread changes in stress-sensitive regions like ventral hippocampus (vHC), prefrontal cortex, hypothalamus or nucleus accumbens, but tissue was collected long after the termination of the last stressor32–35. To understand how the transcriptomic response to the stressor itself adapts when the stressor is presented repeatedly, we chose a widely used restraint stress model and dynamically profiled multi-omic changes across multiple time points after the first and subsequent stress exposures in the mouse hippocampus. We describe, in both sexes, how the transcriptional response to acute restraint stress unfolds over time. We show profound transcriptional habituation and blunted chromatin remodeling across different cell-types as the stressor becomes chronic, and we identify two distinct mechanisms associated with different temporal components of stress habituation.

Results
Transcriptomic dynamics after acute restraint stress in the hippocampus
Acute restraint stress (ARS) is a commonly used stress model, which activates the HPA axis and lead to strong transcriptional responses in various brain regions25,36. Because the dynamic molecular changes triggered by ARS have never been characterized in detail, we exposed adult male and female mice to 1h30 of restraint stress, and collected the ventral hippocampus (vHC) for bulk RNA-sequencing at baseline, 45min, 1h30, 3h, 5h30 and 24h after stress initiation (Fig. 1A). Consistent with previously described transcriptional responses to various acute stressors9,37, we detected strong gene expression changes 45min after stress initiation, which subsequently evolve in waves and resolve within several hours (Fig. 1B). While some differentially expressed genes were still detected 5h30 after stress initiation, none were detected 24h afterwards (Fig. 1D). These dynamics and the overall return to baseline is also visible in a multi-dimensional scaling representation of the gene expression response (Fig. 1C, left). In previous work we showed that after an acute, short swim stress exposure, transcription had fully returned to baseline already 4h after stress exposure9. To test whether the effects at 5h30 represent lingering transcripts from earlier transcription, we looked at unspliced reads as a proxy for active transcription. We no longer observed significant group differences at 5h30, suggesting that transcription had indeed returned to baseline at this point (Fig. 1C, right).

Global transcriptomic habituation in response to repeated stress
To examine how the molecular response to ARS changes after repeated stress exposure, we subjected male and female mice to chronic restraint stress (CRS). CRS is one of the most widely used chronic stress models, it leads to habituation and results in long-term behavioral, molecular, and structural changes in the hippocampus and other brain regions38. As the transition from acute to chronic stress is hard to define, we exposed mice to 10 or 20 days of daily 1h30 restraint stress. The vHC was collected for bulk RNA-sequencing at baseline, 45min, 1h30, 3h, 5h30 and 24h after the last stress initiation (Fig. 1A). Due to large cohorts, the experiments had to be conducted separately. Thus, additional handling control groups were added to both the 10d and 20d CRS experiments, for which tissue was collected at baseline, 45min and 24h after an ARS exposure. This allowed us to compare groups across cohorts, but it also enabled detailed within-cohort comparisons. Each group contained 8 mice, and each mouse served as a biological replicate (192 mice in total, a detailed schematic of the experimental design is provided in Supplementary Fig. 1).

Across both sexes and treatment durations (10 and 20 days of CRS), we observed a profound habituation of the transcriptional response to ARS in repeatedly stressed animals. Virtually all stress-responsive genes showed a reduced response in at least one time-point in CRS animals (Fig. 1E). Surprisingly, we did not observe any overall difference in the magnitude of habituation between 10 and 20 days of CRS, as indicated by a strong correlation between the CRS:ARS interaction fold-changes, which captures how each CRS treatment impacts the response to ARS (Fig. 1F and Supplementary Fig. 2A). This reveals that, at the level of gene expression, 10 days of CRS are sufficient to elicit strong transcriptional habituation of the acute stress response. Figure 1G shows the change in expression of prominent stress-responsive genes across the 3 treatment groups. Finally, we could not observe any significant expression differences between CRS and handling animals at baseline (Fig. 1H). Although surprising given the extensive literature reporting long-term changes following chronic stress exposure38, this indicates that the habituation after CRS is not the result of a change in baseline gene expression. Thus, we next explored which mechanisms could account for the blunted gene expression response.

Different mechanisms are associated with the temporal dynamics of habituation
The profiling of transcriptional changes over multiple time points allowed us to follow the temporal dynamics of stress habituation in unprecedented detail. To quantify these dynamics, we used the angle on the centered multidimensional scaling plot as a pseudo-temporal labeling of the samples across the global stress response (Fig. 2A). This revealed that—in comparison to the first stress exposure—the global transcriptional response in CRS animals (10 or 20 days) was reduced at 45min. At 1h30 the molecular response peaked similarly in all groups, but returned to baseline much faster in CRS animals (Fig. 2B). To further explore these different temporal dynamics, we used spectral clustering to capture groups of genes with similar transcriptional dynamics across time (Supplementary Fig. 2B). Having previously shown that most transcriptional downregulation following acute stress is explained by RNA half-lives rather than direct transcriptional control9, we restricted the clustering to upregulated genes (see downregulated in Supplementary Fig. 2B and F). We resolved three main habituation patterns (Fig. 2C): Cluster 1 captured transcripts that were blunted already at 45min. Cluster 2 captured transcripts that peaked normally in CRS animals, but returned to baseline faster. Cluster 3 captured transcripts that were normally upregulated at 3h but failed to do so in CRS animals. In the unspliced RNA fraction, most genes in cluster 3 show comparable active transcription at 1h30 (Supplementary Fig. 2C), suggesting, also in this case, a faster return to baseline. Enrichment analyses comparing the genes in each cluster to all expressed genes or to all stress responsive genes revealed converging results (Fig. 2D and Supplementary Fig. 2E): The genes blunted already at 45min (cluster 1) were especially enriched for CRH and Orexin receptor signaling pathways, as well as GO terms related to the vasculature (Fig. 2D). Genes in this cluster particularly included neuronal IEGs, such as Fos, but also early response genes from the vasculature such as Apold1 (Fig. 2E, left). Instead, genes that peaked normally at 1h30 but returned to baseline faster in CRS animals (cluster 2) were enriched for response to abiotic stimulus and oxygen levels, as well as response to corticosteroids (Fig. 2D). In line with this, one of the most specific glucocorticoid-responsive genes, Fkbp539,40, was normally upregulated at 1h30, but significantly blunted at later time points (Fig. 1G), and the glucocorticoid regulated kinase Sgk1 showed a similar pattern (Fig. 2E, center). The genes ablated at 3h (cluster 3) did not show enrichment that passed multiple testing correction, yet some exemplary genes are shown (Fig. 2E, right panel). Testing whether the three clusters showed enrichment for curated gene targets of transcription factors (based on the collecTRI regulons41) yielded corroborating results (Supplementary Fig. 2D): Genes in cluster 1 were most prominently enriched for targets of the cAMP response element binding CREB1 (2.3-fold, FDR~5e-4), as well as other factors including the Serum Response Factor (SRF). In contrast, genes from cluster 2 were most strongly enriched for targets of the glucocorticoid receptor (GR) (3-fold, FDR~1e-7), with no target enrichment for cluster 3. To further confirm the involvement of GR signaling, we compared the response to ARS with the hippocampal gene expression response to an IP injection of the GR agonist dexamethasone (DEX)11. The vast majority of genes upregulated 3h after ARS (which are all blunted in CRS animals, whether in cluster 2 or 3) were also upregulated by DEX injection (Fig. 2F and G). Conversely, we had previously shown that in vitro treatment of CA1/CA3 hippocampal slices with the cAMP activator forskolin42 largely recapitulates—in the absence of corticosterone signaling—many of the early gene expression changes seen in response to acute stress exposure1. Indeed, we replicate this in the current dataset, where most of the genes blunted at 45min (cluster 1) were upregulated by forskolin in vitro (Fig. 2H). In contrast, this was true only for a small minority of those genes that return faster to baseline (cluster 2), and only for 1 gene from those genes with ablated expression at 3h (cluster 3). These observations suggest that there are two mechanisms, a rapid, GR-independent early blunting that is related to cAMP signaling, and a slower, GR-dependent blunting that emerges at later time points.

Behavioral and neuroendocrine outcomes of habituation
We were surprised that 10 days of restraint stress led to strong molecular habituation that was indistinguishable on the whole transcriptome level from 20 days of restraint. Thus, we asked whether behavioral changes, typically reported after 3 weeks of restraint stress exposure34, are also apparent after 10 days. Therefore, we compared the behavioral profiles of control, ARS and CRS male and female mice in the open-field test 45min after the end of their first ARS exposure, or after 10 or 20 days of CRS exposure (Fig. 3A and B). In general, we observe an increase in distance moved and supported rearing upon both 10d and 20d CRS, indicating an increase in general activity compared to control animals (Supplementary Fig. 3A and C), with male and female animals generally responding similarly (Supplementary Fig. 3B and D). Additionally, we find a trend in CRS animals spending less time in the center of the arena. This suggests that behavioral changes emerge 10 days after stress exposure, which are similar to the changes observed after 20 days of stress exposure. Given the comparable behavioral response, and to increase power, we combined the animals of both experiments. This global analysis shows strong differences between those animals being stressed for the first time and those being repeatedly stressed. After CRS, mice are hyperactive (more movement and more supported rearing) and spend less time in the center, suggesting increased anxiety (Fig. 3C and D).

Given that habituation of the HPA axis has been well-documented, and because we see a clear involvement of GR in the blunting of certain clusters of genes (Fig. 2C and F), we measured corticosterone (CORT) levels in the blood in response to an ARS challenge presented either for the first time (ARS) or after 10 days of restraint stress (CRS-ARS) (Fig. 3E). In both sexes, CORT levels rose after restraint in both conditions, with only a modest blunting of the response in males and no detectable blunting in females (Supplementary Fig. 4A). To verify these findings, we repeated the experiment in an independent replication cohort in which we also collected the hippocampus after the final blood sampling to assess the impact of CRS on CORT and transcriptional blunting in the same animals (Supplementary Fig. 4B). This replication confirmed that after 10 days of CRS both males and females still mounted a strong CORT response (Supplementary Fig. 4B). When data from both cohorts were pooled, a moderate blunting of CORT levels was detectable only in males at the 45min time point, but CORT levels peaked similarly at 90min (Fig. 3F). RNA-seq confirmed that hippocampal gene expression changes were again blunted after CRS in both sexes (Fig. 3G), showing a disconnect between CORT levels and transcriptional blunting. In line with this, we found no correlation between CORT levels and expression of the GR-sensitive gene Fkbp5 (Supplementary Fig. 4C). To compare CORT levels with the overall strength of transcriptional habituation within each animal, we developed a response coefficient score by regressing each sample’s response on the average logFCs of the differentially-expressed genes at 3h in the original ARS cohort (Fig. 3H and methods). This confirmed the strong habituation effect at the transcriptional level (Fig. 3I), but revealed no significant correlation between CORT levels and the transcriptional response coefficient (Supplementary Fig. 4D). Together, these results strengthen the conclusion that after 10 days of CRS, peripheral CORT levels peak normally in males and females, while the transcriptional response in these animals is strongly blunted.

Whole-tissue chromatin accessibility does not explain molecular habituation
To test whether changes in chromatin prevent stress-induced activation of gene expression after repeated stress exposure, we performed ATAC-seq of male and female mice that were either handled or stressed for 10 days both at baseline and 45min after ARS initiation (Fig. 4A). Notably, analyses were performed in the same animals that were used for RNA sequencing, by using the contralateral hemisphere of the ventral hippocampus. We chose the 45min time point, because earlier work had shown that within an hour after stress exposure clear chromatin remodeling can be observed43. Indeed, we observe differences (particularly increases) in peak accessibility, upon ARS (Fig. 4B and C). A differential promoter accessibility analysis yielded several genes with increased accessibility upon ARS (Fig. 4D), including known stress-responsive genes such as Fkbp5 or Nfkbia, and lesser-known genes like Phyhd1 (Fig. 4D and E). We correlated changes in promoter accessibility upon ARS with changes in RNA, using the union of genes that were significant in either ATAC-seq or RNA-seq. Since RNA changes persist much longer than changes in active transcription, and DNA accessibility relates to the latter, we focused on unspliced RNA at 45min, and observed a weak, but significant correlation (Fig. 4F, left). Of note, some genes (e.g., Hif3a and Fkbp5) show an increase in promoter accessibility already at 45min, while the effect on gene expression can only be observed at later time points, indicating that chromatin accessibility foreshadows changes in transcription. We therefore performed the same comparison with the transcriptome at the next time point (1h30), and observed a slightly higher correlation (Fig. 4F, right).

We next searched for transcription factors whose putative binding sites changed in accessibility (Fig. 4G and H). The most significant motif was of the glucocorticoid receptor (as well as the very similar motifs of other hormone receptors), which showed an increase in accessibility, followed by MEF2-related motifs, showing decreased accessibility. Finally, we also observed a weaker, but significant increase in accessibility of the JUND motif (Fig. 4G and H). While we could identify major changes in accessibility upon ARS, we could not detect any change associated with CRS, either at baseline or in interaction with ARS. This means that 10 days of CRS do not alter ARS-induced chromatin remodeling. Even when focusing on the genes blunted at 45min or 1h30, we could not observe any general pattern in the accessibility of their promoter in response to CRS. Thus, chromatin accessibility does not explain the habituation in gene expression after repeated stress exposure. However, accessibility changes occur faster and are more transient than RNA changes, thus it is possible that we have missed CRS-associated changes in accessibility and TF activity that have resolved at the 45min time point. Additionally, it is possible that bulk ATAC-seq lacks the sensitivity to detect CRS-related changes, especially if they happen in a subset of cells. To address these two caveats, we turned to single-nuclei analysis of both transcriptome and chromatin accessibility across several time points.

Transcriptomic single-cell dynamics of the acute stress response
It has recently been shown that stress-induced changes have cell-type specific effects, and that the analysis of bulk tissue can obscure biologically meaningful changes44,45. To test how habituation to the acute stress response evolves over time across cell types, we exposed female mice to ARS after 10 days of handling (handling-ARS), or after 10 days of chronic restraint stress (CRS-ARS), and collected the vHC at 15min, 45min and 3h after initiation of the last stress exposure (n=3 mice per group) (Fig. 5A). We performed simultaneous single-nucleus RNA-seq and ATAC-seq on the same nuclei (10X Multiome). We sequenced a total of 215’349 nuclei across 23 samples, and after stringent quality filtering (see methods) were retained ~160’000 cells with a median of 7343 UMIs per cell. Non-neuronal cells included oligodendrocytes and oligodendrocyte progenitor cells (OPC), with a fraction of committed oligodendrocyte progenitors in-between (Fig. 5B). Microglia and astrocytes had prominent clusters, and we could also identify perivascular macrophages (PVM). Although lower in number, we could observe clear populations of pericytes, endothelial cells, smooth muscle cells, and some smaller sub-populations of vascular cells with unclear signatures, and we merged them as ‘vascular’ cells for downstream analysis.

The different neuronal subpopulations are shown in Fig. 5C (marker expression is shown in Supplementary Fig. 5A). Inhibitory neurons clustered by known type (Sst, Pvalb, Vip, Ntng1 interneurons, and Lamp5 neurogliaform interneurons), while excitatory neurons clustered by region, first with a clear distinction between granule cells from the dentate gyrus, mossy cells from the hilus, and CA1 and CA3 excitatory neurons. Within the pyramidal layer (i.e., CA1-3), cells then clustered by relatively more dorsal vs ventral regions (Fig. 5C and D) based on signatures from Cembrowski et al.46. Dorsal scores were generally lower, and fewer cells had a high dorsal score, in line with the fact that only ventral hippocampi were analyzed (Supplementary Fig. 5B). This shows that the ventral-dorsal gene expression gradient is preserved–and can be resolved–even within the ventral portion. In addition, we identified a gradient of cells spanning from the CA1 to the pro-subiculum and the distal subiculum (Fn1+), as well as a separate population most likely from the posterior subiculum (with Ly6g6e and Cntnap4 expression). Another gradient especially characterized by high Camk2d expression extended from the ventral CA1 population. Based on correlations with the Allen brain atlas, they appear to be neurons from the hippocampus-tegmental amygdala (HATA) and/or hippocampus-inferior temporal cortex interfaces. These populations were the only ones to (mildly) vary in abundance across samples, consistent with their location on the border of the dissected tissue. Of note, the transcriptome-based cell type labels were also clearly segregated in the ATAC-seq embedding, both for non-neuronal (Fig. 5E) and neuronal (Fig. 5F) populations, corroborating the robustness of the cell type assignment.

To date, no dynamic profiling of the acute stress response on the single-cell level exists in the hippocampus. Therefore, we first examined the transcriptional changes to ARS on a pseudo-bulk level (Supplementary Fig. 6). Globally, the different broad cell types all responded to ARS early on (Fig. 5G). In neuronal cell types, the response peaked at 15min (inhibitory neurons) or 45min (excitatory neurons), before strongly decreasing at 3h. In contrast, the response of all non-neuronal cell types was sustained at 3h. In neurons, a particularly strong response was observed in CA1 pro-subiculum and VIP interneurons, which upregulated a number of major IEGs already 15min after stress initiation, such as Fos/Fosb, Npas4, and Btg2 (Fig. 5H). With the exception of vascular cell types, which also responded very quickly, other non-neuronal populations were slower to respond (Fig. 5I). The significant genes in non-neuronal populations included a number of prominent stress-induced genes, such as Apold1 in vascular cells47, the cell cycle regulator Cdnkn1a, or the GR target Fkbp5 across multiple cell types. Of note, some changes, such as the strong upregulation of the cAMP/cGMP-degrading enzyme Pde10a specifically in astrocytes upon stress (Fig. 5I and Supplementary Fig. 6), had not been detected in the bulk transcriptome, presumably because the gene is also highly expressed in other glial populations at baseline, underlining the value of single-cell profiling.

Transcriptomic habituation to repeated stress occurs across cell-types
Visual inspection of the heatmaps in Fig. 5H and I suggest that across cell types, the response to ARS was generally blunted after CRS, with different temporal patterns emerging for different cell types. To investigate this systematically, we tested, within each cell type, which stress-induced gene expression profiles were altered by CRS at any time point. Several genes were significantly affected by CRS, with a clear pattern of ARS-responsive genes habituating in CRS animals (Fig. 6A). This was also the case for several prominent stress-responsive genes (Fig. 6B). Notably, while Npas4 upregulation appeared to be completely ablated in CRS animals at 45min based on the bulk transcriptomic data, the single-cell data revealed an upregulation at 15min also in CRS animals across multiple excitatory neuron subtypes (Supplementary Fig. 7B). However, in CRS animals this upregulation was weaker, and had already returned to baseline at 45min across cell types, while it remained upregulated in animals who experienced the stressor for the first time (Fig. 6B and Supplementary Fig. 7B), consistent with the bulk data. In contrast to our bulk sequencing data, a number of genes appeared differentially expressed at baseline in CRS animals (Fig. 6C). While some of these candidates showed a variable expression pattern across ARS time-points (e.g., Cacna1e and Vwf), some looked more coherent, such as Pde1a in VIP interneurons, a gene that is involved in cGMP/cAMP homeostasis48. To ascertain these findings, we also specifically tested for differences at baseline, which revealed some statistically significant candidates (Supplementary Fig. 7A). While the identification of rare, cell-type specific baseline gene expression changes emphasize the power of single-cell sequencing, these limited changes cannot explain the widespread habituation of stress-induced transcriptional activity across cell types.

After habituation, fewer cells respond to stress
A blunted transcriptomic response to stress could be due to fewer cells responding, or to the same proportion of cells having a weaker response. cFOS or ARC stainings for instance indicate that only a fraction of hippocampal neurons respond to salient, stressful events36. Leveraging the power of single-cell transcriptomics, we estimated the recent activity of single neurons using an activity-dependent transcription (ADT) score based on genes that are differentially expressed upon ARS at the pseudo-bulk level (see methods). The ADT score indicates, for each cell, whether it looks more like a cell that does not respond to ARS (score of 0) or like the most stress-responsive cells of that type (score of 1). As expected, while only a minority of the cells showed ADT, these were unequally spread across subpopulations (Fig. 7A), with the highest activity in subiculum, pro-subiculum and CA1 neurons. Using an ADT score above 0.5 as a threshold for an activated neuron, we observed that the proportion of activated neurons doubled upon the first ARS exposure, a response that was blunted after CRS (Fig. 7B). In contrast, while the median ADT score of activated neurons did show a significant increase upon ARS, the difference was very small in magnitude (roughly 5%, Fig. 7C). Together, this indicates that most of the blunted response is due to fewer neurons being activated, rather than the same number of neurons being activated less strongly.

We next repeated this analysis at the level of individual cell types, which showed that the habituation after CRS is cell-type specific, with the strongest blunting of ADT observed in VIP interneurons and in several excitatory cell populations (Supplementary Fig. 8). To corroborate these results, we also used an independent method to estimate neuronal activation from single-cell data, neuroestimator49. Although the estimated baseline activation levels across subpopulations varied depending on the methods, the main patterns upon ARS and CRS were largely reproduced for inhibitory cell types, with less congruence in excitatory populations (Supplementary Fig. 8). Figure 7D shows the proportions of each cell type activated based on the average of the single-cell scores from the two methods. Consistently, VIP interneurons showed the highest change in stress-induced activity. Strong responses in ADT were also observed in CA1 and subiculum populations, in mossy cells, as well as in other inhibitory neuron subtypes (Pvalb and Ntng1). Notably, the temporal dynamics varied across cell types, e.g., VIP interneurons were strongly activated already at 15min, while other inhibitory neurons showed a more delayed activation. In all these populations, the strong ARS-induced change was blunted or completely absent after CRS, especially at 45min (Fig. 7D).

Contrary to the neuroestimator method, our ADT approach is applicable to any cell type, and we thus used it to estimate the proportion of cells recruited across non-neuronal populations (Fig. 7E and F). Only a minority of vascular cells reacted to ARS, whereas the majority of astrocytes and oligodendrocytes strongly reacted upon first ARS exposure, with microglia and OPC showing intermediate activation. In all glial cell types, but not in vascular cells, the ARS-induced activation was blunted in CRS animals, especially at 3h. This reduction was particularly strong for astrocytes, dropping roughly by half already from 45min on. In contrast, oligodendrocytes seem to be recruited similarly in response to stress, but while their activity is sustained across 3 hours in handling animals, activity drops sharply at this time point in CRS mice.

Repeated stress shortens GRE accessibility across cell types
To investigate the mechanisms underlying the observed transcriptional changes following acute stress, and the blunting after chronic stress exposure, we next turned to the changes in chromatin accessibility across cell types. The strongest and most significant ARS-induced changes in motif accessibility consistently came from the glucocorticoid receptor and highly similar steroid hormone receptor motifs (Supplementary Fig. 9), with increased accessibility already at 15min and extending to later time points (Fig. 8A). While the ARS-induced change in GR motif accessibility was similar between handling-ARS and CRS-ARS at 15min across all cell types, we observed clear transcriptional habituation in the CRS-ARS group at 45min in many cell types (Fig. 8A). After 3h, GR motif accessibility completely returned to baseline in the CRS-ARS group, but not in handling-ARS animals (Fig. 8A). Notably, the increase in GR motif accessibility upon ARS was much weaker in interneurons, which also did not upregulate its target Fkbp5. Expression of the GR gene (Nr3c1) is low in interneurons, and across cell types we observed that the magnitude of the increase in GR motif accessibility in response to ARS was proportional to the expression of Nr3c1 (Fig. 8B and Supplementary Fig. 10A).

While the accessibility of glucocorticoid response elements (GRE) already increases 15min after stress exposure, GR-responsive genes – like Fkbp5 – show transcriptional changes only at later timepoints. Across many cell types, Fkbp5 expression increased at 45min, and expression was significantly blunted 3 hrs after stress initiation (Supplementary Fig. 10B and Fig. 8B). These results suggest that, in both handling and CRS animals, glucocorticoids reach the brain and activate GREs already within 15min after initiation of stress, but that following CRS, GRE accessibility is shortened, and accordingly GR-dependent transcription shuts down more rapidly. Further, it supports the hypothesis that much of the transcriptional habituation observed at later time points is GR-dependent, while the blunted response observed at 15min is not.

Repeated stress induces differences in cAMP-associated TFs independent of GR activation
Apart from GR, we identified several other TF motifs that strongly respond to ARS (Supplementary Fig. 9; Fig. 8C–E), as well as candidate motifs showing baseline differences in CRS animals (Fig. 8F). In non-neuronal cell types, we particularly observed strong activation of CEBPA in microglia and NFIA in oligodendrocytes (Fig. 8D). In neurons, across stress-responsive neuronal cell types (i.e., especially CA1, pro-subiculum and VIP neurons), we could observe a strong increase in the accessibility of FOS/JUND motifs at all time points following ARS (Fig. 8C and E). However, this was not well coupled to changes in the accessibility of cAMP-response elements (CRE): while the CREB1 motif appeared to increase in accessibility at 15min in Ntng1 and VIP interneurons (Fig. 8E), this was not the case in most excitatory neurons, where it instead showed a trend in the other direction (Fig. 8C). This was especially surprising given the clear recruitment of these neurons and the known rapid CREB-dependent transcription upon depolarization50,51. To confirm the expected involvement of CREB in this paradigm, we tested motif enrichment within regions accessible in neurons in proximity to genes that respond to stress rapidly (i.e., genes that are significantly increased by stress at 15min in single-nuclei data or at 45min in the spliced bulk transcriptome). Amongst several over-represented motifs, CREB1 showed the strongest enrichment (Supplementary Fig. 10C). A similar analysis using experimental binding data highlighted a number of factors (Supplementary Fig. 10D), but CREB1 was the only significant factor in the combined analysis (Supplementary Fig. 10E).

To further investigate the seemingly inconsistent CRE accessibility pattern, we looked at motif accessibility across activated (defined by ADT>=0.5) and non-activated excitatory neurons (Fig. 8G–J). In activated neurons, we observed a strong increase in motif accessibility across several fast responding TFs like CREB1, SRF and JUN/FOS family members (Fig. 8G–I), confirming the expected increase in CRE accessibility in activated compared to non-activated neurons. The mild decrease in CRE accessibility in non-activated neurons (which are the majority) likely masks, in pseudo-bulk of the entire population (Fig. 8C), the increase from activated neurons. The accessibility of the TF’s respective motifs had differing temporal dynamics, which likely relate to how sustained the activity of the different TFs is following activation: SRF strongly spiked at 15min but returned close to baseline already at 45min (Fig. 8H), while JUND appeared more comparably active across activated neurons from different time points (Fig. 8I). Of note, the increase in SRF motif accessibility in activated neurons was significantly blunted in CRS animals, suggesting a difference in activation that cannot be attributed to fewer neurons being activated. This could be due to a difference in SRF response to neuronal activity or, given its fast cAMP-independent. Our single-nucleus ATAC-seq analysis further confirms the involvement of CREB- and SRF-dependent signaling.

What remains unclear, however, is whether after habituation some IEGs are blunted from the start, or whether their initial response is preserved but transcription shuts down faster. Npas4, for example, shows complete habituation 45min after repeated stress in our bulk RNA-Seq, but we still find a significant – albeit weak – increase in Npas4 transcripts at 15min in several cell types in our single-nucleus data. In a few cell types, however, Npas4 seem to be completely silenced across time points (Supplementary Fig. 7B). Due to the lack of a finer temporal resolution within the first 15min of stress initiation, we cannot fully resolve the question whether IEG transcription is similarly initiated between ARS and CRS. However, our method for quantifying activity dependent transcription (ADT) allows us to make two critical observations: First, we find that, across cell types, the blunted transcriptional response to ARS after a history of CRS is chiefly due to fewer cells being activated, rather than the same proportion of cells showing a weaker activation. Second, activated neurons show a weaker increase in SRF motif accessibility at 15min (Fig. 8), which could indicate either a shorter or weaker activation response.

GR-dependent transcriptional habituation
Acute restraint stress strongly induces HPA axis activation. While CORT levels in the blood rise within a few minutes, CORT is detected in the brain after 10min5,6. On the bulk and single-cell level, our transcriptomic analyses reveal genes that peak similarly at 1h30 after stress onset, but shutdown faster and return to baseline levels at the 3h time point in animals with a history of chronic stress, while at that time point still being strongly upregulated after the first ARS exposure. Converging evidence shows that these transcripts are enriched for known GR target genes, such as Fkbp539,40, and they overlap with DEX-responsive genes from an independent analysis. In both our bulk and single-cell ATAC sequencing, the strongest ARS-induced changes were associated with an increase of GR motif accessibility, consistent with previous research43. After CRS exposure, this ARS-induced increase is largely blunted at 45min and 3h across cell types. However, we do observe a strikingly similar increase in GR motif accessibility at 15min (Fig. 8A), thus indicating that indeed GR activation is initiated similarly to ARS, but shortened. Importantly, across neurons this late component of transcriptional habituation is independent of whether the cell was activated or not. CORT levels measured in the blood support the notion that the initial rise is similar in animals with or without a history of prior stress. However, at later time points, CORT dynamics paint a more complex picture, suggesting that CORT levels habituate slightly in males, but not in females. Notably, several reports find that the CORT levels habituate strongly in response to repeated restraint stress18–20, yet not all studies observe this habituation57,58 and we suspect that varying the time of day (unpredictable restraint stress) might have prevented HPA-axis habituation in our setup. While our results support previous reports of sex differences in HPA-axis regulation and a more pronounced CORT response in females59–61, they contrast with the fact that we observe no corresponding sex differences in overall transcriptomic habituation patterns. These data suggest that in males, a blunting of the HPA-axis response could account – at least to a small degree – for the blunted expression of GR-sensitive genes after repeated stress exposure. However, we think this link is weak for several reasons: (1) The CORT-blunting observed in males at a single timepoint is a small effect, whereas the molecular habituation is strong and occurs across time points. (2) We observe no clear correlation between CORT levels and gene expression changes (Supplementary Fig. 4). (3) In females, where the transcriptomic blunting looks very similar to males, we detect no changes in CORT levels after repeated stress exposure. Notably, sex differences in glucocorticoid binding globulin have been reported60, thus it remains to be explored if this could contribute to differences in stress habituation. Another possibility, supported by the observation that the initial rise of GRE accessibility is comparable between groups, but returns to baseline faster in CRS animals, could be that the GR-dependent blunting is due to differences in the retention of GR in the nucleus, for instance through nuclear export or what has been proposed as the chaperone cycle62–64. More work is needed to reveal the mechanisms behind these unexpected effects.

Methods
Animals
C57BL/6J (C57BL/6JRj) mice were obtained from Janvier (France). Mice were maintained in a dedicated animal facility on a 12h reversed light–dark cycle (lights off at 09:30 a.m.) with food and water ad libitum. The temperature in the animal rooms was maintained at 22–24°C with a maximum deviation of ±1°C. The relative humidity was kept at 50%, with a maximum deviation of ±10%. Male and female mice were housed, separated by sex, in groups of 2–5 animals per cage and used for experiments when 2–3 months old. All tests were conducted during the animals’ active (dark) phase. All experimental procedures were carried out in accordance with the guidelines of the Federal Veterinary Office of Switzerland and were approved by the Cantonal Veterinary Office in Zurich.

Stress paradigm
Mice were exposed to restraint stress for 1h30 in a 50ml Falcon tube with a large air hole. All stress exposures were carried out between 10:00 a.m. and 5:00 p.m. Animals were single-housed 24h prior to the final stress exposure to reduce confounding effects induced by cage mates on the day of the experiment61. For acute stress experiments, animals were exposed to one restraint session (1h30). The restraint tube was placed in the home cage for the duration of stress. For chronic stress experiments, animals were subjected to either 10 or 20 consecutive days of daily restraint stress at varying times of the day. Control animals and animals exposed to acute stress were gently handled daily for an equal number of days. Gentle handling involved lifting the animals by the tail, simulating the handling method used before restraint stress.

CORT measurement
To assess CORT levels, blood samples were collected via tail prick at 0min, 15min, 45min, 90min and 180min after the onset of restraint. Blood samples were stored at RT for 30 min or 4°C overnight to allow for clotting and centrifuged at 2500g for 10min at 4°C. The serum supernatant was collected and stored at −80°C until further analysis. Corticosterone levels in the serum were measured using an ELISA kit (Invitrogen, EIACORT) according to the manufacturer’s instructions. A 1.5μl serum sample was diluted 1:100 for the assay. Absorbance was measured at 450nm with a reference at 560nm using a plate reader, and final concentrations were determined based on the provided standards.

Tissue collection
At the indicated time point after initiation of stress, mice were euthanized by cervical dislocation and decapitation. When tissue was collected at 15min, 45min or 1h30 after stress initiation, mice were removed from the restraint tube at the indicated time point and immediately euthanized for brain collection. The brain was quickly dissected on a cool surface and isolated ventral hippocampi were snap-frozen in liquid nitrogen and stored at −80°C until further processing.

Whole tissue RNA extraction
One snap-frozen ventral hippocampus per animal was used to obtain transcriptional profiles. Samples were processed in a multiple block design with every block containing one replicate of each condition and sex. Processing within blocks was randomized. Samples were homogenized in 500µl Trizol (Invitrogen 15596018) using metal beads (Qiagen 69989) in a tissue-lyser bead mill (Qiagen, Germany) at 4°C for 2min at 20Hz. RNA was extracted according to the manufacturer’s recommendations. RNA purity and quantity were determined with a UV/V spectrophotometer. RNA integrity was assessed with high sensitivity RNA screen tape on an Agilent Tapestation/Bioanalyzer, according to the manufacturer’s protocol.

Whole-tissue library preparation and sequencing
For experiment 1 (see Supplementary Fig. 1), library preparation and sequencing was performed at the Functional Genomics Center Zurich (FGCZ) of University of Zurich and ETH Zurich. The TruSeq stranded RNA kit (Illumina Inc.) was used according to the manufacturer’s protocol. For bulk sequencing library preparation, the TruSeq stranded RNA kit (Illumina Inc.) was used according to the manufacturer’s protocol. The mRNA was purified by polyA selection, chemically fragmented and transcribed into cDNA before adapter ligation. Single-end (100nt) sequencing was performed with HiSeq 4000. Samples within experiments were each run on one or multiple lanes and demultiplexed. A sequencing depth of ~20M reads per sample was obtained.

For experiments 2+3 (see Supplementary Fig. 1), library preparation and sequencing was performed at Novogene, UK. Briefly, RNA samples were used for library preparation using NEBNext® Ultra RNA Library Prep Kit for Illumina®. Indices were included to multiplex multiple samples. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. After fragmentation, the first strand cDNA was synthesized using random hexamer primers followed by the second strand cDNA synthesis. The library was ready after end repair, A-tailing, adapter ligation, and size selection. After amplification and purification, insert size of the library was validated on an Agilent 2100 and quantified using quantitative PCR (Q-PCR). Libraries were then sequenced on Illumina NovaSeq 6000 S4 flowcell with PE150 according to results from library quality control and expected data volume. Samples within experiments were each run on one or multiple lanes and demultiplexed. A sequencing depth of ~20M paired reads per sample was obtained.

Bulk transcriptome analyses
Bulk RNA-seq was analyzed as described in ref. 9. Briefly, trimmed reads were used for quantification of the GENCODE M17 transcripts using kallisto 0.44, reads from protein-coding transcripts were summed to the gene level, and differential gene expression analysis was performed using edgeR 4.6.3. In addition to a FDR threshold of 0.05, an absolute logFC threshold of at least 0.15 was used to determine differentially-expressed genes (DEGs). Wherever surrogate variable (SV) analysis (using sva package version 3.50.067) identified more than 5 SVs, we restricted it to 5. To identify habituation patterns, sva-corrected, unit-variance-scaled logFCs of DEGs across samples from time points 45min to 3h were first extracted, and values ceiled between −5 to 5 to restrict outlier influences. Genes were then clustered using the Spectrum 1.1 package with k=10.

The MDS was computed through the plotMDS function of the limma package version 3.58.1 using the sva-corrected expression of the acute stress DEGs across Handling samples, and the CRS samples were then projected onto this space. The MDS was recentered so that all timepoint medians were equidistant to the origin before calculating angles.

Cluster enrichment analysis was done using Fisher’s exact test (followed by FDR correction) on the KEGG, WIKIPATHWAYS and GO gene sets of effective sizes ranging from 5 to 1000 from the MSigDB68, accessed from the msigdbr R package version 7.5.1. For the enrichment of motifs/ChIP peaks in early stress-responsive genes (Supplementary Fig. 9C–E), we first used liftover to translate the human brain TF peaks to mouse coordinates, and then overlapped the lifted neuronal peaks on the single-cell ATAC peaks. We then used Fisher’s exact test to test for enrichment of motifs or peaks in the broad promoter region (±10kb around TSS) of stress-responsive genes, compared to expressed genes.

To make the data from the Cort profiling cohort (Fig. 3J) comparable to the other batches, we used the empirical Bayes batch effect correction (ComBat, as implemented in the sva package) on log(CPM) data, treating the Handling + 3h ARS samples of each batch as coming from the same group. We then computed response coefficients by regressing each sample’s logFC of the genes significantly responding at 3h onto the corresponding average logFCs in the original cohort using a robust linear model without intercept, and used the slope/coefficient to assess overall response.

Nuclei isolation for ATAC-Seq
Nuclei isolation was performed as previously described69. All steps were performed on ice. Briefly, one vHC was homogenized in lysis buffer (26mM Sucrose, 30mM KCl, 10mM MgCl2, 20mM Tricine-KOH (pH 7.8), 1mM DTT, 0.5mM Spermidine, 0.15mM Spermine, 0.3% NP-40, Complete protease inhibitor in nuclease-free water) by 10 strokes with pestle A and 20 strokes with pestle B. The sample was filtered through a 70μm Flowmi strainer and centrifuged for 5min at 400g at 4°C. Supernatant was removed and additionally sis buffer was added to resuspend the pellet. The sample was then mixed with 50% Iodixanol solution (30mM KCl, 26mM Sucrose, 10mM MgCl2, 20mM Tricine-KOH (pH 7.8), 50% Iodixanol solution in nuclease-free water) and underlaid by 30% Iodixanol solution (50% Iodixanol solution in lysis buffer) and 40% Iodixanol (50% Iodixanol solution in lysis buffer) layers. The sample was centrifuged for 20min at 3000g at 4°C with brakes off. The nuclei band was collected and mixed with ATAC buffer (10mM NaCl, 10mM Tris-HCl (pH 7.5), 3mM MgCl2, 0.1% Tween 20 in nuclease-free water). The nuclei count was assessed and 50’000 nuclei were transferred to a new tube with ATAC buffer. The samples was centrifuged for 10min at 500g at 4°C and supernatant was removed.

Preparation of ATAC libraries
Nuclei were resuspended in 50μl of transposition mix (25μl 2× TD buffer, 2.5μl transposase (Diagenode, C01070012), 16.5μl PBS, 0.5μl 1% digitonin, 0.5μl 10% Tween-20, and 5μl nuclease-free water). Transposition reactions were incubated at 37°C for 30min while shaking at 1000 r.p.m. Tagmented DNA was purified using the MinElute Reaction Cleanup Kit (Qiagen, 28204). Purified DNA was amplified using the NEBNext high-Fidelity 2X master mix (New England Biolabs, Ipswich, Massachusetts, United States, M0541) and Diagenode primers (Diagenode, Marlborough, Massachusetts, United States, C01011036) at 72°C for 5min, 98°C for 30s and 5 cycles of 98°C for 10s, 63°C for 30s, and 72°C for 1min. A fraction of the amplification product was used to determine the number of additional cycles using qPCR, resulting in a total of 8–9 amplification cycles for all samples. A double-sided size selection (0.5X and 1.3X) was performed to remove large and very small library fragments. Insert size of the library was validated and quantified on an Agilent 2100 Bioanalyzer. Libraries were then sequenced on Illumina NovaSeq 6000 S4 flowcell with PE150 according to results from library quality control and expected data volume aiming for 50 Million read pairs per sample. Sequencing was performed at Novogene, UK.

Analysis of bulk ATAC sequencing
Reads were trimmed with trimmomatic and aligned to the GRCm38 genome using bowtie2 with --dovetail --no-mixed --no-discordant -I 15 -X 2000. Duplicates were marked using Picard. Peaks were called using Genrich 0.6.1 (https://github.com/jsh58/Genrich) with ATAC mode with a FDR threshold of 0.05, discarding PCR duplicates and excluding the mitochondria and ENCODE blacklisted regions. Quantification was performed by counting the number of insertion events (i.e., shifted fragment starts/ends) within peaks or promoter regions. Two samples (one from each of the two 45min groups) were excluded because they had much lower coverage, lower sequencing quality, longer fragments, a different GC distribution from the others samples, and were very clear outliers on an MDS plot. Differential analysis was performed with edgeR glmFit using a model including Sex, Treatment, ARS timepoint, and a Treatment:TimePoint interaction. For promoter-based analysis, an additional surrogate variable based on the SVA package was included. Insertion coverage tracks were generated (and plotted) with the epiwraps package v0.99.102, using binWidth=10L, shift=c(4L,-5L), type=“ends”, extend=3L, and includeDuplicates=FALSE.

Motif accessibility analysis was performed using chromVAR. For this purpose, peaks were resized to 300bp and insertions were counted as described above. 2000 iterations were used, and the z-scores were quantile-normalized as recommended in ref. 70. Differential analysis was performed with limma using 1 surrogate variable.

Nuclei isolation for single-cell experiments
Nuclei isolation was performed similarly as for ATAC-Seq, with the addition of RNase inhibitors as previously described69. All steps were performed on ice. Briefly, one vHC was homogenized in lysis buffer (26mM Sucrose, 30mM KCl, 10mM MgCl2, 20mM Tricine-KOH (pH 7.8), 1mM DTT, 0.5mM Spermidine, 0.15mM Spermine, 0.3% NP-40, Complete protease inhibitor, 0.2 U/μl RNase inhibitor in nuclease-free water) by 10 strokes with pestle A and 20 strokes with pestle B. The sample was filtered through a 70μm Flowmi strainer and centrifuged for 5min at 400g at 4°C. Supernatant was removed and additionally sis buffer was added to resuspend the pellet. The sample was then mixed with 50% Iodixanol solution (30mM KCl, 26mM Sucrose, 10mM MgCl2, 20mM Tricine-KOH (pH 7.8), 50% Iodixanol solution in nuclease-free water) and underlaid by 30% Iodixanol solution (50% Iodixanol solution in lysis buffer) and 40% Iodixanol (50% Iodixanol solution in lysis buffer) layers. The sample was centrifuged for 20min at 3000g at 4°C with brakes off. The nuclei band was collected and mixed with ATAC buffer (10mM NaCl, 10mM Tris-HCl (pH 7.5), 3mM MgCl2, 0.1% Tween 20, 0.2 U/μl RNase inhibitor in nuclease-free water) before the nuclei count was assessed. The sample was centrifuged for 5min at 500g at 4°C and supernatant was removed. Diluted nuclei buffer (0.1M DTT, 1 U/μl RNase inhibitor, 1X Nuclei buffer (10X) in nuclease-free water) was added to achieve the desired nuclei concentration.

Preparation of single-nucleus Multiome libraries
10X Multiome libraries were prepared following the Chromium Next GEM Single Cell Multiome ATAC + Gene Expression Reagent Kits User Guide (CG000338, Rev. F, 10X Genomics) with a targeted nuclei recovery of 10,000 nuclei (1 sample failed due to clogging). 7 cycles were performed for sample index PCR during ATAC library preparation. 6 cycles were used for cDNA amplification and a total of 13 cycles was performed for sample index PCR for gene expression library preparation. RNA integrity was assessed on an Agilent Tapestation, according to the manufacturer’s protocol. Samples of each library type were pooled. Both libraries (gene expressions and ATAC) were then sequenced on one Illumina X Plus 25B flowcell on individual lanes with an adjusted read configuration (150-10-24-150bp). A sequencing depth of ~50,000 paired reads per nucleus was targeted. Library preparation and sequencing were performed at the Functional Genomics Center Zurich (FGCZ) of University of Zurich and ETH Zurich.

Single-cell analysis
Pre-processing and filtering. Pre-processing was performed using cellranger-arc 2.0.2 using the optimized mm10 transcriptome from ref. 71. A RNA quantification using alevin-fry72 was also performed to estimate splicing proportions. Decontamination was performed using cellbender 0.3.073 with custom initialization, and doublet identification with scDblFinder 1.19.074. For the RNA modality, scDblFinder was run using standard parameters twice (i.e., with two seeds) and the scores averaged, the AMULET method was used on the ATAC modality. Cells were excluded as doublets if they were called as such using both scDblFinder runs, or if their mean doublet score was above 0.2 and their amulet q-value below 0.01. Only cells with at least 500 fragments (after decontamination), less than 60% contamination, less than 1% mitochondrial reads, more than 20% unspliced fragments (before decontamination) and less than 50% of the fragments assigned to the top 10% features were used for further analysis.

Single-nucleus RNA-seq analysis
Cell annotation was performed based on the RNA modality. A first rough Seurat 5.0.1 clustering75 was performed to split neuronal and non-neuronal cells, which were then processed separately using Seurat.

For neurons, 3500 variable features were selected using the disp method. Prominent activity-dependent genes were removed to limit the influence of the activation status on cell annotation, and known markers were added (see repository for the exact list). Batch correction was performed using harmony 1.2.176. The first 30 dimensions of the harmony embedding were then used to identify the 50 shared nearest neighbors used for clustering at resolution 1.3 and UMAP projection. Cell clusters with a median proportion of unspliced reads below 50%, a median doublet score above 0.3, a median proportion of mitochondrial reads above 0.02, a median proportion of contamination above 0.3, or half of their markers being more highly expressed in empty droplets than in cells were removed as low-quality clusters. Clustering was then performed again on the cleaned object. The clusters were then compared to the Allen brain atlas annotation using SingleR77 for the purpose of annotation, in combination with known markers. In addition, we used the CA1/3 marker genes from Cembrowski et al.46 to distinguish ventral and more dorsal populations.

For non-neuronal cells, we first used geo-sketching78 to 6000 cells as implemented in sketchR 1.1.2 to increase sensitivity to rare cell types. We took the union of the top 2000 variable features identified using the disp and vst methods, as well as the top 1000 most expressed features in the sketched data and known markers. We then performed PCA on the entire dataset based on these features followed by harmony. The first 20 dimensions of the harmony embedding were then used for clustering and UMAP projection. Cell clusters with a median proportion of unspliced reads below 35% or a median doublet score above 0.2 were removed, as well as clusters of cycling cells based on the expression of Top2a, Ccnb2, Cdk1 and Ccnf. The remaining clusters were annotated manually based on known markers.

For ADT-based analysis (below) and visualization, only cells with at least 1000 fragments were used.

Pseudo-bulk differential analysis was performed using edgeR with a ~TimePoint*Treatment model. If the median abundance of the cluster across samples was higher than 50 cells, one surrogate variable estimated using the sva package was further included in the model. For each cell type, genes were filtered using filterByExpr with min.count=20, and only genes with an average log(CPM) above 1.5 were considered. We then harnessed our vast bulk RNA-seq dataset to improve power by estimated global FDR in the pseudobulk data with the IHW package v1.30.0 (10.1038/nmeth.3885), using the FDR from the bulk data as covariate (see in-depth demonstration in ref. 79). Specifically, we used the smallest FDR between the spliced and unspliced bulk data from dropping all experimental covariates (for baseline differences between CRS and Handling, the FDR from the baseline bulk comparison was used instead). Genes which were below filtering thresholds at the bulk level, and which consequently had no attached bulk FDR, were assigned a value of 2 to have their weights estimated separately. IHW was then applied independently for each contrast, across cell types of the same class (i.e., neuronal/non-neuronal). Finally, a minimum absolute logFC of 0.5 was used in addition to FDR to identify differentially-expressed genes.

ADT scoring
The Activity-Dependent Transcription (ADT) score was computed separately for broad cell classes (with neurons grouped as either excitatory or inhibitory). For each class, we first defined a set of stress-induced genes with global (i.e., across time points) FDR below 0.15 (0.05 for non-neurons), an absolute logFC above 0.2 and a log(CPM) above 1. Since the latest time point is chiefly due to GR, we additionally required the gene to have a FDR < 0.15 at either 15min or 45min. For each gene, we then computed the weighted average of its pseudobulk logFC at 15min and 45min, weighted by 1-sqrt(FDR) in the respective timepoint. These logFC estimates then served as the expected activity logFC, on which each cell was regressed. Specifically, we used Gamma-Poisson regression (using the glmGamPoi package 1.17.3 with global overdispersion and no normalization80 of the cell’s counts for those genes on two covariates: the scaled sum expression of the genes across all cells of that type, and the expected logFC. The coefficient assigned to the second covariate then represents a raw activity score, which is however strongly influenced by the proportion of cells of that type that are in an activated state. To address this issue, we then scale the score for each cell, setting 0 to the median score in baseline animals, and 1 to the 98th percentile of the score. An ADT score of 0.5 then represents the point where cells are more similar to the most strongly activated cells than to the median baseline cell. For neurons, an alternative activity scoring was performed using the neuroestimator package version 0.149, and to circumvent the biases of both methods, the two scores were averaged for each cell for the purpose of quantifying the proportions of cells activated.

Single-nucleus ATAC-seq analysis
Using the RNA-based clusters, fragments were merged for neuronal and non-neuronal cell types and peaks were called separately using MACS2 version 2.2.9.181 via the CallPeaks function of Signac 1.13.0 version82. Downstream analysis was performed on pseudobulks of the RNA-based clusters obtained using muscat version 1.16.083. For motif accessibility analysis, peaks were resized to 500bp before running chromVAR on each cell type separately and performing differential analysis as discussed for bulk ATAC-seq. To avoid reporting the several sets of highly similar motifs showing a common pattern, we first applied hierarchical clustering of motifs using the Yule similarity of their matches across peaks, and cutting the dendrogram at a height of 0.5. For each cell type, we took the significant motifs from the same cluster and excluded the ones for which the corresponding gene had an expression below 2 log CPM in that cell type. We then reported the remaining ones in decreasing order of expression.

To assess differences in motif accessibility separately for activated and non-activated neurons (Fig. 8G), we split neurons based on an ADT threshold of 0.5 and compiled pseudobulk accessibility profiles. As activated neurons are fewer, we downsampled the counts of non-activated pseudobulk samples so that their median would match that of activated ones.

Open field test
Open-field testing took place inside sound insulated, ventilated multi-conditioning chambers (TSE Systems Ltd, Germany). The open field arena (45cm × 45cm × 40cm [L×W×H]) consisted of four transparent plexiglas walls and a light gray PVC floor. Animals were tested for 10min under dim lighting (4 lux). Video recordings were collected and distance, time in center, supported rears and unsupported rears were analyzed.

Pose estimation using DLC and behavior analyses
DeepLabCut 2.0.7 (DLC)84,85 was used to track 13 body points of each animal. Tracked points included nose, head center, neck, right ear, left ear, body center, body center left, body center right, left hip, right hip, tail base, tail center and tail tip. The four corners of the OFT were tracked additionally to automatically detect the arena boundaries in each recording. X and Y coordinates of DLC tracking data were imported into RStudio (v3.6.1) and processed with the DLCAnalyzer package86. Furthermore, a previously trained supervised classifier86 was applied to quantify supported and unsupported rears in the OFT on a per-frame basis.

Data availability
The sequencing datasets generated in this study have been deposited to GEO under accession GSE292240. We additionally used previously-published datasets GSE190712 and GSE79790. Other data are provided as Supplementary Information/Source Data. Source data are provided with this paper.

Code availability
All code for analyzing sequencing experiments are available on Github https://github.com/ETHZ-INS/CRS_habituation, as well as on Zenodo https://doi.org/10.5281/zenodo.1746518287.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The study “Distinct molecular mechanisms of stress habituation in the mouse hippocampus” investigated the dynamic changes in RNA expression and genomic accessibility after acute or chronic restraint stress to investigate how the molecular responses to restraint stress change with repeated exposures. The authors leveraged both bulk RNA-seq and ATAC-seq and single nucleus multiome sequencing to identify distinct transcriptomic and genomic signatures of habituation. These thorough studies revealed that repeated stress leads to habituation at a molecular level, largely without baseline gene expression changes or new stress-induced responses. This habituation was mediated by fewer cells being activated by stress, an early cAMP-associated blunting of transcription, and a late GR-dependent mechanism that shortened transcriptional response. However, chromatin accessibility does not explain this molecular habituation. There are also several field-relevant findings such as the similarity between 10 and 20 days of CRS (suggesting 10 is sufficiently “chronic”), and that restraint, stress, swim stress, and Dex produce similar molecular responses. Overall, this manuscript is very well-written, the conclusions are justified, it provides clear insights into molecular mechanisms of stress habituation, and will be a valuable resource for the field, particularly with the interactive web app available by the research team. Minor comments only.\nMain suggestion for publication:\n1. The manuscript proposes two main mechanisms for habituation: a cAMP-associated early mechanism and a later corticosterone-dependent mechanism. While the data supports this, the precise molecular events within these mechanisms remain somewhat vague. Additional discussion or clarification would be useful. In addition, it seems like the blunting of stress response due to fewer cells being activated is a stronger conclusion/better supported than the cAMP and cort hypotheses.\n- For the cAMP mechanism, it’s unclear how the downstream effectors of CREB1 and SRF might lead to the blunted response. Is it decreased expression of these transcription factors, altered phosphorylation, or other modifications?\n- For the corticosterone mechanism, provide more detail on how GR activation is hypothesized to lead to a shortened transcriptional response. Is it active repression of transcription, enhanced degradation of transcripts, or a combination of factors? Is it only in males where CORT peak is lower?\n- Suggestion to include cells-activated hypothesis in abstract/main conclusions\n2. For multiome ATAC, one main assertion is that CRS does not change baseline. It would be valuable to see the baseline CRS vs ARS volcano plots for each cluster as supplement to drive this point home.\n- Authors could tie the lower proportion of activated cells to the snATAC data by dividing the cells into ADT > vs < 0.5, re-pseudobulking snATAC-seq, and comparing GCR and CREB motif accessibility. One would expect motifs to have similar accessibility in activated cells across groups, and similarly closed accessibility in non-activated cells.\n3. It is difficult to glean anything related to habituation from the heatmaps in Figs5-6. These are valuable, but perhaps a “summary” of the meaned ARS vs CRS responses below each heat map would help (as in 6B-C, but across genes shown).\nNice-to-know or minor edits:\n1. In the heatmap in Fig4, are data row-normalized across all 0 and 45m samples or only within the matched time bin, or relative to the group’s own (ARS vs CRS) baseline?\n2. Figure 4 caption for G and H are switched.\n3. For the 10 and 20 day CRS cohort of animals; did the authors also investigate “active” transcripts via looking at unsliced reads as they did for the ARS cohort of animals? It would be interesting to see if CRS perhaps stops active transcription at an earlier time point or does active transcription start at a later timepoint. This would also be interesting to look at in the subset of genes within the “damping patters” in figure 3c. How does the active transcription patters for genes within each of these 3 clusters changes across the timepoints?\n4. Fig 8G: Does this include all genes in the transcriptome or was some threshold applied?\n5. Was there any correlation between cort levels and GR-dependent gene expression in both sexes? This could help strengthen the assertion that that \"in males, a blunting of the HPA-axis response could account – at least in part – for the damping of GR-sensitive genes after repeated stress exposure, whereas in females processes downstream of cort release seem to be at play.”\n(Remarks on code availability)\nThere is a README, I have browsed the code, and I have browsed the web app they made to share the data. I have NOT had time to attempt to run the code, but it appears well-annotated.", "reviewer_3": "This manuscript by Waag et al. represents a tour-de-force study to comprehensively characterize hippocampal transcriptomic changes following repeated restraint stress, integrating both bulk and single-nucleus data, and ATAC-seq under specific conditions, across a wide range of time points and in male and female mice. The authors provide a rich molecular dataset that can serve as a valuable resource for the community. The effort to include both sexes and multiple time points is particularly commendable. However, while the manuscript has many strengths, several aspects require revision to enhance clarity, precision, and alignment between the experimental design and the conclusions. As explained below, a key aspect to review is the reframing o habituation as a transcriptomic phenomenon observed under repeated stress, rather than as a mechanism.\nThe following are specific comments for the authors consideration:\n1. The title and manuscript framing should refer to “habituation to stress” rather than “mechanisms of stress habituation.” The current phrasing gives the impression that the study identifies mechanistic underpinnings, while the work predominantly provides transcriptomic data induced by repeated stress exposure. This change is important for accuracy and appropriate communication of the study’s contributions and should be implemented throughout the manuscript. For instance, the Abstract includes the sentence “We identify two mechanisms that contribute to habituation.” As the study does not demonstrate causality, a more accurate phrasing would be “We identify gene clusters whose expression habituates over time, with distinct temporal dynamics.” Similarly, the phrase “a detailed molecular resource characterizing the dynamic acute stress response and the process of habituation across the genome” should be revised to “a detailed resource characterizing the dynamic molecular acute response to stress in the hippocampus across the genome and how it habituates with repeated exposure to stress.”\n2. Several sentences in the Abstract and Introduction are also too general and overstate the breadth of the study. For example: “It remains unclear whether repeated stress exposure leads to habituation or new adaptive responses.” This statement is overly broad and suggests a focus on behavioral or systemic adaptation, whereas the study is specifically focused on transcriptomic changes in the hippocampus. The language throughout should reflect this specific focus. The same applies to the sentence “We investigated how the molecular response to acute restraint stress changes with repeated exposure,” which should more accurately reference the hippocampal transcriptome. These kinds of reframing are necessary throughout the manuscript to ensure that the interpretation of results is well aligned with the study’s actual scope.\n3. The type of sequencing used should be consistently and correctly described. The Abstract and main text should clarify whether bulk, single-nucleus or single-cell approaches were used. The Abstract could better reflect the methodological novelty of the study. While many previous studies have addressed stress-induced transcriptional changes using bulk RNA-seq, the combination of multi-omic approaches and detailed temporal resolution used here is a distinguishing feature and should be emphasized.\n4. The manuscript refers to “mechanisms that drive different temporal components of stress habituation,” but this again suggests causality that is not demonstrated. The authors cluster genes based on their temporal trajectories and possible GR sensitivity, but the data are not mechanistic.\n5. The study design is comprehensive, but the number of conditions and time points makes the paper difficult to follow. This complexity is partially due to the fact that the acute restraint stress (ARS) rather than being integrated into the chronic restraint stress (CRS) timeline, as point 1, ARS is treated as a separate arm, with a different time structure. Consider restructuring or eliminating these data (i.e., the first results subsection that, in fact, has a separate goal on its own).\n6. Throughout the manuscript, and especially in the Results section, speculative interpretations are presented as if they were empirically supported. For example, the statement “Altogether, these results suggest that glucocorticoids reach the brain and trigger GR translocation to the nucleus within 15min after initiation of stress” is not supported by direct evidence. This interpretation could be correct, but in the absence of brain CORT measurements (microdialysis) or GR translocation assays, it remains speculative and should not be mentioned in the Results. It could be moved to the Discussion, as a hypothesis, not a conclusion. Similarly, the reasoning that “changes in CORT release cannot explain the molecular habituation in females, and they likely cannot fully explain the strong molecular habituation observed in males either” is speculative and inappropriate in the Results section. Further, the rationale to explore chromatin accessibility should not be built on an untested assumption that CORT signaling does not explain the observed effects.\n7. In fact, although the authors measured plasma corticosterone levels at different time points, this study was not done systematically, and it appears presented as a separate analysis, perhaps in a different cohort of animals, not the ones in which the gene expression analyses were undertaken? If that is the case, the authors should be cautious as to where the habituation took place. Given their GR-related findings in a cluster of genes, there is the danger (as not well controlled) that the habituation of the stress response takes place peripherally (as also stated in the Introduction) and that could be on its own responsible, to some extent, of the fact that hippocampal gene expression levels modify/habituate with time, questioning the validity of the single/isolated focus on hippocampal gene expression to explain habituation to stress responses over time – perhaps what it habituates is instead a peripheral trigger of part of those changes. This point should be addressed in the Discussion.\n8. The Discussion includes a section entitled “Is habituation an adaptive or maladaptive process?” which is tangential to the aims and findings of the study. While it is an interesting question, it is not addressed empirically here and feels out of place. If the authors wish to comment on this topic, it could be done very briefly in the concluding lines of the Discussion, but not in a dedicated section.\n9. Finally, the use of the term “healthy” in the Conclusions is confusing. The phrase “These data offer a detailed characterization of the ‘healthy’ stress response” is ambiguous. Chronic restraint stress is known to induce behavioral, structural, and physiological changes. If the authors wish to convey that the model used did not induce overt pathology, that should be explicitly stated. As currently written, the term “healthy” is unclear and risks misinterpretation.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Response to Reviewers’ Comments\nGeneral response to reviewers\nWe thank the reviewers for their valuable feedback, which has prompted us to conduct\nadditional experiments and provide novel analyses. We believe that these new results\nmarkedly increase the robustness and impact of our findings. To help the reviewers evaluate\nthese revisions, we summarize them here in three major points, although we will refer to\nthem again in the individual responses to the reviewers comments.\n1) In response to a suggestion by reviewer 2, we decided to simplify the flow of the paper\nand thus removed the comparison of restraint stress to swim stress at the beginning of\nthe manuscript. As a consequence, we have merged the former figures 1 and 2 into a\nnew figure 1, which focuses only on restraint stress and thus gives the manuscript a\nsharper focus.\n2) Both reviewers remarked on the relationship between transcriptomic habituation and\nchanges in corticosterone (CORT) levels, which had been conducted in separate\ncohorts in the original version of the manuscript. We thus added a new experiment in\nwhich we exposed mice to 10 days of restraint stress and then characterized the\nCORT response and performed bulk RNA-seq in the hippocampus. The experiment\nconfirms that – in the same animals – we simultaneously observe robust transcriptional\nhabituation (new Figure 3J), but an intact CORT response even after 10 days of stress\n(new Figure 3I). This corroborates and strengthens our original conclusion that there is\na disconnect between the systemic HPA axis response, which does not habituate in\nour model, and the transcriptional habituation observed in the hippocampus after\nrepeated stress.\n3) As suggested by reviewer 1, we leveraged our multi-omic single-cell data and our\nactivity-dependent transcription (ADT) score, to separately evaluate motif accessibility\nin activated and non-activated neurons (new Figure 8G). This analysis turned out to be\nextremely informative: First, it revealed a clear increase in the accessibility of\ntranscription factor motifs from early response genes upon ARS (an effect that was\noriginally obscured by the combined analysis). Second, these results demonstrate that\nthe accessibility of the glucocorticoid response element - and thus the transcriptional\nhabituation driven by GR signaling - are independent of neuronal activity. This further\ncorroborates our central claim that two temporally and mechanistically distinct\nprocesses underlie the profound transcriptional habituation that is evident across\nhundreds of stress-responsive genes.\nThe study “Distinct molecular mechanisms of stress habituation in the mouse hippocampus”\ninvestigated the dynamic changes in RNA expression and genomic accessibility after acute\nor chronic restraint stress to investigate how the molecular responses to restraint stress\nchange with repeated exposures. The authors leveraged both bulk RNA-seq and ATAC-seq\nand single nucleus multiome sequencing to identify distinct transcriptomic and genomic\nsignatures of habituation. These thorough studies revealed that repeated stress leads to\nhabituation at a molecular level, largely without baseline gene expression changes or new\nstress-induced responses. This habituation was mediated by fewer cells being activated by\nstress, an early cAMP-associated blunting of transcription, and a late GR-dependent\nmechanism that shortened transcriptional response. However, chromatin accessibility does\nnot explain this molecular habituation. There are also several field-relevant findings such as\nthe similarity between 10 and 20 days of CRS (suggesting 10 is sufficiently “chronic”), and\nthat restraint, stress, swim stress, and Dex produce similar molecular responses. Overall,\nthis manuscript is very well-written, the conclusions are justified, it provides clear insights\ninto molecular mechanisms of stress habituation, and will be a valuable resource for the\nfield, particularly with the interactive web app available by the research team. Minor\ncomments only.\nRE: We thank the reviewer for the positive assessment of our work.\nMain suggestion for publication:\n1. The manuscript proposes two main mechanisms for habituation: a cAMP-associated early\nmechanism and a later corticosterone-dependent mechanism. While the data supports this,\nthe precise molecular events within these mechanisms remain somewhat vague. Additional\ndiscussion or clarification would be useful. In addition, it seems like the blunting of stress\nresponse due to fewer cells being activated is a stronger conclusion/better supported than\nthe cAMP and cort hypotheses.\nRE: We thank the reviewer for this thoughtful comment and we appreciate the opportunity to\nexpand on this important point. In our view, the two mechanisms proposed (cAMP vs GR)\nare complementary to the observation that fewer cells are activated. Indeed, the simplest\nexplanation for the blunting of the early response is that fewer cells are being activated,\nalthough new analyses suggest that even activated CRS neurons show a lower increase in\nCreb1/Srf motif accessibility (newly added Figure 8G-H). Regarding the late, GR-dependent\nresponse, a completely different picture emerges, as this response cannot be explained\nthrough fewer cells being activated, but rather due to a shortened accessibility of GR to the\nDNA. We acknowledge, however, that this was not very clearly conveyed in our manuscript,\nand that the newly added analyses (in part in response to the suggestions of this reviewer)\nhave strengthened and clarified some of these observations. These results are now\npresented in Figure 8G-J and are described in more detail in our answer to question 2.1\nbelow.\n1.1 For the cAMP mechanism, it’s unclear how the downstream effectors of CREB1 and SRF\nmight lead to the blunted response. Is it decreased expression of these transcription factors,\naltered phosphorylation, or other modifications?\nRE: As discussed in the previous response, one simple possibility is that there are no\nchanges in the downstream effectors within single cells, but simply changes in how many\n2\nneurons are activated, which are then mediated by those pathways. The fact that individual\nactivated neurons have a similar ADT score across conditions (Figure 7C) suggests a similar\nactivation of the pathways in activated neurons, although new data (now in Figure 8G-H)\nindicate that there might be differences in the kinetics of their activation. Such changes (in\nCREB1/SRF activity) are certainly not transcriptionally regulated, however it remains\nunknown what kind of post-transcriptional regulation underlies these effects.\n1.2 For the corticosterone mechanism, provide more detail on how GR activation is\nhypothesized to lead to a shortened transcriptional response. Is it active repression of\ntranscription, enhanced degradation of transcripts, or a combination of factors? Is it only in\nmales where CORT peak is lower?\nRE: We conducted an additional experiment to directly compare CORT levels in the blood to\ntranscriptional habituation in the hippocampus (new Figure 3E-I). This experiment confirmed\nthe observation that both males and females are still able to mount a robust CORT response\nafter 10 days of chronic restraint stress. The CORT response in males is slightly weaker\nafter chronic restraint, but significant only at one time point and the peak CORT response is\nidentical. In contrast, in the very same animals, the blunted gene expression response after\nchronic stress is again very strong in both sexes (Figure 3G-I). Further, we observe no\nsignificant correlation between CORT levels and gene expression changes (Supplementary\nFigure 4). Finally, the shortened response of GR-dependent transcripts is very pronounced\nin both sexes. Altogether, we thus think it is very unlikely that the absolute levels of\nperipheral CORT play a critical role in the GR-dependent transcriptional habituation,\nalthough the mechanism remains unknown. Importantly, we see the faster return to baseline\nalso on the level of chromatin GRE accessibility, a fact that we emphasize more in the\nre-written results section (Figure 8A).This argues against a role for RNA degradation, and\noffers new potential mechanisms, although testing these possibilities goes beyond the scope\nof our work. We now added the following passage to the discussion section to address these\npoints:\n“These data suggest that in males, a blunting of the HPA-axis response could account – at\nleast to a small degree – for the blunted expression of GR-sensitive genes after repeated\nstress exposure. However, we think this link is weak for several reasons: 1) The\nCORT-blunting observed in males at a single time-point is a small effect, whereas the\nmolecular habituation is strong and occurs across time points. 2) We observe no clear\ncorrelation between CORT levels and gene expression changes (Supplementary Figure 4).\n3) in females, where the transcriptomic blunting looks very similar to males, we detect no\nchanges in CORT levels after repeated stress exposure. Notably, sex differences in\nglucocorticoid binding globulin have been reported (Galea et al., 1997), thus it remains to\nbe explored if this could contribute to differences in stress habituation. Another possibility,\nsupported by the observation that the initial rise of GRE accessibility is comparable\nbetween groups, but returns to baseline faster in CRS animals, could be that the\nGR-dependent blunting is due to differences in the retention of GR in the nucleus, for\ninstance through nuclear export or what has been proposed as the chaperon cycle\n(Carrigan et al., JBC 2007; George, Schiltz and Hager, 2009; Conway-Campbell et al., Mol\nEndocrino 2011). More work is needed to reveal the mechanisms behind these unexpected\neffects.”\n1.3 Suggestion to include cells-activated hypothesis in abstract/main conclusions\nRE: Thank you, we implemented this change.\n3\n2. For multiome ATAC, one main assertion is that CRS does not change baseline. It would\nbe valuable to see the baseline CRS vs ARS volcano plots for each cluster as supplement to\ndrive this point home.\nRE: We have now added volcano plots comparing gene expression at baseline in the CRS\nvs ARS conditions (line labelled “CRS.baseline”) to Supplementary Figure 9. This\ncorroborates that we do not detect significant differences at baseline.\n2.1 Authors could tie the lower proportion of activated cells to the snATAC data by dividing\nthe cells into ADT > vs < 0.5, re-pseudobulking snATAC-seq, and comparing GCR and\nCREB motif accessibility. One would expect motifs to have similar accessibility in activated\ncells across groups, and similarly closed accessibility in non-activated cells.\nRE: We thank the reviewer for this excellent suggestion. When splitting activated and\nnon-activated (CA1/3 excitatory) neurons based on the ADT score as proposed, we now\nclearly see the expected increase in CREB1 accessibility in activated cells (new Figure 8G,\nreproduced below). This increase was likely overshadowed by a small decrease in CREB1\naccessibility in non-activated cells, which make up the majority of cells. In contrast, as\ncorrectly predicted by the reviewer, GCR motif accessibility was increased upon stress in all\ncells, independently of their activation status (new Figure 8J, reproduced below). Importantly,\nwe again observed that stress-induced GR binding to DNA is shortened after chronic stress,\nbecause the significant reduction in GR accessibility after chronic restraint stress at 45min\nand 3h was observed in both activated and non-activated neurons (Figure 8J). Of note,\nacross conditions and time-points, GCR motif accessibility was generally lower in activated\ncells. This is not GR-specific: a large number of TFs, in particular lineage-specific ones,\nshow an overall lower accessibility in activated cells. This could be due to competition with\nactivity-dependent transcription, and would be in line with similar observations after kainic\nacid treatment (Fernandez-Albert et al, Nat Neurosci, 2019). These results demonstrate that\nthe accessibility of the glucocorticoid response element - and thus the transcriptional\nhabituation driven by GR signaling - are independent of neuronal activity. This further\ncorroborates the central claim that emerges from our work: Two temporally and\nmechanistically distinct processes underlie the profound transcriptional habituation that is\nevident across hundreds of stress-responsive genes in response to repeated stress\nexposure.\n3. It is difficult to glean anything related to habituation from the heatmaps in Figs5-6. These\nare valuable, but perhaps a “summary” of the meaned ARS vs CRS responses below each\nheat map would help (as in 6B-C, but across genes shown).\nRE: We understand the reviewer’s point, and have now added such summary responses in\nFigure 5 (reproduced below). Unfortunately, it would not be meaningful to do so for Figure\n4\n6A, due to the fact that the genes displayed have completely different dynamics. Rather than\nconveying a unified overall pattern, we think of Figure 6A as providing a set of candidate\nCRS-sensitive genes in individual cell types. We hope the reviewer agrees with our\nimplementation of the suggested changes.\nNice-to-know or minor edits:\n1. In the heatmap in Fig4, are data row-normalized across all 0 and 45m samples or only\nwithin the matched time bin, or relative to the group’s own (ARS vs CRS) baseline?\nRE: The normalization is across the entire row, and log2FC are relative to the average of the\nHandling 0min samples. Thank you for flagging this lack of clarity, we have now specified\nthis in the figure legend.\n2. Figure 4 caption for G and H are switched.\nRE: We thank the reviewer and have now corrected this.\n3. For the 10 and 20 day CRS cohort of animals; did the authors also investigate “active”\ntranscripts via looking at unsliced reads as they did for the ARS cohort of animals? It would\nbe interesting to see if CRS perhaps stops active transcription at an earlier time point or\ndoes active transcription start at a later timepoint. This would also be interesting to look at in\nthe subset of genes within the “damping patters” in figure 3c. How does the active\ntranscription patters for genes within each of these 3 clusters changes across the\ntimepoints?\nthe dynamics of the damping patterns, and we now include it in Supplementary Figure 2C\n(reproduced below). Note that such analyses rely on the small proportion of unspliced reads\n(in the order of 10% in such samples), and are thus bound to be noisy. Nevertheless, the\ngeneral trend shows that genes with blunted expression in clusters 2+3 have comparable\ntranscriptional activation at prior timepoints. This supports the hypothesis of a similar\nactivation, but faster return to baseline. We now state this in the main text. However, we\nwould like to point out a few complexities: 1) Not observing a similar increase in active\ntranscription (as is the case for a few genes of these clusters) does not represent evidence\nagainst this hypothesis, as it could equally be interpreted as capturing the moment at which\ntranscription has stopped in CRS, but not handling animals. 2) The faster return to baseline\ncannot be clearly established for all genes, as for some genes the window of active\ntranscription appears to fall between our sampling timepoints – a natural limitation of any\ntime-course sampling. This limitation is particularly important for the genes of cluster 1 (i.e.\nblunted already at 45min) given the lack of an earlier timepoint, and we therefore refrain from\ndrawing a conclusion regarding these genes.\n4. Fig 8G: Does this include all genes in the transcriptome or was some threshold applied?\nRE: Unfortunately we could not determine which figure/panel this question refers to, as there\nwas no Figure 8G in the original submission, nor did this question seem to apply to any other\npanel G in the other figures of the manuscript. (While the question could be relevant to the\noriginal Figure 3G, it was specified in the legend that the plot was based on ‘the union of\ngenes changed after 4h DEX IP or 3h ARS’).\nGenerally speaking, RNAseq analysis uses standard filtering of lowly-expressed genes\n(based on edgeR’s filterByExpr function). To determine DEGs, an absolute logFC threshold\nof at least 0.15 was used in addition to standard FDR thresholding at FDR<0.05.\nIn several instances, comparisons were made subsetting to the relevant genes (for instance\nthe 3h response coefficient relies on the DEGs identified at 3h ARS), however we have tried\nto mention this in the text. Thanks to this comment, we did identify however one place where\n6\nthis was not explicitly mentioned: the comparison of the interaction logFCs between 10 and\n20 days CRS (now Figure 1F) uses the genes with a significant interaction at either\ntimepoint. We have now clarified this in the text.\n5. Was there any correlation between cort levels and GR-dependent gene expression in both\nsexes? This could help strengthen the assertion that that \"in males, a blunting of the\nHPA-axis response could account – at least in part – for the damping of GR-sensitive genes\nafter repeated stress exposure, whereas in females processes downstream of cort release\nseem to be at play.”\nRE: This is an important question, but since we had not performed sequencing and CORT\nanalyses in the same animals, we have now repeated the experiment to address this point.\nMale and female mice were exposed to 10 days of chronic restraint stress (CRS) or daily\nhandling (n=5 mice/treatment/sex), and on day 10 mice from all groups were exposed to\nacute restraint stress (ARS). Throughout the restraint session on day 10, blood was sampled\nfor CORT measurements and 3 hours after stress initiation all mice were killed for brain\nharvest. Ventral hippocampi were collected and processed for bulk RNA-seq. This confirmed\nthat transcriptional damping also happens in the very samples in which we fail to observe a\ndifference in cort activation (new Figure 3G-I). To address the reviewer’s question, we tested\nfor a correlation between CORT levels at every timepoint and either 1) the overall\ntranscriptional response coefficient (at sacrifice, i.e. 3h after ARS initiation) or 2) the\nexpression of Fkbp5 (being one of the most specific targets of GR). Given the substantial\ndifferences in cort levels between sexes, we did this for each sex separately. The results,\nall shown in the new Supplementary Figure 4, show the lack of any statistically significant\ncorrelation (nor was there any significance when combining sexes). We now describe these\nresults in depth in the newly added Figure 3 and in the corresponding text in the results\nsection.\nReviewer #1 (Remarks on code availability):\nThere is a README, I have browsed the code, and I have browsed the web app they made\nto share the data. I have NOT had time to attempt to run the code, but it appears\nwell-annotated.\nRE: Thank you. We routinely use the app in our lab and we integrate user feedback to\nimprove navigation and functions in response to user feedback.\nReviewer #3 (Remarks to the Author):\nThis manuscript by Waag et al. represents a tour-de-force study to comprehensively\ncharacterize hippocampal transcriptomic changes following repeated restraint stress,\nintegrating both bulk and single-nucleus data, and ATAC-seq under specific conditions,\nacross a wide range of time points and in male and female mice. The authors provide a rich\nmolecular dataset that can serve as a valuable resource for the community. The effort to\ninclude both sexes and multiple time points is particularly commendable. However, while the\nmanuscript has many strengths, several aspects require revision to enhance clarity,\n7\nprecision, and alignment between the experimental design and the conclusions. As\nexplained below, a key aspect to review is the reframing o habituation as a transcriptomic\nphenomenon observed under repeated stress, rather than as a mechanism.\nRE: We thank the reviewer for the positive assessment of our work, and for the constructive\nfeedback provided below.\nThe following are specific comments for the authors consideration:\n1. The title and manuscript framing should refer to “habituation to stress” rather than\n“mechanisms of stress habituation.” The current phrasing gives the impression that the study\nidentifies mechanistic underpinnings, while the work predominantly provides transcriptomic\ndata induced by repeated stress exposure. This change is important for accuracy and\nappropriate communication of the study’s contributions and should be implemented\nthroughout the manuscript. For instance, the Abstract includes the sentence “We identify two\nmechanisms that contribute to habituation.” As the study does not demonstrate causality, a\nmore accurate phrasing would be “We identify gene clusters whose expression habituates\nover time, with distinct temporal dynamics.” Similarly, the phrase “a detailed molecular\nresource characterizing the dynamic acute stress response and the process of habituation\nacross the genome” should be revised to “a detailed resource characterizing the dynamic\nmolecular acute response to stress in the hippocampus across the genome and how it\nhabituates with repeated exposure to stress.”\nRE: We appreciate these thoughtful comments and the constructive suggestions on how to\nimprove the clarity of our reporting. We have adopted the last sentence suggested by the\nreviewer, and implemented similar changes throughout the manuscript, to emphasize the\nlargely descriptive nature of our study design and the focus on transcriptional habituation to\nstress in the mouse hippocampus. We also changed the title to emphasize that our study\nfocuses on the transcriptomic changes that occur in the mouse hippocampus as an animal\nhabituates to the repeated exposure to a stressor.\nDuring the revision process we have added important new experimental data and additional\nanalyses, which further strengthen our conclusion that we identified two temporally and\nmechanistically distinct processes that are involved in this transcriptional habituation. Thus\nwe still refer to these distinct mechanisms in the revised title, which we believe is an\nimportant finding which ought to be emphasized. We hope the reviewer agrees with this\ndecision.\n2. Several sentences in the Abstract and Introduction are also too general and overstate the\nbreadth of the study. For example: “It remains unclear whether repeated stress exposure\nleads to habituation or new adaptive responses.” This statement is overly broad and\nsuggests a focus on behavioral or systemic adaptation, whereas the study is specifically\nfocused on transcriptomic changes in the hippocampus. The language throughout should\nreflect this specific focus. The same applies to the sentence “We investigated how the\nmolecular response to acute restraint stress changes with repeated exposure,” which should\nmore accurately reference the hippocampal transcriptome. These kinds of reframing are\nnecessary throughout the manuscript to ensure that the interpretation of results is well\naligned with the study’s actual scope.\nRE: We understand the reviewer’s concern and we have carefully adapted the wording\nthroughout the paper. 1) Both the title and the abstract of the paper now explicitly state that\nour work is conducted in the mouse hippocampus, thus clearly delineating the scope of the\nstudy. 2) We now consistently refer to “transcriptional habituation” and “transcriptional\n8\nblunting” throughout the paper, and we added reminders throughout the introduction and\nresults that the results relate to the mouse hippocampus. 3) The statement “It remains\nunclear whether repeated stress exposure leads to habituation or new adaptive responses.”\nis the second sentence in the abstract, and it is meant to provide the bigger framework in\nwhich our work should be understood. While we kept this sentence, we qualified the specific\nfocus of our work in the next sentence: “Here, we investigated how the transcriptomic\nresponse to acute restraint stress changes – in the mouse hippocampus – with repeated\nexposure.”. We have re-written the rest of the abstract to in response to additional comments\nfrom this reviewer (below), comments from reviewer 1 (about the importance of number of\nactivated cells), and in light of the new findings that emerged from the additional experiments\nconducted during the revision process (a clearer dissociation between early and late\ntranscriptional blunting, which emerge as GR-independent and GR-dependent processes,\nrespectively). We hope that these changes adequately address the reviewer's concerns.\n3. The type of sequencing used should be consistently and correctly described. The Abstract\nand main text should clarify whether bulk, single-nucleus or single-cell approaches were\nused. The Abstract could better reflect the methodological novelty of the study. While many\nprevious studies have addressed stress-induced transcriptional changes using bulk\nRNA-seq, the combination of multi-omic approaches and detailed temporal resolution used\nhere is a distinguishing feature and should be emphasized.\nRE: We appreciate this helpful suggestion, we have re-written the abstract as described in\nour previous response.\n4. The manuscript refers to “mechanisms that drive different temporal components of stress\nhabituation,” but this again suggests causality that is not demonstrated. The authors cluster\ngenes based on their temporal trajectories and possible GR sensitivity, but the data are not\nmechanistic.\nRE: We agree with the reviewer that we do not provide fully elucidated mechanisms,\nhowever we contend that we do establish that two distinct mechanisms are involved – an\nearly GR-independent, and a late GR-dependent one. In addition to the clustering and\nbioinformatic evidence, and the fact that early activity-dependent changes are too rapid to be\nGR-mediated, our study provides the following mechanistic insights:\n1) An injection with the GR-agonist dexamethasone is sufficient to recapitulate the late\ntranscriptomic changes, whereas forskolin-induced neuronal activity does not (Figure 2F,G).\n2) Conversely, forskolin-induced cAMP activation recapitulates the early changes, but not\nthe late GR-dependent changes (Figure 2H).\n3) Single-cell ATAC-seq shows that chromatin accessibility early on (15min) rises similarly in\nmice that experience restraint stress for the first or for the tenth time (in line with the similar\npeak in CORT levels), but at later time points (45min and 3h) GR accessibility drops rapidly\nin chronically stressed mice. This supports a mechanism by which much of the late\ntranscriptional habituation is driven by shortened GR-binding to DNA, while the early\ntranscriptional response is independent of GR.\n4) In the revised version, we have added an analysis in which we split the single-cell\nATAC-seq data based on the activity-dependent transcription score (ADT) from the same\ncells into active and inactive cells (Figure 8G-J, reproduced below). We show a clear\nmechanistic dissociation between cAMP effectors (CREB1 and SRF) and GR: The\ncAMP-signaling cascade implicated in early blunting of gene expression after stress\nhabituation is primarily linked to a reduced number of activated cells. In contrast, the\n9\nshortened transcriptional response of corticosterone-associated genes, which is associated\nwith a protracted blunting of the gene expression response after stress habituation, affects\ncells independently of their activation status.\nIn light of this evidence, we hope that the reviewer agrees with the following change to tone\ndown the sentence in question (last sentence of the introduction) as follows: “[...] we identify\ntwo distinct mechanisms associated with different temporal components of stress\nhabituation.”\n5. The study design is comprehensive, but the number of conditions and time points makes\nthe paper difficult to follow. This complexity is partially due to the fact that the acute restraint\nstress (ARS) rather than being integrated into the chronic restraint stress (CRS) timeline,\nas point 1, ARS is treated as a separate arm, with a different time structure. Consider\nrestructuring or eliminating these data (i.e., the first results subsection that, in fact, has a\nseparate goal on its own).\nRE: We are grateful for this suggestion, it helped us change the manuscript structure and\nsharpen the focus of our study. We had originally started the paper with a comparison\nbetween the transcriptional response triggered by an acute swim vs. acute restraint stress,\nbecause our own legacy data on the acute stress response had primarily focused on swim\nstress, a stressor that does not lend itself to daily repetition (for reasons of animal welfare\nand stark behavioral changes related to floating). However, this section only showed that the\ntranscriptional response to these acute stressors is highly conserved in the mouse\nhippocampus, despite the fact that these stressors are very different. We agree with the\nreviewer that this comparison is distracting and not essential to the main focus of this\nmanuscript. Therefore, we deleted this comparison and only mentioned briefly that the\nresponse to restraint stress is very similar to our previously reported effect of swim stress.\nWe subsequently merged former Figures 1 and 2 into a new Figure 1, which shows the full\nstudy design upfront and thus creates an easier flow (Figure 1A).\n6. Throughout the manuscript, and especially in the Results section, speculative\ninterpretations are presented as if they were empirically supported. For example, the\nstatement “Altogether, these results suggest that glucocorticoids reach the brain and trigger\nGR translocation to the nucleus within 15min after initiation of stress” is not supported by\ndirect evidence. This interpretation could be correct, but in the absence of brain CORT\nmeasurements (microdialysis) or GR translocation assays, it remains speculative and should\nnot be mentioned in the Results. It could be moved to the Discussion, as a hypothesis, not a\nconclusion. Similarly, the reasoning that “changes in CORT release cannot explain the\nmolecular habituation in females, and they likely cannot fully explain the strong molecular\nhabituation observed in males either” is speculative and inappropriate in the Results section.\nFurther, the rationale to explore chromatin accessibility should not be built on an untested\nassumption that CORT signaling does not explain the observed effects.\ndiscussion. However, we believe that the dense presentation of results benefits from small\nsummary statements that help the reader interpret individual sections. To strike a balance\nbetween these needs, we revisited the results section to restrict interpretations to a minimum\nand instead offer concise summaries to guide the reader and ease the transition between\nindividual sections. The two specific examples highlighted by the reviewer have been\ndeleted as part of the general restructuring and newly added data (ATAC-seq motif\naccessibility in both activated and non-activated cell populations, and within-animal\ncorrelations between peripheral CORT levels and hippocampal transcriptional habituation).\n7. In fact, although the authors measured plasma corticosterone levels at different time\npoints, this study was not done systematically, and it appears presented as a separate\nanalysis, perhaps in a different cohort of animals, not the ones in which the gene expression\nanalyses were undertaken? If that is the case, the authors should be cautious as to where\nthe habituation took place. Given their GR-related findings in a cluster of genes, there is\nthe danger (as not well controlled) that the habituation of the stress response takes place\nperipherally (as also stated in the Introduction) and that could be on its own responsible, to\nsome extent, of the fact that hippocampal gene expression levels modify/habituate with time,\nquestioning the validity of the single/isolated focus on hippocampal gene expression to\nexplain habituation to stress responses over time – perhaps what it habituates is instead a\nperipheral trigger of part of those changes. This point should be addressed in the\nDiscussion.\nRE: The reviewer is correct that in the initial experiments, CORT measurements and\ntranscriptional profiling was done in separate cohorts of animals. To address this important\nconcern, we have now repeated this experiment. As described in the “general response to\nreviewers”, these analyses corroborate our initial observation that also after chronic restraint\nexposure the peripheral CORT response is largely maintained and animals still mount the\nsame peak CORT response even after 10 days of chronic restraint, while in the same\nanimals we clearly observe strong transcriptional blunting in the hippocampus. These new\ndata also allowed us to perform correlational analyses between CORT levels and both 1)\nindividual GR-responsive genes, and 2) the overall strength of transcriptomic blunting.\nThese results are described in detail in the newly added Figure 3 and in the corresponding\nresults section. Additionally, as suggested by the reviewer, we have also expanded our\ndiscussion to evaluate the disconnect between peripheral CORT levels and transcriptional\nactivity in the hippocampus.\n8. The Discussion includes a section entitled “Is habituation an adaptive or maladaptive\nprocess?” which is tangential to the aims and findings of the study. While it is an interesting\nquestion, it is not addressed empirically here and feels out of place. If the authors wish to\ncomment on this topic, it could be done very briefly in the concluding lines of the Discussion,\nbut not in a dedicated section.\nRE: We agree with this suggestion and we have removed this section from the discussion.\n9. Finally, the use of the term “healthy” in the Conclusions is confusing. The phrase “These\ndata offer a detailed characterization of the ‘healthy’ stress response” is ambiguous.\nChronic restraint stress is known to induce behavioral, structural, and physiological changes.\nIf the authors wish to convey that the model used did not induce overt pathology, that should\nbe explicitly stated. As currently written, the term “healthy” is unclear and risks\nmisinterpretation.\nagree this was not the optimal choice of words. We instead focus on a more descriptive\nsummary of our findings.\nReviewer #4 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This\nis part of the Nature Communications initiative to facilitate training in peer review and to\nprovide appropriate recognition for Early Career Researchers who co-review manuscripts.\nRE: We thank the reviewer for their contributions to the review process and their feedback.\n12"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have substantially revised the manuscript, including additional experiments and analyses, and have clarified text and added nuance to discussion. The new analyses provide important mechanistic information and increase excitement and impact of this work, which was already set to be a valuable contribution to the field.", "reviewer_3": "The authors have done an excellent job in reviewing this manuscript. It is an important contribution to the stress literature. I do not have any further issues or comments.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": ""}]

Decision letter