Distinct regulation of hippocampal neuroplasticity and ciliary genes by corticosteroid receptors

Glucocorticoid hormones (GCs) — acting through hippocampal mineralocorticoid receptors (MRs) and glucocorticoid receptors (GRs) — are critical to physiological regulation and behavioural adaptation. We conducted genome-wide MR and GR ChIP-seq and Ribo-Zero RNA-seq studies on rat hippocampus to elucidate MR- and GR-regulated genes under circadian variation or acute stress. In a subset of genes, these physiological conditions resulted in enhanced MR and/or GR binding to DNA sequences and associated transcriptional changes. Binding of MR at a substantial number of sites however remained unchanged. MR and GR binding occur at overlapping as well as distinct loci. Moreover, although the GC response element (GRE) was the predominant motif, the transcription factor recognition site composition within MR and GR binding peaks show marked differences. Pathway analysis uncovered that MR and GR regulate a substantial number of genes involved in synaptic/neuro-plasticity, cell morphology and development, behavior, and neuropsychiatric disorders. We find that MR, not GR, is the predominant receptor binding to >50 ciliary genes; and that MR function is linked to neuronal differentiation and ciliogenesis in human fetal neuronal progenitor cells. These results show that hippocampal MRs and GRs constitutively and dynamically regulate genomic activities underpinning neuronal plasticity and behavioral adaptation to changing environments.

G lucocorticoid hormones (GCs) are of critical importance for central nervous system functioning both during daily activities and after a stressful challenge. They play a pivotal role in stress resilience and behavioral adaptation [1][2][3][4][5] . Impaired GC secretion and function has been associated with stress-related mental disorders like major depression and posttraumatic stress disorder 3,5,6 .
GC secretion from the adrenal glands varies over the circadian cycle, peaking at the onset of the active phase. Exposure to a stressful event results in a transient surge of GC secretion that is superimposed on the circadian variation. In the brain, GCs bind to two types of receptors, the mineralocorticoid receptor (MR) and the glucocorticoid receptor (GR) 7 , which show striking differences in localization and function in the brain. Whilst expression of MR is predominantly in neurons within specific limbic structures such as the hippocampus and lateral septum, GR is widely distributed throughout the brain, including the hippocampus 8 . MRs and GRs co-localize in hippocampal neurons 8,9 . Once bound and activated by GCs, MRs and GRs translocate to the cell nucleus to act as ligand-dependent transcription factors. It is thought that genomic actions in the hippocampus underlie the distinct roles of MR and GR in the control of circadian and stress-related physiology, cognition and behaviour, however, the molecular underpinnings of these receptormediated actions are still largely unresolved 10,11 .
Comprehensive knowledge about genome-wide MR and GR interactions within the hippocampus under physiological conditions is currently lacking. Specifically, which genes are regulated by MR and/or GR after stress or across the circadian cycle is still unknown. Previously, we conducted chromatin immunoprecipitation (ChIP) analyses to assess the binding of MR and GR to well-known GC-target genes in the hippocampus of rats killed under baseline and stress conditions 12 . This study showed for the first time that MRs and GRs bind to GC-response elements (GREs) within the GC-target genes Fkbp5 (FK506 binding protein 5), Per1 (Period circadian regulator 1) and Sgk1 (serum/glucocorticoid regulated kinase 1) as homo-or heterodimers in a geneand GRE-dependent manner in vivo 12 . Despite providing novel insight into the complexity of MR/GR interaction, the conclusions of this study were limited to those genes. Therefore, a comprehensive genome-wide study into GC receptor interaction with the DNA was required to identify MR-and GR-regulated genes and clarify the possible existence of selectively MR-or GRregulated genes. Moreover, as long as the identity of all MR-and GR-targeted genes is undetermined, the full implications of GC signaling in the hippocampus will be unknown.
Therefore, in this study, we investigated the interaction of MR and GR with the entire genome of the rat hippocampus after an acute stress exposure (forced swimming (FS)) as well as under baseline early morning (BLAM) and baseline late afternoon (BLPM) conditions. Moreover, we integrated these data with genome-wide gene transcriptional responses occurring in the rat hippocampus under the same conditions. This approach revealed a comprehensive map of the genomic loci bound by MR and/or GR, the transcription factor binding motifs involved, as well as the associated gene transcriptional responses providing novel insights into the function of MR and GR in the hippocampus.

Results
MR and GR ChIP-seq reveals overlapping and distinct features in genome regulation in the hippocampus. Genome-wide changes in hippocampal MR and GR binding after acute stress, or in response to the circadian rise in corticosterone levels, were quantified by ChIP followed by next-generation sequencing (ChIP-seq). MR and GR ChIP-seq was performed on hippocampal chromatin samples (n = 4 independent samples per group) from rats killed under baseline conditions (BLAM, 7.00-9.00 am; BLPM, 5.30-7.00 pm, representing the trough and peak of circadian GC secretion, respectively), or 30 min after the start of an acute, 15-min FS challenge, i.e. at the peak of the stress-evoked plasma corticosterone rise (FS30; Supplementary  Fig. 1a). Regions of the genome with significant peaks of MR or GR binding over input were identified (Supplementary Data 1 & 2) and summarized in Fig. 1. Under BLAM conditions, MRs bound to >300 sites on average across all samples which increased significantly following acute stress (>5-fold, p < 0.05) or at the circadian peak (>7-fold, p < 0.005) (Fig. 1a). In contrast, GRs bound to <100 sites on average under BLAM conditions but this number increased ∼20-fold (p < 0.0001) following acute stress. The circadian peak resulted in GR binding at over 750 sites on average, representing an >8-fold increase (p < 0.005) compared with BLAM GR binding levels (Fig. 1a). On top of each bar, Fig. 1a also provides the number of unique genes showing an MR or GR peak in any of the 4 replicates. Binding of MR and GR to well-characterized target genes, such as Fkbp5, Per1, and Sgk1, was consistent with our previously reported findings (Supplementary Data 1 and 2) 12 . A representative image of MR and GR binding peaks within intron 5 of Fkbp5 is provided in Fig. 1b. MR and GR are clearly binding to the same region within this intron and therefore this location constitutes an overlapping MR/GR peak, which is consistent with previously observed MR-GR coincident binding (as revealed by tandem ChIP) at this genomic site 12 . To determine to which extent MR and GR were binding to the same or distinct loci in each data set (BLAM, FS30, BLPM), Venn diagrams were created to show the number of selective vs. overlapping peaks for each experimental condition. To ensure maximal stringency, these data only included peaks which exclusively bound one receptor in all four replicates but not the other receptor in any replicate ( Analysis of the genomic location of MR and GR binding peaks from all samples (Supplementary Data 1 & 2) revealed the highest levels of binding within distal intergenic regions (MR, 40-50% of peaks; GR, 50-70% of peaks, depending on the experimental condition) (MR, Fig. 1f; GR, Fig. 1g). The remaining peaks were predominantly located in the promoter or intronic regions. Acute stress and circadian influences led to significant increases in the percentage of both MR and GR peaks within intronic regions and a corresponding relative reduction in the percentage of peaks at promoter (MR) or intergenic regions (GR), representing a change in the relative distribution of peaks due to the experimental condition ( Fig. 1f, g). Supplementary Figure 2 provides an overview of the distribution of MR and GR peaks across all chromosomes for each of the four independent replicates and all physiological conditions studied. These peak distribution overviews confirm the increase in MR and GR peaks in the BLPM and FS30 replicates compared with BLAM samples. The increases in MR and GR peaks occur across all chromosomes, including the X-chromosome but not the Y-chromosome where only one peak was found in only a few samples under any condition.
The top-500 peaks, based on statistical significance (false discovery rate, FDR), from the MR and GR datasets were subjected to de novo motif analysis using the Multiple Expectation maximization for Motif Elicitation (MEME) software 13 to identify the binding motif with the highest occurrence and statistical significance. As expected, with high statistical significance, the most prevalent motif was matched to a palindromic GRE (Nr3c1, Nr3c2 and AR motifs) for both MR ( Fig. 1h; E = 1.2e−157) and GR ( Fig. 1i; E = 8.7e−688) binding datasets. Whilst a GRE was present in 100% of the top-500 most significant GR peaks, only ∼67% of the top-500 most significant MR peaks contained a GRE, prompting the question about other non-GRE MR binding motifs (addressed later in this paper).
Genome-wide MR and GR binding is regulated by acute stress and circadian influences. To determine the effect of acute stress and circadian influences on the genome-wide binding of MR and GR, the FS30 and BLPM data were compared with the BLAM data through differential binding analysis, which was performed with the R package DiffBind. This analysis identified genomic regions in which MR binding and GR binding was significantly (FDR ≤ 0.05) upregulated or downregulated in the FS30 or BLPM conditions when compared with the BLAM condition. A correlation heat map of the read count data shows a robust clustering of biological replicate samples obtained under each of the three experimental conditions, as well as showing a clear differentiation between samples collected under BLAM or BLPM conditions, or after stress for both MR (Fig. 2a) and GR (Fig. 2b). DiffBind first used MACS2-called MR or GR peaks to derive a consensus peak set for each transcription factor (MR or GR). A total of 1753 MR peaks and 1066 GR peaks were subsequently analyzed by the DiffBind algorithm and annotated to the nearest gene (1505 genes for MR; 932 genes for GR) using the Bedtools suite. We found that the majority of MR and GR peaks responded to both acute ). This analysis also found that a substantial number of MR peaks did not respond to stress or circadian changes (Fig. 2c, e), whereas this was less the case for GR peaks (Fig. 2d, f). Furthermore, follow-up analyses showed that most of the non-responsive MR peaks (i.e. 490 peaks, 28% of all MR peaks) did not respond to either stress or circadian changes (Fig. 2g). We termed these peaks as being constant. Similarly, we assessed that 59 GR peaks (5.5% of all GR peaks) are constant, thus not responding to either physiological change (Fig. 2h). The constant MR and GR peaks including their gene annotations are listed in Supplementary Data 4 ('CON' tabs). Together, these findings indicate that the interaction of GR with the genome in response to stress and circadian changes is dynamic and corresponds with this receptor's differential GC occupancy pattern under these physiological conditions versus the BLAM situation 7 . Regarding MR, the situation is more complex. Most MR interactions with the genome respond dynamically to stress and circadian changes, which is in line with our previous findings regarding a few classic GC target genes 12 . At the same time, more than a quarter of genomic MR binding is constant. The latter observation dovetails with the long-standing notion that MRs play a key role in mediating constitutive ('tonic' 14 ) effects of GCs. Clearly, however, MRs have a dynamic role as well.
MR and GR peaks encompass distinct patterns of transcription factor binding motifs. We showed that MEME analysis on the top-500 MR and GR peaks had identified the GRE consensus sequence as the predominant transcription factor binding motif (see Fig. 1h Table 1). It also showed that 68-74% of these peaks encompass Krüppel-like Factor (KLF; 68%), Specificity Protein (SP; 70%) and Zinc Finger (ZNF; 74%) motifs. In the MR peaks, a relatively lower incidence of GREs (77%) was observed (Supplementary Table 1). Moreover, we found that the incidence of other motifs such as early growth response protein (EGR; 64%), KLF (70%), SP (74%) and ZNF (76%) was almost similar. These observations support the notion that MRs and GRs, in addition to interacting with GREs, may also interact with other motifs and/or interact with other transcription factors at other motifs.
The high incidence of non-GRE binding motifs for MR prompts the question: if MR or GR peaks do not contain a consensus GRE motif, which would be the predominant binding motifs present in such cases? Among the total of 1753 MR and 1066 GR peaks we found 401 MR peaks and 74 GR peaks that did not encompass a GRE motif. FIMO analysis identified SP (80%), KLF (78%), ZNF (76%), EGR (74%) and RFX (Regulatory Factor X; 73%) as the predominant motifs within the non-GREcontaining MR peaks. With the exception of RFX, these are also amongst the most prominently present motifs if all MR peaks were considered indicating that MRs may interact with such motifs in the presence or absence of a GRE. Within the non-GREcontaining GR peaks, the most frequently observed motifs were Nuclear Receptor (NR; 65%) and Signal Transducer and Activator of Transcription (STAT; 62%). Interestingly, these motifs were not present among the top-5 motifs after FIMO analysis on all GR peaks.
Hippocampal gene transcription changes after acute stress and circadian variation. RNA extracted from the hippocampi of rats killed under baseline conditions (BLAM or BLPM) or at various time points after acute FS stress was sequenced following ribosomal RNA (rRNA) depletion (Ribo-Zero RNA-seq). Supplementary Figure 1b shows the plasma corticosterone responses in these rats which mirror levels obtained in animals used to generate ChIP-Seq samples. After quality control steps and bioinformatics analyses, count tables for both intronic and exonic sequence reads were generated (Supplementary Data 6 and 7, respectively). The majority of intronic reads could be attributed to the nuclear RNA pool and, as such, are a good representation of the heteronuclear RNA (hnRNA) content 15 . hnRNA comprises the unspliced, direct transcript from a gene, therefore its levels provide a more direct index of gene transcriptional activity than the (spliced) mRNA levels, which are a composite measure of RNA synthesis and breakdown. Exonic reads include count data from both mRNA and hnRNA species regardless of polyadenylation status, which is an advantage over the more traditional poly-A+ selection approach 15 . To investigate the effects of acute stress and circadian changes on gene transcription, expression levels at each time point after stress and in the BLPM group were compared with levels at BLAM using differential Fig. 1 Identification of significant MR and GR binding peaks within the rat hippocampus. Genomic regions significantly bound by MR and GR under early baseline conditions (BLAM), 30 min following forced swimming (FS30), or under late baseline conditions (BLPM) were identified. a The average number of significant MR and GR peaks [± SEM, n = 4 biologically independent samples] under BLAM (MR = 345, GR = 86), FS30 (MR = 1849, GR = 1660) and BLPM (MR = 2610, GR = 767) conditions. Above each bar, we depicted the number of unique genes within 5 kb of a peak occurring in any of the four replicates. One-way ANOVA: MR: F(2,9) = 10.06, p = 0.0051; GR: F(2,9) = 59.01, p < 0.0001; *p < 0.05 compared with the respective BLAM group, posthoc Dunnett's multiple comparison test. b Representative images of significant MR and GR peaks within intron 5 of Fkbp5 (FK506 binding protein 5) at BLAM, FS30 and BLPM. Images representative of 4 biologically independent samples. c-e Venn diagrams representing the number of selective (no crossreactivity with opposite receptor in any samples) and overlapping MR and GR peaks at each time point and present in all four biologically independent samples (c BLAM; d FS30; e BLPM). We performed a Fisher's exact test for each Venn diagram and found all overlaps to be highly significant: p = 1.0874e −24 for the BLAM condition and P = 0 for the FS30 and BLPM conditions (zero indicating a value too small to be reported). f, g Genomic location of all MR and GR peaks at each time point represented as % of peaks intersected from all 4 biologically independent samples (Mean ± SEM (n = 4), Two-way ANOVA; f (   To assess experiment quality, we conducted principal component analysis (PCA) using the 1000 most variable genes based on the analysis of the intronic and exonic read counts. Figure 3a-c show the PCA analysis on the intronic data of the BLPM and the FS30 and FS60 time points compared with the BLAM intronic data. In each comparison, the experimental groups separate along PC1 indicating consistent transcriptional changes within groups. PCA of all comparisons including the FS120-, FS180-and FS360-BLAM comparisons are shown in Supplementary Fig. 3. While PCA on intron counts showed strong separation between experimental groups, PCA on the exonic read data only revealed clustering in the BLPM-and, to some extent, the FS30-BLAM comparisons ( Supplementary Fig. 4). This further affirms our hypothesis that intronic reads reflect a more direct measure of transcriptional activity at the measured time points. Using the top 50 differentially expressed genes from each comparison (BLPM, FS30, FS60 vs. BLAM) we obtained a further view on the transcriptional response in these conditions (Fig. 3d). These top differentially expressed genes present a clear distinction between the respective experimental groups and the BLAM group. We can observe the overall similarity between the FS60 and BLPM groups (top segment), a more distinctive BLPM signal (second from the top), genes with high expression at FS30 (second from bottom), as well as genes with high expression at BLAM that decreases in other conditions (bottom segment). Fkbp5, Fgf2 (Fibroblast Growth Factor 2), Zbtb16 (Zinc Finger and BTB Domain Containing 16), Glul (Glutamate-Ammonia Ligase) and Hif3a (Hypoxia Inducible Factor 3a) are some key genes that showed an increased expression both after stress and in BLPM (Fig. 3d). Interestingly, one of the top differentially expressed genes with significantly decreased expression at FS60 is Nr3c1 (which encodes GR); an observation corresponding with our earlier findings 16 .
To compare inRNA expression with exRNA expression across physiological conditions, we compiled the stress time-course data and the BLPM data and present them as constant (CON) RNA levels and responsive (RES) RNA levels to either physiological condition (Fig. 4a, b). Stress-and circadian-induced changes in the RNA expression of individual genes can be found in Supplementary Data 8 and 9. Furthermore, we determined the maximal fold change in intronic and exonic reads (regardless of time point and direction) in response to acute stress for each gene along with its FDR values. These stress-responsive genes were collated and labelled as 'FS merged table' in Supplementary Data 8 and 9. As presented for both intronic (Fig. 4a) and exonic data (Fig. 4b), most genes expressed in the hippocampus showed constant (CON) expression, i.e. no significant change in expression across all physiological conditions (FDR > 0.05; Fig. 4a, b; excluding low expressed genes with less than 10 reads mapped to each sample in the group). We found, however, >4times more responsive (RES) genes showing significant changes (FDR < 0.05) in intronic (cf. exonic) RNA expression (Fig. 4a, b). This difference between differentially expressed genes identified using intronic vs. exonic reads is also illustrated in the volcano plots (Fig. 4c, d; number of genes: intronic reads: 1979; exonic reads: 549).
Figure 4e and f show a direct comparison of the magnitude of the stress-induced versus circadian drive-associated changes identified using intronic and exonic read quantification. This comparison was conducted on those RNAs with changes after both stress and circadian variation. The analysis revealed significant positive (Spearman) correlations for both intronic (r s (470) = 0.9281, p < 0.0001; Fig. 4e) and exonic (r s (171) = 0.8909, p < 0.0001; Fig. 4f) genes, indicating a shared regulatory mechanism for stress and circadian variation. Interestingly, a few genes showed opposing responses to acute stress and the circadian drive (Triangles, Fig. 4e and f) indicating that these genes are likely to be regulated by different mechanisms under these distinct physiological states.
Integration of ChIP-seq and RNA-seq data reveals correlations between MR and GR binding and the magnitude of the gene transcriptional response. Subsequently, we addressed the question: to what extent is MR or GR binding to genes under the different physiological conditions linked to the respective RNA response of such genes in the hippocampus? We found that, of the 1505 MR-binding genes and the 932 GR-binding genes used in the DiffBind analysis (Supplementary Data 4), the expression of 1198 (∼80%) and 718 (∼77%) genes, respectively, could be detected based on intronic and/or exonic RNA count data (Supplementary Data 10). Next, we asked whether these genes would bind both MR and GR, only MR, or only GR. 48% of the genes (i.e. 662 genes) showed MR binding (MR-only), ∼39% (536 genes) both MR and GR binding, and ∼13% (182 genes) GR binding (GR-only). To further scrutinize the specificity of these genes in terms of exclusive MR or GR binding across all four independent replicates, we examined the binding data from individual samples. Of the 662 genes designated as MR-only, 366 were exclusively associated with MR binding, whilst only 28 of the 182 designated as GR-only were exclusively associated with GR binding, i.e. no evidence of opposite receptor binding in any of the individual samples (Supplementary Data 10). Thus, of all RNA expressing, MR/GR binding genes across all three physiological conditions, 366 genes exclusively bind MR (no GR in any sample), 28 genes bind exclusively GR (no MR binding in any sample), and 536 genes bind both MR and GR in all replicate samples.
Subsequently, we determined whether changes in MR and/or GR binding, occurring after stress or during the circadian rise, were correlated with the transcriptional RNA responses in the respective gene. We performed correlation analysis on those genes that showed both a significant change in MR or GR binding in either physiological condition and a significant change in gene expression (in terms of inRNA or exRNA) in these physiological conditions. We found significant, positive correlations between the fold-change in MR binding and the changes in associated inRNA after stress (Fig. 5a) and during the circadian drive (Fig. 5b). Furthermore, a significant, positive correlation was found between changes in GR binding and inRNA during the circadian rise ( Fig. 5d) but not after stress (Fig. 5c). Regarding correlations between exRNA changes and changes in MR or GR binding, we only found significant correlations between MR and GR versus exRNA changes as a result of the circadian rise, but not after stress ( Supplementary  Fig. 5).
Pathway analysis of MR and/or GR bound genes shows individual and overlapping roles for these receptors in the hippocampus. Regarding pathway analysis, we initially used the Gene Ontology (GO) pathway analysis on all RNA-expressing genes with associated MR-only binding (n = 366, Fig. 6a), MR and GR binding (536 genes, Fig. 6b) or GR-only binding (n = 28 genes) (Supplementary Data 11). Given the low number of RNA-expressing GR-only genes, this analysis did not result in any significant pathway associations. In contrast, GO pathway analysis generated an association of MR-only binding with 59 genes involved in cilium structure and function ( Fig. 6a; Supplementary Data 11). Strikingly, most (>70%) of these genes showed constant MR binding across all conditions. In a high proportion (>90%) of these (42) constant genes, MR was binding to genomic regions without a GRE motif but with RFX motifs present (38) (Supplementary Data 11). Such ciliumrelated genes included those encoding transition zone/ciliary gate proteins (e.g. Tmem17, B9d2), centrosomal proteins (e.g. Cep41, Prkar1a), intraflagellar transport proteins (e.g. Ift43, Ift57, Ift122), axonemal dyneins (e.g. Dnah2, Dync2h1) as well as cilia and flagella associated proteins (Cfap; i.e. Cfap20, Cfap298).
In addition to GO enrichment analysis, we used the Ingenuity Pathway Analysis (IPA) suite to analyze both the FS30 and BLPM MR and GR ChIP-Seq data sets, initially, at FDR≤1, thus including all MR and GR bound genes. IPA confirmed the selective association of MR with genes encoding cilia-related structures and (dys)functions (Fig. 6c, d). As expected, due to the constant nature of the majority of these cilia-related genes, the FS30 and BLPM data showed very similar results.
Subsequently, we conducted IPA analysis on the MR-and GRbound genes, which showed significant differential binding (FDR < 0.1) under FS30 and BLPM conditions (Fig. 7). Figure 7 presents various categories grouping a substantial number of genes involved in cell morphology, cell-to-cell signalling, and cell development and movement, as well as genes associated with behavior and neuropsychiatric diseases. Significant associations of MR-and GR-bound genes were found with neuromorphological , d, f) under different physiological conditions compared with BLAM conditions. Differential expression (edgeR) analysis identified genes of which the intronic or exonic expression was upregulated or downregulated (referred to as responsive (RES)) or remained stable (constant (CON)) following acute stress or the circadian rise. a, b Pie charts displaying the proportion of genes of which the intronic (a) and exonic (b) RNA expression is significantly altered in response to stress or circadian influences (RES) or genes of which the expression does not change between conditions (CON). c, d Volcano plots showing the fold-change in RNA expression (log 2 FC) versus the statistical significance of the change (-log 10 FDR). RES genes are shown in red circles whereas CON genes are depicted in blue circles. c Intronic RNA, d exonic RNA). e, f Scatter plots displaying the highly significant Spearman rank correlations between changes in gene expression following stress (FS) and those under circadian influences (BLPM). Intronic (e); n = 470, r s = 0.9295, p < 0.0001; Exonic (f); n = 171, r s = 0.8198, p < 0.0001. Genes responding in the same direction to both stress and the circadian drive are depicted by red circles, whilst genes showed opposing responses to acute stress and the circadian drive are depicted by red triangles (e, Intronic RNA; f, exonic RNA). Source data are provided in the Source Data File.
processes like neurito-and dendritogenesis, inter-neuronal communication like synaptic transmission including long-term potentiation (LTP) and long-term depression (LTD), and neuronal developmental and migratory processes (Fig. 7). In terms of behavior, MR-and GR-associated genes were found to be particularly involved in learning and memory, as well as anxiety-related behaviors. Corresponding with hippocampus function, spatial learning and memory was highlighted. In many categories, a high (>2) Z Score was found indicating a significant activation of these cellular and cognitive functions. IPA also revealed associations of MR-and GR-bound genes with mental, cognitive and neurodegenerative disorders. In view of the significant negative Z Score, it appears that MR-associated genes exert an inhibitory action on neurodegenerative processes (Fig. 7). The analysis of MR-and GR-associated genes revealed involvement in largely parallel functional pathways with the exception of neuronal differentiation, spatial learning, and neurodegeneration. There was considerable similarity between the outcome of the pathway analyses conducted on MR-and GRassociated genes in the FS30 versus the BLPM condition. This observation was supported by the high degree of similarity (86.5-100%) in receptor-associated genes between these physiological conditions as determined by Jaccard analysis (Supplementary Data 12).
Supplementary Data 12 also shows the top-15 canonical pathways linked to MR-and GR-associated genes as analyzed within IPA based on p value and Z score (>2 or <-2). Comparison of MR-and GR-associated pathways presents a substantial level of overlap, but differences as well. For instance, the synaptic LTP and LTD pathways, calcium signalling, Wnt/Ca 2+ (although noncanonical), CREB signalling, and the TEC kinase signalling pathway could be found among both MR-bound and GR-bound data bases of both FS30 and BLPM conditions (Supplementary Data 12). In contrast, AMPK signalling, Rac signalling, and ERK/ MAPK signalling were among the top-15 of MR-bound genes, but not so for the GR-bound genes. Conversely, GR-associated canonical pathway analysis featured micropinocytosis, Sphingosine-1-phosphate, IL-6, and ErbB4 signalling which could not be found among the top-15 list of MR-associated pathway analysis (Supplementary Data 12). Nevertheless, based on this analysis, MR and GR clearly bind to genes participating in inter-neuronal communication (LTP, LTD), neuroplasticity (e.g. nNOS, calcium, CREB, ERK/MAPK, Wnt/Ca 2+ , integrin signalling), cellular metabolism (AMPK signalling), as well as differentiation, growth, migration and survival mechanisms (e.g. TEC, Rac, phospholipase C, sphingosine-1-phosphate signalling). This canonical pathway analysis revealed a number of pathways showing a high degree of overlap between the BLPM and FS30 conditions and all Spearman rank correlation analysis between the fold-change in receptor binding and the associated fold-change in RNA expression under stress and circadian conditions. The results of the correlation analysis are depicted in the subfigures. A significant positive correlation was found between the fold-change in MR binding and the fold-change in inRNA levels after exposure to acute stress (a) and during the circadian rise (b). The change in GR binding post-stress did not significantly correlate with the change in inRNA levels (c). The increased GR binding and the changes in inRNA levels as a result of the circadian rise, however, showed a significant positive correlation (d). Source data are provided in the Source Data File.
indicating an activated status. The overlap between the BLPM and FS30 databases was particularly striking for the pathway analysis of the GR ChIP data: the top-15 list was identical with almost identical p values and Z-scores (Supplementary Data 12). This result suggests that GCs via GR engage the same signalling pathways irrespective of whether GC secretion was brought about by acute stress or the circadian drive.
In Fig. 8, we highlight two typical cilium-related genes and two typical neuroplasticity genes from the MR/GR ChIP-seq and RNA-seq data sets. Their gene maps are shown in Supplementary  Fig. 6. The chosen two cilium MR-selective genes are centrosomal protein 41 (Cep41) 17 and intraflagellar transport 57 (Ift57) 18 ( Fig. 8a-c). The gene promoter of Cep41 showed constant binding by MR in all 12 samples across all three physiological conditions and there was no significant GR binding to this region in any of the samples. Likewise, we found constant MR binding in the gene promoter of Ift57 in 10 out of 12 samples across conditions with no significant binding of GR in any samples. As illustrated in Fig. 8a-c, neither acute stress nor the circadian drive induced a significant change in MR binding in either gene (Fig. 8a) and expression of both intronic and exonic Cep41 and Ift57 RNA was stable under all experimental conditions (Fig. 8b, c).
Protein-tyrosine kinase 2-beta (Ptk2b) 19 and Camk2a 20 are neuroplasticity-related genes that both contain multiple MR and GR binding sites ( Supplementary Fig. 6). Acute stress and BLPM conditions both resulted in significantly increased MR and GR binding to both Ptk2b and Camk2a GREs (Fig. 8d). The increase in MR and GR binding to Ptk2b was associated with increased intronic RNA expression following stress and during BLPM conditions (Fig. 8e). Regarding Camk2a, however, only stress but not BLPM conditions resulted in a significant increase in intronic RNA expression compared with BLAM levels (Fig. 8e). No changes in Ptk2b and Camk2a exRNA counts were found under these conditions (Fig. 8f).
Experimental replication in a distinct cohort of rats validates findings from ChIP-seq and RNA-seq experiments. To further validate the findings from the ChIP-seq and RNA-seq experiments described above, a similar experimental design was applied to an independent cohort of rats with the addition of an extra experimental group which exposed rats to FS at the PM time period. These rats were either killed at 30 min (for ChIP analysis -FS30(PM)) or 60 min (for RNA analysis-FS60(PM)) after the start of stress. Circulating corticosterone levels in BLAM, FS30 and BLPM rats of this experiment were very similar to those in the original sequencing experiments ( Supplementary Fig. 1a-d). Animals that were exposed to FS in the PM showed augmented corticosterone peak levels (at FS30(PM); Supplementary Fig. 1c), most likely due to enhanced HPA drive at this time of day, but relatively lower hormone levels at 60 min post-stress (FS60(PM)) compared with BLPM levels (Supplementary Fig. 1d). The latter finding can be explained by the increased feedback of stressinduced corticosterone at this time of day.
MR and GR ChIP-qPCR analysis of Ptk2b, Camk2a, Ift57 and Cep41 on hippocampus tissue of BLAM, FS30 and BLPM groups of this cohort confirmed our findings in the ChIP-Seq study ( Supplementary Fig. 7). We also included Hif3a in this analysis as one of the highest responding genes in terms of both MR and GR binding and RNA expression. We found significantly increased binding of MR and GR at Ptk2b, There were not enough GR-only bound genes expressing RNA (n = 28) to produce a significant output. GO cellular components associated with MR-only binding were predominantly cilia-related whereas genes binding both receptors were significantly associated with terms relating to synaptic plasticity and spine morphology. c, d Ingenuity Pathway Analysis (IPA) of the FS30 and BLPM MR ChIP-Seq data at FDR≤1. This analysis confirmed the selective association of MR-bound genes with cilia-related structures and (dys)functions. As expected, due to the constant nature of the majority of these cilia-related genes, the FS30 (c) and BLPM (d) associations showed very similar results. The number of genes within the dataset relating to each GO term (a, b) or IPA category (c, d) is provided next to the corresponding bar. Source data are provided in the Source Data File. n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a  Fig. 7a, b, d, e, g, h) and no changes in MR binding to Ift57 and Cep41 ( Supplementary Fig. 7j, m). As expected, GR showed very low enrichment at the cilia genes (B/ I=1-2 ( Supplementary Fig. 7k, n), i.e. similar to IgG negative controls 12 ( Supplementary Fig. 8e-h) but did show a small, but statistically significant increase in binding at both cilia genes due to the circadian rise ( Supplementary Fig. 7k, n). RNA analysis by RT-qPCR of these genes in hippocampus tissues collected under BLAM, FS60 (peak time of responsive hnRNAs under study) and BLPM conditions revealed significant increases in Ptk2b, Camk2a and Hif3a hnRNA levels, but not Ift57 and Cep41 hnRNA levels ( Supplementary Fig. 7c, f, i, I, o). Despite the substantially higher corticosterone levels in rats killed at FS30 in the PM, there was no additional enhancement in MR and GR binding at Ptk2b and Camk2a GREs ( Supplementary  Fig. 7a, b, d, e), possibly indicating attainment of a ceiling effect. The hnRNA levels obtained under these conditions corresponded with the MR and GR binding responses ( Supplementary Fig. 7c, f). Our observations regarding Ptk2b and Camk2a, however, were not universal for all stress-and circadian-responsive MR-and GR-interacting genes. For instance, GR binding, but not MR binding, to Hif3a GRE was significantly higher after FS in the PM than after stress in the AM ( Supplementary Fig. 7g, h), which corresponded with the enhanced stress-evoked Hif3a hnRNA responses in the PM as compared with the AM ( Supplementary  Fig. 7i). Thus, it appears that stress responsiveness of genes over the circadian cycle is gene dependent.
Gene-dependent effects of MR and GR antagonist treatment. Subsequently, we determined if treatment with the MR antagonist spironolactone (Spiro) and the GR antagonist RU486 would affect MR and GR interaction with genes and gene expression. Under BLAM conditions, RU486 increased binding of GR to Ptk2b and Camk2a ( Supplementary Fig. 9b, f), which may be the result of the antagonist-GR complex still being capable of binding to these genes, as previously described in vitro 21,22 . The GR antagonist neither affected stress-induced GR and MR interaction with Ptk2b and Camk2a nor influenced the hnRNA responses of these genes ( Supplementary Fig. 9c, d, g, h). Spiro had no effect on either binding or RNA expression of Ptk2b and Camk2a (Supplementary Fig. 9a-h). In contrast, the antagonist treatments revealed distinct responses in the case of Hif3a verifying efficacy ( Supplementary Fig. 9i-l). RU486 blocked the stress-evoked interaction of both GR and MR with Hif3a as well as the stress-induced hnRNA response of this gene. Spiro treatment led to a decrease in MR binding to Hif3a after stress but did not affect GR binding under this condition. Unexpectedly, Spiro treatment resulted in enhanced FS-induced Hif3a hnRNA responses (Supplementary Fig. 8l) which correspond with the elevated plasma Cort levels in these animals at FS30 (Supplementary Fig. 1e-g) possibly enhancing GR-mediated effects on Hif3a. In addition, it cannot be excluded that MR binding exerts an inhibitory effect at this gene with Spiro administration resulting in diminished MR binding and increased FS-induced hnRNA production. As expected, in view of MR's constant binding to ciliary genes, the antagonists had no effect on MR binding to and hnRNA  responses of Ift57 and Cep41 ( Supplementary Fig. 9m-r). The results regarding Ptk2b, Camk2a and Hif3a show that the impact of antagonist treatment is gene-dependent. This gene dependency may be due to differences in the interaction of antagonistreceptor complexes and (epigenetic) co-factors with the local chromatin, as well as possible compensatory mechanisms involving other transcription factors.
RFX3 binds to the MR-binding site in Cep41 but not to the MR/GR-binding site in Ptk2b. The cilium-related gene Cep41 showed only (constant) MR binding and no GR binding in any sample ( Fig. 8 Figure 9 shows that RFX3 is highly enriched at the 'MR-binding' locus in Cep41 (B/I: 20-30; Fig. 9a) but inconsequential at the 'MR/GRbinding' locus in Ptk2b (B/I: 1-1.5; Fig. 9b). Thus, despite the presence of RFX motifs, RFX3 did not bind to Ptk2b. RFX3 binding to Cep41 was constantly high as it did not change after acute stress (Fig. 9a). These data confirm that at MR-binding loci within ciliary genes, it is not just MR that is binding but RFX3 and/or possibly other RFX factors are binding as well. These observations substantiate a possible joint role of MR and RFX transcription factors in the regulation of ciliary genes.
Critical role of MR in neuronal differentiation and ciliogenesis in human foetal neuronal progenitor cells. We used human foetal neuronal progenitor cells (hfNPCs) to investigate if the MR was playing a functional role in ciliogenesis, an important part of neuronal differentiation. As ciliogenesis has been shown to be linked to neuronal differentiation 27,28 and our GO and IPA analyses revealed a strong link between MR and ciliary genes (Fig. 6), we postulated that hfNPCs would show cilia formation during the neuronal differentiation process. Figure 10a shows that proliferating hfNPCs (as indicated by incorporation of the thymidine derivative EDU (5-Ethynyl-2′-deoxyuridine)) express Nestin, a filament protein expressed in stem cells, but not MRs. After neuronal differentiation, however, the young neurons (stained by Tuj1) express MR protein abundantly (Fig. 10b). Typically, differentiation of hfNPCs is induced by mitogen (EGF/ FGF-b) withdrawal and plating on a suitable substrate. Many years ago, it was empirically determined that the steroid hormone progesterone is required for this process 29 and, nowadays, it is a typical component of commercially available neuronal supplements. Progesterone, however, can act as an agonist or antagonist of the MR 30,31 . Therefore, given that our IPA analysis had revealed a marked role of MR in neuronal differentiation, we postulated that progesterone acts through MRs on this cellular process. To test this postulate, proliferating hfNPCs were cultured using commercial 'N2' supplement (containing 20 µM progesterone) or in home-made supplement of identical composition as the N2 supplement (indicated as 'Prog' in Fig. 10c-e) or with home-made N2 in which progesterone was either omitted (no steroid; 'N.S.') or replaced by the MR agonists, the glucocorticoid hormone corticosterone ('Cort', 0.1 µM) or the mineralocorticoid hormone deoxycorticosterone (DOC, 0.01 µM) (Fig. 10c-e). Figure 10c-e show that administration of Cort or DOC mimicked the effect of progesterone enabling neuronal differentiation, whereas omission of any of these hormones completely blocked this cellular process. Representative images of Tuj1-and DAPIstained neurons cultured under N2, N.S. or DOC conditions are shown in Fig. 10c. The formation of cilia (as visualized by adenylate cyclase 3 (AC3) staining (Fig. 10g)) mirrored the pattern of neuronal differentiation: ciliogenesis occurred in the N2, progesterone, Cort and DOC conditions but not if steroids were omitted (Fig. 10e). These experiments indicate that agonist action via MR is required for neuronal differentiation and ciliogenesis in hfNPCs.
To provide additional evidence for a critical role of MR in neuronal differentiation and ciliogenesis, hfNPCs were differentiated in the presence of commercial N2 supplement and the MR antagonist Spiro. Addition of Spiro (0.1 or 10 µM) resulted in a dose-dependent decrease in nuclear MR staining in the young neurons ( Supplementary Fig. 10). Figure 10f shows representative images of the Tuj1-stained neurons differentiated under Spiro or vehicle conditions. Quantification of neurons and cilia revealed that inclusion of Spiro dose-dependently inhibited neuronal differentiation and cilia formation underlining the requirement of MR action in these processes (Fig. 10h, i).

Discussion
We used a direct approach to study physiologically relevant effects of GCs on the rat hippocampus by delineating the MRand GR-interacting genes and their transcriptional activity under BLAM, BLPM and acute FS stress conditions. Almost all GR peaks responded to the PM and stress conditions whereas only approximately three-quarters of MR peaks were responsive. The other MR peaks remained unchanged. Elaborate pathway analyses clearly indicated functional differences between constant MR-associated genes and stress/circadian-responsive MR/GRassociated genes. We discovered that the constant MR peaks were associated with genes that encode ciliary proteins. Using hfNPCs, we could demonstrate that MRs are indeed functionally involved in neuronal differentiation and ciliogenesis. In contrast, the responsive MR/GR-associated genes are strongly involved in neuroplasticity processes, learning and memory, and several neuropsychiatric disorders. Our analyses also show that there is similarity between the cellular and behavioural processes regulated by MRs and GRs after stress and during the PM phase of the day. Thus, our study provides novel insights into the genomic regulatory mechanisms and pathways engaged by GCs under baseline and stress conditions in the hippocampus. Previous work aiming to gain insight into the effects of GCs on the brain often used indirect, endocrine (e.g. adrenalectomy and GC substitution) approaches. In addition, MR and GR gene deletion strategies have been applied. These studies have been important in contributing to our understanding of cellular and behavioral processes in which GCs could be involved. These approaches, however, have not resulted in comprehensive knowledge about the genes, thus the molecular mechanisms, underpinning the MR-and GR-mediated GC effects and their gene transcriptional implications under physiological conditions. Therefore, we studied these effects through assessment of MR and GR interaction with the hippocampal genome in combination with gene transcriptional responses under physiologically relevant GC conditions (BLAM, BLPM and acute stress). This approach has delivered direct, parallel read-outs of MR-and GR-bound genes as well as RNA expression under these biologically relevant conditions. Although both hippocampal receptors bind corticosterone, our data reveal that, at the genomic level, MRs and GRs behave quite differently. This difference goes beyond their wellknown difference in ligand binding affinity 7,32 . We found that whereas almost all GR peaks respond to the altered physiological condition, this only applied to 72% of MR peaks. Whilst the responsive GR peaks appear to act predominantly via GREs (>95% of stress/circadian-responding peaks contain a GRE), this is less the case with the responsive MR peaks (~70% of these peaks contain a GRE). It should be added that not all GREs within the genome are equally accessible 12 . In general, GRs have been behaving predictably: With rising circulating corticosterone levels due to stress or circadian drive, GR occupancy by ligand increases resulting in increased binding of GRs to accessible GREs within the genome. Regarding MR the situation is more complex. Due to their very high ligand binding affinity, MRs are always highly occupied 7,32 . Nevertheless, despite this high occupancy and their intra-nuclear localization 33 , we showed that MR interaction with the GC target genes Fkbp5, Per1 and Sgk1 is low under BLAM conditions and only increases substantially at BLPM or post-stress 12 . We postulated that the BLPM and stress conditions render GREs and other sites accessible for MR, most likely as a result of epigenetic and other mechanisms 10,12 and is likely to be a factor in the responsiveness of 72% of MR peaks. The other, constantly bound, MR peaks did not respond to BLPM or stress conditions and were less likely to contain GRE motifs. We demonstrated that RFX3 binds the RFX motif in the (constant) MR binding locus in Cep41. This finding indicates that constant MR binding at loci may involve interaction with other transcription factors like RFX3. Moreover, the RFX site within the MR/GR-binding, GRE-containing locus within Ptk2b was not accessible for RFX3. These observations show that control of accessibility of transcription factors to their target motifs in the genome is regulated in a gene-, transcription factor-, and motifdependent manner. Therefore, each transcription factor interaction with the genome is unique and may require subjective experimental verification.
The incidence of MR and GR binding to loci that do not contain a GRE element prompts the question whether these receptors bind directly or indirectly to such non-GRE motifs. An indirect action of MR and/or GR at a non-GRE motif would require the local binding of the intrinsic transcription factor to such motifs. As mentioned, RFX3 binds to the RFX site in the promoter of Cep41, thus MR may be binding in this region through interaction with RFX3. Co-binding of MR or GR can be envisaged at other sites such as EGR, KLF, STAT and FOX as our RNA analyses revealed the induction of the respective transcription factor RNAs EGR1, EGR2, EGR4, KLF2, KLF4, KLF9, KLF13, KLF15, STAT3, STAT5b and FOXO3 after stress and at BLPM (Supplementary Data 8 and 9). Previously, we showed that FS induces EGR1 in distinct hippocampal neurons; a process requiring a confluence of NMDA receptor-mediated ERK1/2, Elk-1, MSK1/2 and (non-genomic) GR signaling, formation of the dual histone mark H3K9ac-S10p at the gene promoter, as well as DNA demethylation of the promoter and 5′-untranslated regions of the Egr1 gene [34][35][36] . Moreover, Chen et al. recently provided evidence that GR and EGR1 interact at the EGR motif within the Bdnf gene 37 . Such interactions have also been reported for GR and STAT3/STAT5b, and GR and KLF2/KLF4/KLF9 at STAT and KLF/GRE binding motifs, respectively, in various cell types [38][39][40][41] . We found that a number of these transcription factor genes (i.e. Klf2, Klf4, Klf7, Klf9, Klf13, Klf15) showed increased MR and GR binding after stress or at BLPM (Supplementary Data 4). Thus, MR/GR action involves direct interaction with GRE motifs and indirect interactions via other transcription factors with their intrinsic binding motifs, whereby GCs can enhance the latter mode of interaction through induction of the respective transcription factor.
In parallel to assessment of genome-wide MR and GR binding, we determined in-and exRNA responses to acute FS stress and the BLPM condition in the rat hippocampus. This approach would reveal if changes in MR and/or GR binding to genes could be linked to changes in transcription of these genes. Previous ChIP-seq studies adopting an endocrine approach did not include assessment of RNA responses and were analysed based on earlier rat genome versions (Rn4 42,43 ). We found significant correlations between circadian-responsive changes in MR and GR binding and the in-and exRNA responses of these MR-and GR-binding genes. Regarding responses to stress, only responsive changes in MR correlated significantly with altered inRNA expression. The higher occurrence of correlations among the inRNA data may be expected as these data are more reflective of dynamic transcriptional activity. The observations regarding the higher correlations in MR and circadian data as compared with GR and stress data are, however, surprising. Apparently, the stress paradigm caused greater variation in the receptor and RNA data than circadian variation. This could be because the stress response is induced by an unpredictable external environment and, therefore, likely to mobilize additional genomic effectors, such as transcription factors other than MR or GR and/or epigenetic mechanisms, which impact on gene transcriptional output. It is likely that multiple transcriptional waves of regulation occur after stress making it hard to correlate the magnitude of all the transcriptional responses back to the initial MR and GR binding events. Conversely, the circadian drive is an endogenous process occurring in a controlled environment, therefore, transcriptional responses may be more predictable. It is unlikely that the optimal time point of the response to stress was missed as we conducted an extensive time-course analysis of the response. Presently, we are only at the beginning of mapping the factors and mechanisms which work alongside MR and GR in the execution of these molecular and physiological processes. Our previous work points to the possible involvement of the transcription factor CREB and the dual histone modification H3K9ac-S10p in acute FS-induced responses 35,44 . Our present results show, however, that there are still many unknowns in the genomic action of MR and GR in physiological responses that require clarification.
The complexity of physiological gene regulation is exemplified in the studies in which we compared FS stress-induced effects between the AM and PM time points or investigated the effects of MR and GR antagonists at the individual gene level. In these studies, we compared MR and GR binding and hnRNA changes in the stress-and circadian responsive genes Ptk2b, Camk2a and Hif3a. In the AM and PM stress experiment, the MR and GR binding pattern to these genes was largely similar as was the hnRNA responses of Ptk2b and Camk2a. The hnRNA response of Hif3a was, however, clearly different as this gene, in contrast to Ptk2b and Camk2a, showed an augmented stress response in the PM, compared with stress effects in the AM. Regarding Hif3a, its PM hnRNA stress response corresponds with the amplified plasma corticosterone and GR to Hif3a GRE binding responses. Such alignment is not observed for Ptk2b and Camk2a. These observations indicate that the hippocampal Hif3a gene is differently regulated than the other two genes. The enhanced responsiveness of Hif3a in the PM may be due to an increase in the number of hippocampal cells in which the Hif3a GRE would be available for GR binding, but apparently not MR binding. According to our data, there seems to be no change in the availability of Ptk2b and Camk2a GREs. Nevertheless, these data show, irrespective of underlying mechanisms, that, at the individual gene level, transcriptional activity is controlled in a genedependent, and most likely cell-dependent, manner. Our data also underline that assessment of a GC-responsive gene is not representative for all GC-responsive genes and plasma GC responses are not a universal predictor of MR/GR binding and RNA responses in GC target genes. The latter notion strengthens our conclusions made in a previous study based on a comparison of different stressors 12 .
The receptor antagonist experiments provided additional evidence that stress-evoked gene regulation is gene dependent. The observation that the GR antagonist RU486 inhibited FS-induced MR binding to Ptk2b, Camk2a and Hif3a seems contradictory. Previously, however, we have shown, that MRs and GRs can heterodimerize at GREs 12 ; therefore, the diminished MR binding after stress in the presence of RU486 may be due to an inability of the RU486-bound GR to heterodimerize with MR 45 . In case of Ptk2b and Camk2a, but not Hif3a, the RU486-bound GR is able to bind to GREs in both the BLAM and FS conditions. As the local chromatin structure and molecular co-factors are known to modulate GR action 46 , the observation of distinct GRE-binding properties of RU486-GR at genes further highlights that such molecular mechanisms are engaged in a highly gene-dependent way. Except for an inhibition of stress-induced MR binding to Hif3a GRE and enhanced Hif3a hnRNA responses, Spiro had little effect on MR/GR binding and hnRNA levels. This may be expected given that hippocampal MRs are always highly occupied by endogenous corticosterone 7,32 preventing any Spiro from receptor binding. The enhanced stress-induced Hif3a hnRNA response in Spiro-treated rats may be due to elevated plasma GC levels in these animals acting through GRs, but this needs to be clarified by future work. Both RU486 and Spiro resulted in increased FS-evoked plasma GC levels, as expected 33,47 , presumably as a result of feedback inhibition. In conjunction, the results of the AM and PM stress experiment and receptor antagonist experiment on selected genes underline the genedependency of GC-regulated transcription.
IPA analysis revealed the activational role (based on >2 Zscores) of stress and PM responsive MR-and GR-associated genes in neuronal morphology, development, and cell-to-cell communication. Genes were found with well-known roles in neuroplasticity and synaptic processes (e.g. dendritogenesis, LTP, LTD) like Ptk2b, Gria1, Grin2b, and Dlgap1; all genes previously unknown to be directly regulated by GC receptors. We also observed an association of responsive MR and GR peaks with genes known to be involved in cognition, learning and memory, and, interestingly, -hippocampus-associated-spatial learning and memory processes. It has been known for many years that GCs play a critical role in hippocampus-associated adaptive behaviors such as observed in the FS test, Morris water maze learning, and contextual fear conditioning 1,4,12,35,[48][49][50][51] . Based on GC effects on (electro-)physiology 52-54 and behavior 1,49,50,55,56 , the notion has been developing that these hormonal effects are the result of molecular and cellular neuroplasticity changes. Presently, however, there is only limited data available regarding the GC-target genes involved in these neuroplasticity changes. Our integrative sequencing methodology reveals potential target genes which may be mediating the link between glucocorticoid action and the consequences for neuroplasticity. Our work provides evidence that GCs affect physiology and behavior through MRs and GRs acting on genes thereby precipitating the molecular and cellular changes underpinning the observed biological changes.
MR-and GR-bound genes also show enrichment of terms related to mental disorders including depressive disorder, bipolar disorder and schizophrenia spectrum disorders (SSDs). It has been known for many years that depressive and bipolar patients often show disturbances in GC secretory patterns. They exhibit impaired negative GC feedback efficacy as shown in the dexamethasone (DEX) suppression test and DEX/corticotropin-releasing hormone (DEX/CRH) test [57][58][59][60] . A disrupted circadian secretion of GCs has been observed as well 61 , particularly in bipolar patients 62,63 , which may be of interest in view of our observation that we only found significant associations in BLPM MR and GR ChIP-Seq data. In fact, corticosteroid receptor dysfunction is thought to be a key factor in the development of major depressive disorder 64 . Hence, hippocampal MR and GR have been shown to be targets of antidepressant drugs 65,66 . A role of GCs in SSDs is less established but evidence has been accumulating indicating a role of GCs both in early life as well as in adulthood in increasing vulnerability for developing schizophrenic symptoms 67,68 . Evidently, IPA unearthed a substantial number of genes, like Disc1 69 , Syn3 70 , Synpo 71 , Grin2a, Gria1, Foxo3, and many others which reportedly have been shown to be linked with depressive disorder, bipolar disorder and/or SSDs. Of key importance is that most of these genes have never before been associated with MR and/or GR interaction, hence were unknown to be potentially regulated by GCs. Therefore, our work provides the opportunity and rationale for the investigation of the genomic basis of GC-evoked vulnerability of these mental disorders.
Both the FS30 and BLPM conditions show enrichment of similar pathways. These include the cellular, behavioral and neuropsychiatric functions, the number of associated genes, and the Z-score output. With few exceptions, there is also considerable parallelism between MR-and GR-associated functions. Our IPA results correspond with the high degree of overlapping MR and GR peaks responding to the FS and BLPM conditions. This outcome is surprising as stress and circadian variation involve very distinct physiological and behavioral responses. Nevertheless, our data show that GCs secreted during these responses act on very similar genes. Therefore, these hormones appear to play a generally supportive or facilitatory role after stress and during the active phase of the circadian cycle: adaptation and learning and memory post-stress, and anticipatory preparation for the active circadian phase. Both situations appear to require a highly similar GC-regulated genomic response. Therefore, based on this outcome, GCs should not be synonymously called 'stress hormones' but rather 'activity-associated hormones'. This designation would correspond with the circumstance that GC secretion occurs during both adverse (physical and psychological stress 72,73 ) and pleasurable (e.g. sex 74 , voluntary exercise 75 ) circumstances, hence circumstances involving physical and/or mental activity.
Our study discloses unique roles for MR in the hippocampus regarding neuronal differentiation and ciliary structure and function. IPA analysis revealed that, both under stress and BLPM conditions, MR, but not GR, has an activational role in neuronal differentiation. Previously, MRs have been reported to exert neuroprotective and pro-survival effects in dentate gyrus neurons in vivo 76,77 . Recently, an additional role for MR in controlling the neuronal phenotype of CA2 neurons was reported 78 . Furthermore, our analyses uncovered an inhibitory role (as indicated by a significantly negative Z-score) of MR-targeted genes in 'neurodegeneration'. These observations agree with a role of MRs in neuroprotection and neuronal differentiation. Moreover, our results provide a deeper understanding of the MR-target genes participating in these effects. Our study also discovered the association of MR with ciliary genes encoding proteins involved in cilium structure (e.g. axoneme, transition zone, centrosome), cilium assembly, and intraflagellar transport. Cilia are hair-like protrusions from the cell body, thought to be involved in neuronal development and communication, but their function is still largely unknown 26,79 . Primary cilia, however, have been shown to be critical for hippocampus-associated functions 80 . Hence, ablation of these organelles from hippocampal neurons through deletion of certain ciliary genes results in severe disruption of dentate gyrus neurogenesis and hippocampus-dependent behaviours (contextual fear conditioning, spatial object recognition) 27,81-83 . Our study shows that, in the majority of cases, MR was constitutively bound to ciliary genes in genomic regions devoid of GREs but with RFX motifs present. We report here that the RFX motif in the MR binding peak within ciliary gene Cep41 indeed binds RFX3 supporting this MR-RFX association. These observations are of interest because: 1. They indicate that MR binding within this gene and possibly other ciliary genes occurs in conjunction with RFX transcription factors, like RFX3; and 2. The RFX motif in Cep41 is accessible for RFX transcription factors (in contrast to the one in Ptk2b) which are known to be of critical importance for cilia structure and function [23][24][25][26] . Future work should investigate how MR and RFX molecules interact at RFX sites within ciliary genes.
Although our work identified an association of MRs with genes involved in neuronal differentiation and structure and function of cilia, a direct functional link had not been demonstrated, neither by us nor by others. Therefore, we pursued the functional role of MR in neuronal differentiation and ciliogenesis in hfNPCs. In their proliferating state, these cells do not express MRs which corresponds with the absence of these receptors in progenitor cells in the subgranular zone of the dentate gyrus 84,85 . Triggering of the differentiation process in our hfNPCs resulted in the expression of MRs in the young neurons. Occupancy and activation of these receptors was a critical necessity as omission of an MR agonist completely blocked neuronal differentiation as well as ciliogenesis. Many years ago, it had been empirically determined that progesterone is a critical constituent of the neuronal supplement ('N2' 29 ). We show that the human MR agonists corticosterone and DOC mimicked the effects of progesterone on the differentiation of hfNPCs into neurons and the expression of cilia in these cells. Moreover, Spiro strongly inhibited the differentiation and ciliogenesis processes further strengthening the functional link between MRs and these processes. Gass et al. 76 have previously shown that mice with MR gene deletion have a decreased number of granular neurons in their dentate gyrus. This was not observed in GR gene deleted animals 76 . Due to the MR deletion, progenitor cells and/or neuroblasts may be unable to differentiate into neurons. Inability to generate cilia may be a critical factor in the differentiation process as deletion or impairment of ciliary genes results in major disturbances in embryonic neurodevelopment and adult neurogenesis 27,28,82,86,87 . Furthermore, dentate granular neurons require agonist-occupied MRs for their viability 88,89 . Together, these MR-dependent mechanisms explain the loss of dentate granular neurons in the MR knock-out mice. Our study provides the identity of the neuronal differentiation-and ciliary genes targeted by MR allowing follow-up studies to gain deeper insight into these mechanisms. Moreover, given the nature of the hfNPCs used, our findings may trigger research on the role of MR in neuronal development during embryogenesis. Our IPA analysis also points to a role of GR, in addition to MR, in maturation, development and morphology (e.g. dendrite-and neuritogenesis) of neurons.
Thus, it appears that MRs play a critical role in the early stages of neuronal differentiation whereas both MRs and GRs orchestrate the morphological development and adaptation of (young) neurons in response to life-imposing activities.
GCs are known to promote resilience, adaptability and cognitive performance to cope with challenges from the environment, but the underpinning genomic mechanisms driving the required physiological and behavioral changes were largely unknown. The present work shows that these hormone actions are mediated by hippocampal MRs and GRs acting at the genomic level on multiple molecular and cellular aspects of neuronal plasticity, in response to an acute stressful challenge as well as anticipatory at the start of the active phase. Our work also reveals a role for GCs acting predominantly and constitutively through MR on ciliary genes which likely involves interactions of this receptor with RFX transcription factors. The functional relevance of the MR-ciliary gene link was observed in differentiating hfNPCs which aligns with the involvement of primary cilia in embryonic and adult neurogenesis. Collectively, our work underlines the multifarious complexity of GC action in the brain and provides the molecular basis for the study of the functional role of these hormones under physiological and pathological conditions.

Methods
All chemicals used below were from Sigma, UK unless otherwise specified.
Animals and stress procedures. All animal work was approved by the University of Bristol Ethical Committee and the Home Office of the United Kingdom (Animal Scientific Procedures Act, 1986, UK).
Adult male Wistar rats (150-175 g on arrival) were purchased from Envigo, UK and group-housed (2-3 per cage). Rats were housed under standard lighting (lights on 05:00-19:00 h, approximately 80-100 Lux) and environmentally controlled conditions (temperature 21 ± 1°C; relative humidity 40-60%) with food and water available ad libitum. Rats were randomly assigned to cages and these were labelled according to a random sequence (https://www.random.org/sequences/) to avoid subjective bias. Due to the nature of the animal work, blinding was not possible at this stage. All rats were handled daily for at least 5 days prior to experimentation to reduce non-specific handling stress on the day of the experiment. Baseline rats were killed straight from their home-cages between either 7-9am (circadian trough, BLAM) or 5-7 pm (circadian peak, BLPM). Alternatively, rats were forced to swim for 15 min in 25°C water, towel dried, returned to their home cage and killed 30, 60, 120, 180 or 360 minutes after the start of swimming (FS30, FS60, FS120, FS180, and FS360, respectively). Sample size was determined based on previous work; n = 4 per group (BLAM, FS30, BLPM) for ChIP experiment (tissue of 2 rats pooled per n; 24 rats in total) and n = 9 per group for RNA (tissue of 1 rat per n; 63 rats in total). Following quick isoflurane anaesthesia (<10 sec) and decapitation, the brain was removed and the hippocampus dissected on an ice-cold steel box. Hippocampus tissues were snap-frozen in liquid N 2 and stored at -80°C. Trunk blood was also collected at the time point of decapitation and stored on ice in tubes containing EDTA and aprotinin for determination of the plasma corticosterone concentration by radioimmunoassay (RIA).
Subsequent validation experiments were performed in the same manner as above but on tissue collected from distinct cohorts of animals. The first experiment was designed to determine if the observed responses to acute stress and the circadian rise were a function of CORT levels or a qualitative difference between the two conditions. Rats were killed under BLAM (8-11am) or BLPM (5-7 pm) conditions or 30 min (for ChIP) or 60 min (for RNA) after FS exposure performed during the AM or PM. Sample size was determined based on previous work 12 and initial sequencing experiments; n = 4 per group for ChIP (tissue of 2 rats pooled per n; 32 rats in total), n = 8 per group for RNA (32 rats in total)). To investigate the effect of MR and GR receptor antagonists, spironolactone (MR antagonist, 50 mg/kg) or RU486 (GR antagonist, 100 mg/kg), respectively, were administered to rats via voluntary dosing which allows the virtually stress-free administration of drugs. Voluntary dosing was conducted by offering rats the vehicle (50% condensed milk, 1 ml) in 1 ml syringes through the bars of the homecage for 5 days prior to day of experimentation. On each day, rats were subsequently handled as described above to reduce non-specific handling stress. On the experimental day, antagonists were mixed into diluted condensed milk/H 2 O (50:50 v/v, 1 ml) and given to rats in the same manner (100% uptake). 60 min after drug (or vehicle) administration rats were subjected to a 15-min FS challenge as detailed above or left in their home cage as non-stress baseline controls. Rats were then killed either 30 min (for ChIP) or 60 min (for RNA) after the start of FS or at an equivalent timepoint (90 min for ChIP; 120 min for RNA) after drug/vehicle administration for non-stress controls. Sample size was determined based on initial sequencing experiments and pilot drug administration tests; n = 4 per group for ChIP (tissue of 2 rats pooled per n; 48 rats in total), n = 6 (for Spiro) -8 (for RU486) per group for RNA (56 rats in total). Finally, given that FIMO analysis identified the RFX binding motif in a high percentage of MR-only binding peaks which lacked a GRE, we performed a validation experiment to determine if RFX3 was also binding in these regions. Groups of rats were killed at BLAM, FS15, FS30, FS60 and FS180 as described above (n = 4 per group (tissue of 2 rats pooled per n); 40 rats in total)). This experiment provided the chromatin to perform the RFX3 ChIP timecourse (Fig. 9).
These experiments were included to strengthen the conclusions drawn from this paper by providing independent biological confirmation of the sequencing results, whilst also providing additional new information regarding the influence of stress during the circadian rise and the effect of specific receptor antagonists. In total, across all experiments included in this publication 295 rats were used. Confounder effects were minimised by experimental design.
Measurement of plasma corticosterone by radioimmunoassay (RIA). Blood samples were centrifuged at 1500 × g and 4°C for 30 min to separate plasma, which was collected into fresh tubes and stored at -80°C until analysis. Plasma corticosterone concentration was determined using a commercial RIA kit (MP Biomedicals) following manufacturers guidelines.
ChIP-Seq sample preparation, library construction, sequencing, mapping and data analysis. Sample preparation -Four independent biological replicates were sequenced per condition; each replicate consisting of the dissected hippocampi from two animals. Frozen hippocampus tissues were cross-linked in a buffered 1% formaldehyde solution containing inhibitors, sonicated and chromatin prepared. ChIP was performed by incubation of antibodies (10 µl/ChIP sample) against either MR (MR H-300, sc11412x, Santa Cruz) or GR (GR H-300, sc8992x, Santa Cruz) with 200 µl chromatin. The specificity of these antibodies for their target has been confirmed by pre-absorption tests and Western analyses previously 12 . Inputs were prepared from 20 µl of the original chromatin sample by reversing crosslinks (NaCl, final conc. 200 nM, at 65°C overnight) and treatment with RNase A (60 µg/ ml, 1 h at 37°C) and proteinase K (250 µg/ml, overnight at 37°C). DNA was purified using Qiagen PCR purification kit as per manufacturer's instructions (Qiagen, Germany).
ChIP-seq library construction and sequencing-DNA was quantified using Qubit (Invitrogen) and the size profile analysed on the 2200 or 4200 TapeStation (Agilent, dsDNA HS Assay). Input material was normalised to 5 ng prior to library preparation. Automated library preparation was performed using the Apollo prep system (Wafergen, PrepX ILMN 32i, 96 sample kit) and standard Illumina multiplexing adapters following manufacturer's protocol up to pre-PCR amplification. Libraries were PCR amplified (18 cycles) on a Tetrad (Bio-Rad) using the NEBNext High-Fidelity 2X PCR Master Mix (NEB) and in-house (Oxford) single indexing primers 91 . Individual libraries were normalised using Qubit, and the size profile was analysed on the 2200 or 4200 TapeStation. Following normalisation, individual libraries were pooled accordingly. The pooled library was diluted to~10 nM for storage. The 10-nM library was denatured and further diluted prior to loading on the sequencer. Samples were pooled and multiplexed libraries were sequenced on 5 lanes of an Illumina HiSeq4000 System to generate 75 basepair (bp) paired-end reads. All samples were assessed for a number of QC metrics, including Q30 and acceptable Passing Filter (%PF). These metrics are in agreement with the values recommended by Illumina for this sequencing machine (HiSeq4000) and read length (75 bp); see https://www. illumina.com/systems/sequencing-platforms/hiseq-3000-4000/specifications.html for more details. In particular, the Q30 value was >96% for read 1 and >92% for read 2 for each of the five sequencing lanes (Recommendation from Illumina: >80%), whereas the %PF was >77% of the raw cluster data in each of the five lanes. Moreover, all ChIP-seq samples in our study had a sequencing depth greatly in excess of the levels recommended for point-source factors in mammalian cells by the ENCODE Consortium (i.e. minimum of 10 million uniquely mapped reads per biological replicate), with an average of >25 million read pairs per sample. All library contruction and sequencing steps were performed by experimenters unaware of group allocation (blinded).
Genome mapping and bioinformatic analysis of ChIP-Seq data. Sequencing data were basecalled with the Illumina software bcl2fastq2 and trimmed to remove any PCR adapters or index primers using Skewer 92 , with settings: -f sanger -x AGATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNN NNNNATCTCGTATGCCGTCTTCTGCTTG -y AGATCGGAAGAGCGTCGTG TAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT -m pe -l 25 -r 0.1 Trimmed reads were aligned using Burrows-Wheeler Aligner (BWA-MEM) 93 to the Rattus norvegicus (rn6) genome build (Ensembl) using default parameters.
Reads with mapping quality < 20, not properly paired, or with insert size > 1000 bp were filtered out using Bamtools 94 . BAM files from individual sequencing lanes were merged and duplicate reads were removed with PICARD (http:// broadinstitute.github.io/picard). The degree of enrichment for the proteins of interest was evaluated through the Normalized Strand Cross-correlation (NSC) and the Relative Strand Cross-correlation (RSC) coefficients, which were calculated for each ChIP sample; in all cases, NSC and RSC values were > 1.05 and >0.8, respectively, in agreement with the ENCODE guidelines 95 . To identify significantly enriched regions, the peak caller Model based Analysis of ChIP-seq (MACS2) was used 96 with parameters -f BAMPE -g 2.5e9 for GR and MR narrow peaks. Peaks with a false discovery rate (FDR) lower than 5% were reported and annotated to their nearest genes (Ensembl release 81) using the Bedtools suite 97 .
The difference in binding levels between conditions -measured by differences in read densities -was performed using the R Bioconductor package DiffBind (version 2.4.8) 98,99 . DiffBind analysis used MACS2-called MR or GR peaks to derive a consensus peak set for each transcription factor (MR or GR). Each consensus peak set contained all peaks that were significantly over input in at least 4 out of 12 samples. The significantly differentially bound sites between sample groups resulting from this analysis were subjected to motif scanning using FIMO (MEME suite, version 4.11.2) 100 and the position weight matrices for vertebrates from the JASPAR CORE database 101 . All the bioinformatic analyses were executed using in-house scripts written in the Python and R programming languages. The pipeline of the bioinformatics analysis of our ChIP-Seq data is depicted in Supplementary Fig. 11. Due to the nature of the bioinformatic analysis blinding was not possible at this stage.
RNA-seq sample preparation, library construction, sequencing, mapping and data analysis. RNA sample preparation-RNA was extracted from nine independent biological samples per group (BLAM, FS30, FS60, FS120, FS180, FS360 and BLPM) using TRI Reagent ® (Sigma, Poole) following manufacturers guidelines. Each biological sample comprised RNA from hippocampi of one rat. Purified RNA was treated with DNase I (Qiagen, Germany) and incubated for 10 min at room temperature. DNase was inactivated by the addition of EDTA (Final conc. 5 mM) for 10 min at 70°C and RNA was purified over MinElute spin column as per manufacturer's instructions (Qiagen, Germany). Three biological samples were randomly selected from the FS60 and FS120 groups to send for initial sequencing to confirm the possibility of conducting Ribo-Zero RNA-Seq. Following the success of this initial test, a further three randomly chosen samples from these groups (FS60 & FS120), along with five randomly selected samples from each of the other condition (BLAM, FS30, FS180, FS360 and BLPM) were subjected to RNA-seq.
RNA-seq library construction and sequencing-RNA was quantified using RiboGreen (Invitrogen) on the FLUOstar OPTIMA plate reader (BMG Labtech) and the size profile and integrity analysed on the 2200 or 4200 TapeStation (Agilent, RNA ScreenTape). RIN (RNA integrity number) estimates for all samples were between 7.4 and 8.5. Input material was normalised to 1 μg prior to library preparation. Total RNA was depleted of ribosomal RNA using Ribo-Zero rRNA Removal Kit (Epicentre/Illumina, Human/Mouse/Rat) following manufacturer's instructions. Library preparation was completed using TruSeq Stranded Total RNA kit (Illumina) following manufacturer's instructions. Libraries were amplified (15 cycles) on a Tetrad (Bio-Rad) using in-house unique dual indexing primers 91 . Individual libraries were normalised using Qubit, and the size profile was analysed on the 2200 or 4200 TapeStation. Individual libraries were normalised and pooled together accordingly. The pooled library was diluted to~10 nM for storage. The 10 nM library was denatured and further diluted prior to loading on the sequencer. Paired-end sequencing was performed using a HiSeq4000 75 bp platform (Illumina, HiSeq 3000/4000 PE Cluster Kit and 150 cycle SBS Kit), generating a raw read count of >45 million reads per sample. All library construction and sequencing steps were performed by experimenters unaware of group allocation (blinded).
Genome mapping and bioinformatic analysis of RNA-Seq data. Quality control analysis was performed on data generated from individual samples as per ChIP (see above). RNA-Seq read pairs were aligned to Rattus norvegicus reference genome, Rn6, using a splice-aware aligner, Hisat2 version-2.0.4 102 . Gene annotation files were downloaded in GTF format from Ensembl, release 81 103 . GTF files with intron coordinates were prepared by subtracting annotated exon features from fulllength transcripts using subtractBed implemented in Bedtools suite v2.26.0. Read fragments mapping to annotated exon features, or intron features were separately quantified to obtain exon, or intron count tables, respectively. All reads were counted with featureCounts 104 , part of subread-v1.5.0 105 , using the specified parameters ("-M" "-O" "-fraction" "-p"). Values for duplication rates and median 3′ bias were estimated using MarkDuplicates.jar, and CollectRnaSeqMetrics.jar implemented in Picard tools v1.92 (https://broadinstitute.github.io/picard/), respectively. Normalized read counts and count based metrics were obtained using in-house R scripts 106 , R core tools, v 3.1.0.
Differentially expressed genes were identified for each comparison using edgeR v3.20.4 107,108 , R core tools, v3.4.2. Each comparison was analyzed separately using subsets of the exon and intron count tables. First, the count tables were selected to samples of interest and expressed genes, keeping genes with greater than 10 reads in each biological replicate in a group. Then, count table subsets were normalized and data were fit using generalized linear models as implemented in edgeR.
Normalized fitted expression values were tested for differential expression between groups using likelihood ratio tests, and finally corrected for multiple comparisons using Benjamini-Hochberg correction 109 . One biological replicate from the BLAM group, RN067, showed a different gene expression profile as compared to other samples in the experiment. We excluded this potential outlier from all downstream analysis as this sample consistently separated from the data and especially the BLAM group in hierarchical clustering as well as principal component analysis (PCA). The pipeline of the bioinformatics analysis of our RNA-Seq data is shown in Supplementary Fig. 12. Due to the nature of the bioinformatic analysis blinding was not possible at this stage.
Pathway analysis. GO pathway analysis-Gene Ontology (GO) pathway analysis was performed using PANTHER (http://geneontology.org/). RNA-expressing gene sets with associated MR-only binding (n = 366, Fig. 6a), MR-and-GR binding (536 genes, Fig. 6b) or GR-only binding (n = 28 genes) were uploaded individually and pathway analysis performed against Rattus norvegicus (all genes in database). Given the low number of RNA-expressing GR-only genes, this analysis did not result in any significant pathway associations. Unmapped genes were manually checked and, where possible, gene IDs were replaced with UniProt IDs to ensure inclusion. Genes which remained unmapped (27 for RNA with MR-only binding; 40 for RNA with MR-and-GR binding) were generally non-coding RNAs and subsequently excluded from pathway analysis.
IPA analysis-Ingenuity Pathway Analysis (IPA; Qiagen) software was used to interpret data from ChIP-seq and RNA-seq experiments. A core analysis was performed on each dataset by IPA which identified canonical pathways, upstream regulators, diseases and functions, and networks predicted to be associated with genes present in the dataset. The settings for all core analyses were as follows: the reference set of genes was taken from the Ingenuity Knowledge base, interaction networks were considered, all node types and data sources were included, experimentally observed predictions were considered, all species were considered, all tissues and cell lines were considered, and all mutations were considered. A ChIP-seq dataset comprising all genes annotated to an MR or GR peak included in DiffBind analysis was uploaded to IPA. For each gene, the fold-change in receptor binding (vs BLAM), p value and false discovery rate (FDR) was included. In instances where multiple peaks of the same category were annotated to a single gene, the largest fold-change value was included in the dataset. In instances where a peak was annotated to multiple possible genes, all gene options were included in IPA. A total of 1996 genes were included in the dataset, with IPA mapping 1495 gene IDs and 501 gene IDs left unmapped. 1303 genes with MR binding and 740 genes with GR binding were analysed. Differentially bound genes at FDR ≤ 0.1 were included in the analysis. An FDR cut-off of 1.0 was also used to include constantly (CON) bound genes.
Core Analysis identified canonical pathways predicted to be influenced by the genes in the datasets based on the directional changes in receptor binding.
Upstream regulators of the genes in the reference datasets were also predicted by interactions by correlating literature-reported effects with genes that had changes in MR/GR binding. Core analysis also identified diseases and functions predicted to be influenced by genes in the reference datasets. A directional change was also predicted on cellular processes and biological functions by correlating observed expression with reported experimental gene effects. Networks were predicted involving non-directional gene interactions.
An 'overlap p value' was calculated using a Right-Tailed Fisher's Exact Test and reflects the likelihood that the association or overlap between the set of significant molecules from the datasets and a given pathway/interaction/network is due to random chance. A smaller p value indicates a lesser likelihood that the association is random. The p value does not consider the directional effect of one molecule on another, or the direction of change of molecules in the dataset.
A 'Z-score' was applied in some analysis types and provided predictions about upstream or downstream processes, considering the directional effect of one molecule on another molecule or on a process, and the direction of change of molecules in the dataset. A Z-score will only have an associated predicted activation state if the Z-score is above 2 or below -2.
Validation and other (RFX) experiments. ChIP-qPCR-Four biologically independent replicates were used per condition for the validation studies described previously (see Animals and Stress procedures), each replicate consisting of the dissected hippocampi from two animals. Chromatin extraction, input preparation and MR/GR ChIPs were performed as described for the ChIP-seq study. Due to discontinuation of all Santa Cruz polyclonal antibodies, validation studies were performed using anti-MR ab64457 or ab97834 (Abcam, UK) and anti-GR 24050-1-AP (Proteintech, UK). This further enhances the validity of our findings in accordance with ENCODE guidelines 95 . These replacement antibodies were checked for target specificity by western blotting (Supplementary Fig. 8a-c). Furthermore, the validation experiments also included IgG ChIPs using recommended IgG isotype control antibodies from the new suppliers (ab171870 (Abcam) and 30000-0-AP (Proteintech)) to determine the non-specific binding ( Supplementary  Fig. 8e-h).
In a separate experiment, ChIP was conducted for RFX3 binding on hippocampal chromatin obtained from rats killed under BLAM conditions or at different times after FS. The RFX3 ChIP was conducted as described for the MR and GR ChIPs using a rabbit polyclonal anti-RFX3 antibody (14784-1-AP, Proteintech; western blot validation in Supplementary Fig. 8d). DNA from input and ChIP samples was quantified using a Qubit 2.0 fluorometer (Invitrogen, UK). All samples were diluted to a standardized concentration with nuclease-free water and analysed using Taqman primers and probes listed in Supplementary Table 2 by qPCR against a standard curve generated from serial dilutions of rat brain genomic DNA (MPBiomedicals, UK). Enrichment of genome binding was expressed as quantity of DNA in ChIP samples (Bound) divided by quantity of DNA in input samples (Input).
RNA-RT-qPCR-RNA was extracted from 8 biologically independent samples per group (see above). Each sample comprised RNA from the hippocampi of one rat. Total RNA (1000 ng) was reverse-transcribed into cDNA using the QuaniTect Reverse Transcription Kit (Qiagen, UK) as per manufacturer's instructions (15 min at 42°C; 5 min at 95°C, BioRad T1000 thermocycler). cDNA was diluted 4-fold and subjected to qPCR analysis using the Taqman primers and probes included in Supplementary Table 2. Results were normalised to expression of two housekeeping genes (Beta-actin (Actb), hypoxanthine phosphoribosyltransferase 1 (Hprt1) and/or tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein, zeta (Ywhaz)) and expressed as fold-change over BLAM levels.
hfNPC tissue collection and culture. Human fetal cortex tissue suspended in phosphate-buffered saline (PBS) was obtained from the SWIFT human foetal tissue bank at Cardiff University (Cardiff, UK) through medical terminations of pregnancy with full donor consent and with ethical approval of the project from the local research ethics committee (SWIFT-RTB 46). PBS was removed and tissue lysed with Accutase ® for 20 min at 37°C to break down into single cells. Single cells were counted using a hemocytometer and seeded at 2 × 10 5 /ml in proliferation media containing high-glucose Dulbecco's Modified Eagle Medium (DMEM):F12 (3:1; Invitrogen), 1% penicillin/streptomycin (Sigma), 1% Gluta-Max (Invitrogen) and supplemented with 2% B-27 serum-free supplement (ThermoFisher), 20 ng/ml fibroblast growth factor (FGF-2), 20 ng/ml epidermal growth factor (EGF) and 5 µg/ml heparin, and grown in a humidified incubator at 37°C in 5% CO 2 . Cells were fed with half the volume taken out and replaced with fresh medium every 4 days until they had grown into neurospheres of~1 mm in diameter. Neurospheres were passaged, returned to proliferation media and incubated as before (37°C, 5% CO 2 ) with regular feeding every 4 days and passaging every 12-16 days from here on. Some cells were incubated with the thymidine analog EDU (10 μM) for 24 h to visualize proliferating cells. After fixation and permeabilization, EDU detection was conducted using the Click-iT EDU imaging kit (Invitrogen, C10337). Immunocytochemistry was performed as described below.
Differentiation of neurospheres. Chamber slides were coated with poly-D-lysine (0.1 mg/ml, Sigma) for 24 h at 4°C, washed thrice with MilliQ water and air dried. 100 µl of laminin solution (0.01 mg/ml laminin, Sigma, in differentiation media, pre-warmed) was added, incubated at 37°C for 1 h before replacement with 100 µl differentiation media (pre-warmed, see below for composition). 2-4 neurospheres were seeded into each well and slides incubated at 37°C for 1 h. An additional 200 µl differentiation media was added to each well after the hour incubation and the cells left to grow in a 37°C incubator with 5% CO 2 for 7 days. Every 2-3 days, cells were fed by removing half the volume out of the wells and replacing with fresh medium. During the differentiation process, cells were transferred to one of the differentiation conditions outlined below. All media were made up in a 3:1 ratio of DMEM:F12 and a volume of 100% ethanol was added to equalize the concentration of solvent in all experimental conditions (final ethanol concentration 0.44%). The composition of the differentiation media was: • Control N2 (N2): 1% N2 supplement 100x (ThermoFisher) containing 0.02 µM progesterone, 5 nM forskolin, 1 mM bovine serum albumin (BSA), 3.68 µM N-acetyl-cysteine (NAC), 1% penicillin/streptomycin (Pen/Strep), 1% Glutamax (Invitrogen). Fixation and immunocytochemistry of proliferating cells and differentiated neurospheres. As described above, proliferating cells were incubated with EDU before fixation and permeabilization and further processing for immunocytochemistry.
Statistical Analysis of correlation data, corticosterone data, and ChIP and RNA validation data. The statistical and graphical package GraphPad Prism 8 was used to generate figures and perform Spearman rank correlation and other statistical analyses (GraphPad Software, San Diego, CA, USA). The independent biological sample sizes (n-values) are indicated in the legends to the figures. The plasma corticosterone, ChIP-qPCR and RNA RT-qPCR data were analyzed with One or Two-way ANOVA followed by the post-hoc Dunnett's (for one-way) or Bonferroni (for two-way) multiple comparison test, as appropriate. Quantification of cilia and neurons, and percentage of neurons over nuclei were statistically analysed using One-way ANOVA and post-hoc Dunnett's test, comparing the mean of each group to the control 'N2' or vehicle group, as appropriate. No statistical analysis was conducted on observations from cells differentiated in N.S. (i.e. the nosteroid condition) as the data did not fit the assumptions of ANOVA (zero cells to count). Testing was always two-sided, as appropriate. P < 0.05 was considered statistically significant.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The data generated in this publication have been deposited in NCBI's Gene Expression Omnibus (GEO) 90 and are accessible through GEO Series accession number GSE126706. Source data on the ChIP-qPCR, RNA RT-qPCR, plasma corticosterone, Tuj1-positive cell numbers and cilia numbers generated in this study are provided in the Source Data File. Source data are provided with this paper.