Introduction

Activity-dependent gene expression is a fundamental mechanism that neural cell types use to respond to experience and coordinate adaptive cellular responses. For canonical immediate-early genes (IEGs) such as Fos, Egr1, and Npas4, transcription begins to rise within minutes and peaks within an hour, highlighting the transient kinetics of these programs^1–6. Large-scale transcriptomic atlases have provided precise molecular definitions of cell types in the brain^7,8, but interpretation of IEG expression within these foundational datasets has been fraught with known artifacts produced during cell capture, a process which itself induces IEGs^9,10. A concerted effort is therefore necessary to characterize the intersection of activity states with a standardized taxonomy of cell types.

The hippocampus is a central structure for learning and episodic memory¹¹, and displays strong cell-type-specific network activity patterns during exploration and rest^12–17. Behavioral states that drive plasticity—such as exploration during arousal and sharp-wave ripples during rest or sleep—predominate at different points in the circadian cycle, positioning the hippocampus as a model system for developing conceptual approaches to the integration of transcriptomic cell types and multidimensional cell states. Prior bulk RNA-seq studies in mice have shown circadian cycling of clock genes in the hippocampus^18,19, but no single-cell data exist to determine how neurons and glia engage these programs. Circadian phase influences excitability^20–23, seizure susceptibility^24,25, and memory performance^18,26, making the interaction between activity-dependent and circadian-dependent transcription biologically important, yet poorly characterized. To date, transcriptional rhythms across the circadian cycle have not been described in the hippocampus using single-cell RNA-seq approaches, leaving open how different cell types entrain their gene expression to this important biological variable.

In this study, we leverage state-of-the-art cell-type resources and a multifactorial experimental design to create an atlas of hippocampal cell states across activity conditions and circadian timepoints, which we call the Activity and Circadian Transcriptomes Dependent on Physiological and Pathological activity (ACT-DEPP) dataset. (Note: Initial analysis of these data was published in a dissertation by J.A.O.²⁷) With this resource, we show that enriched-environment (EE)- and seizure-induced neuronal activity engage qualitatively different transcriptomic programs within all cell types of the hippocampus. We find that within the Allen Brain Cell taxonomy of cell types, there is a collapse of cell type and activity state into a unitary identity among a subset of dentate granule cells. We describe a novel suite of marker genes that are transiently induced by activity in subpopulations of excitatory, inhibitory, and glial cells in the hippocampus. We report a gradient of IEG induction along the dorsal-ventral axis of CA1 in response to EE exposure. Finally, we show here that circadian rhythm coordinates widespread transcriptomic programs in both hippocampal neurons and glia at baseline. By comparing EE exposure to baseline states across the diurnal cycle, we systematically describe how activity- and circadian-dependent transcriptional programs intersect within hippocampal cell types and show that these variables influence endocannabinoid signaling, cell excitability, and chromatin remodeling.

By integrating transcriptional signatures of activity and circadian time into a taxonomy of cell types within the ACT-DEPP dataset, we provide both a conceptual framework for understanding state-dependent gene regulation and a resource for probing how activity patterns shape hippocampal plasticity across the circadian cycle.

Results

Generation of the ACT-DEPP dataset

Single-nucleus RNA sequencing (snRNA-seq) using SPLiT-seq²⁸ was performed on hippocampi from animals exposed to a standard environment (SE); an enriched environment (EE) for 30 minutes or 6 hours; or kainic acid-induced seizures (KA) for 30 minutes or 6 hours. The resulting dataset includes 93,066 nuclei from 48 biological replicates, capturing gene expression in baseline, physiological, and pathological activity conditions across circadian timepoints (Fig. 1A-B, Supp. Fig. 1).

Figure 1.

Quality-controlled nuclei were annotated using MapMyCells and integrated with the Allen Brain Cell (ABC) taxonomy, assigning each nucleus to hierarchical ranks of Class > Subclass > Supertype > Cluster (Fig. 1C inset; ref ⁷). Note: the Subclass rank generally corresponds to conventionally understood cell types within the hippocampus, e.g., dentate granule cells, CA1, PV-expressing basket cells, etc. The UMAP embedding recapitulated relations of canonical hippocampal cell populations at the Subclass level (Fig. 1C), with additional clear transcriptional distinctions for both KA and EE conditions within neuronal and glial Subclasses (Fig. 1D).

To confirm that cell-type annotation was independent of activity-induced transcriptional states, a set of known activity-regulated genes² (ARGs) were removed, and annotations were remapped. Cell-type re-annotations at the Class and Subclass rank were highly stable (Fig. 1E-F), with minimal reassignment even at the lower ABC taxonomic ranks of Subclass and Cluster (Fig. 1G-H). The ACT-DEPP dataset thus robustly registers hippocampal cell types to the ABC taxonomy and reliably distinguishes transcriptional states driven by physiological and pathological stimuli.

The activity-dependent transcriptome is stimulus-specific and obscures cell identity in DG

Seizures have been used extensively to study ARG expression in the past, but direct comparisons to naturalistic stimuli are limited^29,30. Moreover, the near-complete segregation of KA populations within some Subclasses in our uniform manifold approximation and projection (UMAP) plots led us to hypothesize that seizures constitute a categorically different transcriptional state from physiological manipulations. We therefore analyzed differentially expressed genes (DEGs) in CA1 and DG, the most abundant excitatory neuron Subclasses in our samples, at 30 minutes and 6 hours. KA induced an order of magnitude more DEGs than EE at both time points. Approximately one-third of EE-induced DEGs at 30 minutes were distinct from those induced by KA (Fig. 2A). At 6 hours, CA1 responses were fully shared between EE and KA, while DG maintained a unique DEG set specific to EE (Fig. 2B). Thus, novel experiences and seizures trigger categorically distinct transcriptional responses—with respect to both magnitude and specific genes—in these hippocampal subregions.

Figure 2.

To look more closely at the effects of activity on different excitatory-neuron populations, we re-embedded UMAPs and color-coded nuclei by condition or Subclass/Supertype (Fig. 2C). Supertypes were clearly segregated within Subclass clusters, and strikingly, distinct Supertypes emerged as seizure-specific populations. This effect was particularly pronounced in DG, identified by its dramatically altered Supertype distributions between activity conditions (Fig. 2D). In DG, nearly all KA cells clustered within Supertype 138-3, with notable dissociation into Clusters 508-3 (KA 30m; green) and 509-4 (KA 6hr; purple, Fig. 2E). Analysis of Cluster-defining genes revealed six marker genes strongly associated with seizure timing and cluster identity (Fig. 2F-G): Egr2, Egr4, Cenpa, Bhlhe41, Lct, Npy, the latter four not characterized as canonical ARGs. Thus, Supertype DG 138-3 labels recently active granule cells, with additional discrimination of early and late transcriptional responses to activity at the Cluster level. This would likely not be detectable without both EE and KA conditions, given the sparsity of DG activation during naturalistic activity^31–35.

Discovery of cell-type-specific ERGs induced by novel experience

To understand how different cell populations respond to naturalistic activity, we compared early-response gene (ERG) expression across major groups of hippocampal Subclasses (DEGs were categorized as either ERGs or LRGs based on their peak-expression timing; see Methods and ref ¹⁰). Both excitatory and inhibitory neurons share broad upregulation of canonical ERGs such as AP1-, Egr-, Nr4a-, and Dusp-family genes, Arc, and Npas4 across nearly all Subclasses examined (Fig. 3A, E); these genes unsurprisingly exhibit some of the strongest induction genome-wide across nearly all Subclasses, including astrocytes and microglia (Fig. 3B, F, I).

Figure 3.

Next, we were interested in determining which upregulated ERGs exhibited specific activity-dependent induction in one or a few Subclasses. We calculated a tau specificity score previously developed for bulk RNA-seq count data³⁶ to filter for specific changes in expression within EE 30 min vs SE across Subclasses (Methods). When doing so, we observe that canonical ERGs have some of the lowest cell-type specificity (Fig. 3C, G, J). Meanwhile, screening for high-tau ERGs reveals genes whose induction is specific to one or a few related cell types. For example, Hs3st2 and Rcan2 are both genes whose expression is (1) specific to spatially restricted subiculum Subclasses, (2) acutely upregulated following activity, and (3) decaying to baseline expression at 6 hours (Fig. 3D). Likewise, Tac1 is an exemplar activity-dependent marker gene among PV-expressing GABAergic cells and Crhbp in medial-ganglionic-eminence-derived interneurons broadly (Fig. 3H). Finally, we make special note of the astrocytes, surprising in their robust Subclass-specific response to acute EE exposure (Fig. 3I-J). Among the DEGs, Map3k19 and Slco1c1 stand out as genes highly specific to activated astrocytes (Fig. 3K). Thus, comparison of activity-dependent ERG expression following novel experience yields novel marker genes with highly cell-type-specific induction.

Distribution of IEG expression across CA1-ProS Supertypes

Because CA1-ProS neurons were the largest Subclass in our dataset, we used the opportunity to explore experience-induced activity responses within sub-populations of these principal cells. First, we characterized activity in the entire population by adapting a previously used method¹⁰ to classify individual neurons as “active” or “inactive” based on thresholded induction of immediate early genes (Fig. 4A; Methods). Crucially, our library preparation captured only nuclear transcripts before they were exported to the cytoplasm, increasing the likelihood that active cells are classified as such specifically in response to the stimulus^3,9. At baseline (SE), only 8% of cells were active. After 30 minutes of EE, activation increased to 33% of CA1 cells, then declined to 20% by 6 hours. KA seizure robustly activated nearly all cells at 30 minutes, decreasing to about half at 6 hours, with a bimodal distribution possibly reflecting ongoing versus resolved seizures³⁷ (Fig. 4B).

Figure 4.

Transcriptionally defined Supertypes correspond to gradations of anatomically distinct sub-populations of CA1 pyramidal neurons, evidenced by visualization within the ABC MERFISH atlas⁷ (Fig. 4C). We therefore hypothesized that we would observe differences between CA1 Supertypes that may underlie functional specialization within the hippocampus along anatomical axes^38,39. We find that EE strongly activates Supertypes 69-1 and 70-2 compared to others—notable given the dorsal distribution of these Supertypes within hippocampus (Fig. 4C) and the engagement of this part of the hippocampus during exploration^40–42. In contrast, KA seizure manipulations demonstrate that all Supertypes are competent at inducing IEGs following 30 minutes of global activity, although variation in active cells at 6 hours suggests differences in ongoing neural activity, ongoing transcription, or mRNA export among CA1 Supertypes (Fig. 4D-E).

Prior single-cell studies have analyzed IEG expression in different neuronal cell types, including in the hippocampus, to discern whether there are IEG “modules,” i.e. genes that are co-induced together^8,10,30. No such analyses have revealed any interpretable structure of IEG activation. We sought to determine if the unprecedented cell-type resolution among ABC taxa in our dataset would reveal patterns of IEG expression in CA1 subpopulations that were aggregated in prior analyses. If a pattern to stimulus-dependent IEG induction existed, it would emerge in recently active cells, appearing as highly correlated expression modules among IEG subsets, with anti-correlation between subsets. When we tested this hypothesis among active cells in our most abundant Supertypes after 30 minutes of EE, we observed patterns of modest positive correlation among nearly all canonical IEGs across all Supertypes (Fig. 4F). We therefore conclude there is minimal evidence for structured or modular expression of IEGs within subsets of CA1-ProS neurons.

Circadian-dependent transcription among major hippocampal cell types in SE

Biological rhythms entrained to the diurnal light cycle are ancient cellular processes driven by transcription-translation feedback lopps that occur in many tissues, including the mouse hippocampus^23,43. Yet, circadian cycling of transcription remains poorly understood outside the suprachiasmatic nucleus, with no single-cell studies characterizing the extent to which this occurs among different cell types in the hippocampus. We therefore sought to characterize transcriptional changes occurring at four Zeitgeber time points (ZTs) throughout the diurnal cycle under baseline conditions using our snRNA-seq approach, focusing on the four largest Subclasses to reduce noise. Canonical clock genes Per1, Per2, and Cry2 peaked at ZT12 (start of the dark cycle), whereas Clock and Bmal1 had opposite periodicity at ZT0 (Fig. 5A), aligning with previous hippocampal studies^{18,19,23,44–48}.

Figure 5.

We next sought to characterize transcriptional changes taking place genome-wide using DGE analysis. When doing so, we capture the differential expression of canonical clock genes in every subclass tested (Fig. 5B-E, yellow annotations). We note again that Per-family expression peaks at ZT12 in every case but one. In addition to the Per family, its heterodimer partners Cry1, Cry2—as well as transcriptionally correlated Clock-gene regulators Nr1d1, and Nr1d2—peak at ZT12 in Subclasses where they are differentially expressed (Fig. 5B-D). Conversely, although we do not capture differential expression of Clock and Bmal1 in every Subclass, their peak expression in such cases is not concurrent with Per-, Cry-, or Nr1d-family peaks. These phase-shifted expression peaks are functionally consistent with the canonical transcription-translation feedback loop that undergirds circadian rhythms⁴⁹. Thus, the ACT-DEPP dataset captures circadian oscillations of transcription with high fidelity in multiple neuronal and glial cell types.

We next created UMAP embeddings of single-nucleus transcriptomes in circadian-DEG space (2725 genes across all Subclasses) to determine whether there is a transcriptomic structure tracking circadian time. We observe a clear separation of nuclei collected at ZT12 in every Subclass examined, with additional gradations of other timepoints in some cell types (Fig. 5B-D). High tau-specificity analysis confirms statistically what is visually apparent: In all Subclasses, ZT12 is a unique circadian-transcriptomic state with broad upregulation of many genes specifically at that time. Interestingly, in CA1-ProS cells, 12 of 17 of high-tau genes encoded lncRNAs, in contrast to DG granule cells (0/1), astrocytes (4/13) and oligodendrocytes (4/40). This points to a unique recruitment of the non-coding transcriptome in CA1-ProS neurons during this diurnal transition, the functional significance of which is unclear. Finally, among previously characterized ARGs² (ref ⁵⁰ for Adcy genes), several of these varied throughout the circadian cycle, suggesting that there may be an interaction between ZT and activity within the hippocampus.

Interactions between activity- and circadian-dependent transcription in hippocampus

Previous studies have shown differences in canonical ARG expression at different times of day, although these studies could not make direct claims about stimulus-specific ARG induction due to their experimental design^18,18,19. To directly test whether circadian time influences activity-dependent gene expression, we performed EE exposure at different ZTs. We hypothesized, based on these prior studies, that ARG induction would vary in at least some cell types across circadian time; that is, fold-change in ARG expression in EE vs SE would not be constant across the circadian period. We restricted our initial analyses to IEG and clock gene sets used in Figures 4 and 5, respectively. Surprisingly, there were no significant ZT:Activity interactions for IEGs in any of the Subclasses tested (Fig. 6A-D). This was unexpected, given previous reports of circadian-dependent IEG expression at baseline and upstream signaling cascades using non-RNA-seq methods^44,50,51. We nevertheless note strong IEG induction in CA1, irrespective of ZT, unlike other Subclasses.

Figure 6.

Having established that the expression of clock genes varies substantially in baseline conditions (Fig. 5A), we hypothesized that, if present, activity-dependent induction of these genes would not be uniform across ZTs. Indeed, we find that not only are Per1 and Per2 inducible in response to EE exposure, but also that Per1 exhibits a significant ZT:Activity interaction in all Subclasses except DG, with inducibility restricted to times excluding ZT12—notable because ZT12 corresponds to peak expression in SE (Fig. 5A). This observation can be interpreted as a “ceiling effect” of EE-induced expression mimicking the peak circadian expression.

Finally, we tested for ZT:Activity interaction genome-wide. Surprisingly, we only found a small number of genes whose expression was significantly modulated by ZT:Activity interaction, and these genes only appeared in the neuronal Subclasses (Fig. 6I-J). Among the cell types tested, we observe 13 genes that have a statistically significant interaction between ZT:EE exposure, none of which are canonical IEGs. Functional annotation of these genes reveals that two are involved in endocannabinoid synthesis (Daglb, Napepld), three in cell signaling or excitability (Fcho2, Itpkb, Kcnn2), and four in chromatin remodeling or gene expression (Cecr2, Zfp420, Zpf523). We therefore conclude that changes in gene expression following novel experience are relatively constant across the circadian cycle, although special mechanisms may be in place to selectively change induction mechanisms for genes which modulate endocannabinoid synthesis and cell excitability within neurons of the hippocampus.

Discussion

The ACT-DEPP dataset is an atlas of single-nucleus transcriptomes from the hippocampus, capturing transcriptional responses to both physiological and pathological brain activity. The creation of this atlas accomplishes two goals: (1) mapping of experience-dependent gene expression in a high-resolution cell-type taxonomy, and (2) characterization of circadian-dependent transcriptional oscillations in both neurons and glia of the hippocampus. All cell types are registered to a gold-standard brain-wide taxonomy, providing rich annotation that we leverage to demonstrate that stimulus-responsive cell states are embedded within broader cell-type identities, such as hippocampal sub-region (Fig. 3) and anatomical location (Fig. 4), with DG serving as an exception to this rule (Fig. 2). We identify the light-to-dark transition as a key inflection point for diurnal gene regulation in the hippocampus (Fig. 5) and show that activity-dependent induction does not vary widely across the circadian period (Fig. 6). We envision this resource as a roadmap for analyzing cell-type-specific, stimulus-responsive states in the brain using modern transcriptome-registration tools. More broadly, the ACT-DEPP dataset is a resource for understanding how physiological and pathological activity patterns drive experience-dependent transcriptional programs to support neural plasticity across the circadian cycle.

Limitations

We take care to point out a potential conceptual limitation of our study: Dissociating neuronal gene expression driven by circadian phase from that driven by acute locomotion is inherently challenging because the two are tightly coupled in behavior and physiology. Locomotor activity co-occurs with changes in arousal, hormones, and neuromodulator release, many of which independently influence transcription, while the propensity for locomotor activity varies across the circadian cycle^23,52. One argument that favors the distinction we make between activity- and circadian-dependent expression is that canonical IEGs do not overlap with circadian DEGs (Fig. 5), suggesting orthogonal regulatory logic between activity programs and circadian programs. Second, in contrast to neurons, glial transcription is less sensitive to rapid fluctuations in circuit activity and more influenced by systemic cues, making their rhythmic gene expression less confounded by locomotor events. Taken together, these observations support the validity of interpreting our findings as reflecting both activity- and circadian-driven transcriptional programs, rather than a single merged signal.

We also disclaim notable variance in cell number and read depth among a small number of our biological replicates (Supp. Fig. 2 and 3). Although our pseudobulk approach accounts for both of these factors by computing size factors for gene-count normalization before performing DGE analysis, it is nevertheless likely that estimates of variance are noisier in these conditions. However, for all of our DGE analyses, we set DESeq2 gene-filtering parameters that would bias our analysis toward excluding genes from DEG lists based on unacceptably high levels of within-condition variance for that gene model. Pseudobulk DGE analysis has been shown to be the most robust against Type I errors for single-cell data while performing the best in balancing Type II error rate^53,54. Thus, the results reported in this paper should be protected from unacceptably high false-positive rates.

Cell states and cell types

Beyond our investigation of circadian transcription, our study has three main contributions to our understanding of activity-dependent gene expression. First, we extensively leverage registration of our dataset with the ABC atlas—a state-of-the-art ontology of cell types in the mouse brain⁷—to query activity-dependent states among high-resolution cell types. The ability to observe the effects of activity manipulations at the Supertype and Cluster level allows us to infer activity distributions along functionally relevant anatomical axes in the hippocampus (Fig. 4). Second, we distinguish several candidate genes that are highly specific not only to cell types, but also activity states (Fig. 3). These genes may serve as guideposts for the development of simpler next-generation tools to study activity-dependent processes with cell-type granularity going beyond IEG-based permanent labeling systems using Fos, Arc, or Npas4^55–58—all of which are broadly upregulated across multiple cell types in our data, including non-neuronal cell types (Fig. 3). Such tools would have the advantage of being able to permanently access task-specific ensembles within specific cell types without the need for complex intersectional genetic strategies. Third, our descriptions of activity in baseline and “ceiling” conditions such as seizure provide strong evidence for the annotation of specific Supertypes and Clusters as highly temporally specific signatures of recent activity (Fig. 2).

Within DG granule cells, we observed a collapse of cell type and activity state at the supertype level (Fig. 2), and we also demonstrate that the Cluster rank in this population further resolves activity states into both acute and subacute activity histories, on the order of 30 minutes to 6 hours within our dataset (Fig. 2E). This example case in DG begs the question: Does the annotation of activity state—or any other condition-specific state, e.g. cellular manipulation or disease—at the Supertype or Cluster rank merit the collapse of these taxa with closely related inactive populations within the ABC taxonomy? More broadly, to what extent should reference atlases balance heterogeneity of cell states with the classification of cell types? On the one hand, the development of genetic tools to interrogate a particular Supertype or Cluster may produce unexpectedv results if that taxon represents a distinct cell state not present, or overrepresented, within an experimental dataset. On the other hand, it is easy to imagine how the annotation of active DG cells, for example, through a simple compositional analysis of Supertypes or Clusters would be useful for high-throughput readout of activity. It is our opinion that the correlation of specific cell states with lower ranks in the ABC atlas is acceptable, as long as these cell states are ascertained using robust experimental methods whose results are clearly disclosed to the community.

Circadian rhythm changes baseline gene expression in the hippocampus

It has long been established that cells in the suprachiasmatic nucleus (SCN) use transcription-translation feedback loops entrained to ambient-light sensation to keep a “molecular clock” that tracks the diurnal cycle (Patton & Hastings 2018; Patton 2023). However, there are few investigations of hippocampal circadian-dependent transcription in the omics era, and no snRNA-seq studies that we are aware of. Therefore, we sought to characterize how transcription changes in the hippocampus across the circadian cycle using snRNA-seq. Our results extend previous work by confirming the presence of clock-gene cycling in hippocampal neurons and glia across the circadian cycle. The expression patterns of canonical clock genes we observe are consistent with earlier bulk hippocampal studies^18,19,23,44, lending confidence that our sampling strategy captures genuine circadian processes rather than stochastic variation. Importantly, we find that ZT12 represents a distinct transcriptional inflection point across multiple cell types, coinciding with the transition from the light to dark period. This time is behaviorally salient for nocturnal animals such as mice, as locomotor activity begins to rise sharply before and during the onset of darkness^23,47,52,59. Our data therefore support the interpretation that circadian gene expression in hippocampal cells is tightly coupled to this environmental and behavioral transition.

The similarity of transcriptomic profiles at ZT16 and ZT0 further suggests that circadian-driven transcription is concentrated within a narrow temporal window. Rather than a gradual progression across the cycle, many genes appear to undergo sharp, temporally restricted bursts of expression near the light–dark transition. This “bursty” model of circadian regulation is consistent with prior reports of pulse-like transcription in single cells^60–63 and may reflect a strategy by which hippocampal cells track salient transitions without maintaining high levels of rhythmic transcription throughout the entire dark period.

Among the cell types tested, we observed a striking specificity in circadian-dependent transcription: oligodendrocytes, which exhibit relatively few activity-dependent DEGs (Fig. 3), display the largest number of circadian DEGs. This dissociation underscores that activity-dependent and circadian-dependent transcription represent partially non-overlapping regulatory axes. For oligodendrocytes, the high number of circadian DEGs suggests that diurnal regulation may play a dominant role in processes such as myelin turnover or sheath maintenance, at least during the developmental period tested (approximately P28 mice). In contrast, neuronal populations such as CA1 exhibit fewer circadian DEGs overall, with enrichment of lncRNAs at ZT12 compared to other cell types, hinting at separate mechanisms of regulatory control. Together, these results suggest that circadian transcription exerts heterogeneous but highly structured influence across hippocampal cell types.

Circadian cycle has limited effect on activity-dependent transcription

We found no circadian differences in IEG induction (Fig. 6A–D), despite expecting higher expression in the dark period due to increased locomotor activity⁵², state-dependent hippocampal network events such as sharp-wave ripples during rest¹⁴, and prior reports of circadian oscillations in cAMP, MAPK, and CREB signaling^50,52,64,65. Instead, IEGs were induced uniformly across the diurnal cycle in both neuronal and glial populations, exemplified by CA1 cells (Fig. 6). Prior work reported dark-phase upregulation of IEGs in rats under chronic EE⁶⁶ or in cortex at ZT20⁵¹, but those studies differ in species, brain region, EE paradigm, stimulus, timing, and bulk vs. snRNA-seq, making direct comparison to our results inappropriate.

It has been previously shown that Per genes are inducible in the hippocampus in response to circuit activity^18,19,48, a finding that we recapitulate after 30 minutes of EE. This finding, in conjunction with Per’s circadian oscillations in the absence of activity, raises important conceptual questions. What does it mean for there to be circadian-dependent expression when these genes are also induced by activity? We observe that at the highest peaks of EE-responsive induction, Per expression is comparable to its circadian peak. One speculation that could explain this phenomenon is that in the absence of activity, daily peaks in Per induction “guarantee” a dose of permissive regulatory-element access for PER-CRY target genes which would otherwise not be transcribed but are nevertheless necessary for cell function; activity-dependent Per expression would act as an adjuvant to predictable circadian bursts, possibly at genes important for cellular adaptations to circuit activity. Clearly, additional studies are necessary to understand the multifactorial causes and functional consequences of Per expression.

When testing for ZT:activity interaction genes, only 13 DEGs were found, all in neuronal subclasses. The functional classes of these DEGs suggest that there is a special role for fluctuations in intercellular endocannabinoid signaling throughout the circadian cycle, an enticing idea that may be relevant for neurological phenotypes. For example, despite the untold number of molecular etiologies for different epilepsies, circadian rhythms to seizure attacks are common in drug-resistant epilepsy, occurring in approximately 90% of patients²⁵. This observation implies that there are circadian changes taking place in the brain that are not themselves pathogenic, but may affect seizure threshold across the diurnal cycle. In summary, our examination of the intersection of activity- and circadian-dependent gene expression is an important first step toward understanding the circadian transcriptome in the hippocampus, as it may canalize certain disease phenotypes.

Materials and Methods

Animals

Animal experiments were approved by the University of California San Diego Institutional Animal Care and Use Committee and followed guidelines according to the National Institutes of Health’s Guide for the Care and Use of Laboratory Animals. All experiments used wild-type C57BL6/J mice from The Jackson Laboratory (JAX#000664). Mice were approximately four weeks old, with balanced proportions of males and females distributed across control and experimental groups (Supp. Table 1).

Experimental procedures

Standard environment and circadian timepoints

Mice were group-housed in a vivarium on a 12-hour light/dark cycle and provided food ad libitum. To avoid non-experimental stimulation artifacts resulting from cage handling or transport, mice were transported from the vivarium to the standard environment (SE) at least 18 hours before experiments. The SE consisted of a plywood box paneled with sound-dampening foam, air-circulation fans, and LED lights whose power was controlled by a timer mimicking the vivarium’s light/dark cycle. Circadian timepoints were measured and reported on the 24-hour Zeitgeber time scale, with ZT0 marked as “lights-on” (at 06:00 PST) and ZT12 as “lights-off” (at 18:00 PST).

Enriched environment

For enriched-environment (EE) experimental conditions, animals were transferred from their home cage in the standard environment to a larger cage in a different room containing plastic and metal toys, tunnels, vertical scaffolds, and a running wheel. Mice were placed in the EE with at least one cagemate and left to explore continuously for the durations reported (i.e., for 30 minutes or 6 hours).

Kainic-acid seizures

Pharmacological seizures were induced using peritoneal injection of 15 mg/kg of kainic acid (KA; Santa Cruz, #SC-200454) as described previously^37,67. Seizure onset was measured as the onset of rearing and accompanying bilateral forelimb clonus, approximately corresponding to Racine stage 4⁶⁸. Mice were sacrificed either 30 minutes or 6 hours after seizure onset.

snRNA-seq sample collection

For sacrifice and sample collection, mice were deeply and rapidly anesthetized with isoflurane in an induction chamber. Following induction, the brain was removed from the skull into a dissection dish containing ice-cold dPBS (Gibco Cat. No. 14040133). Hippocampus was rapidly extracted from one hemisphere, and the approximate CA1 area was microdissected using a microblade. Although the purpose of this microdissection was to de-enrich samples of dentate granule cells, we still note that these cells nevertheless comprise the 4th largest subclass in our dataset (Fig. 1B). The dissected tissue sample was then placed in a microtube and flash-frozen. The other hemisphere was placed in a cryomold as a whole, embedded in OCT compound, and rapidly flash-frozen in isopentane chilled with dry ice. Left and right hemispheres were alternated for the collection of CA1 and whole-hemisphere collection, across sexes and conditions, to avoid any effects of lateralization (Supp. Table 1). Both samples from each animal were stored at -70°C until use.

Nucleus isolation and fixation

Flash-frozen hippocampal dissectates were stored at -70°C until all experimental samples had been collected. Before nucleus isolation, all buffers, tubes, and instruments were pre-cooled on ice for 15-20 minutes before tissue was removed from -70°C. Reusable instruments such as the Dounce tubes and pestles (Millipore cat. #D8938-1SET) were rinsed with RNAse-free PBS and then with RNAse Away detergent before being chilled on ice. All pipette tips (filtered tips only) and disposable tubes contacting samples were RNAse-free. All buffers used during nucleus isolation contained recombinant RNAsin to preserve RNA integrity (Promega, cat. #N2111), and actinomycin D to prevent artifactual transcription stimulated by the isolation protocol (Tocris, cat. #1229). Two main buffers were used for nucleus isolation: nucleus-isolation medium (NIM) and NIM-Douce-pestle buffer. NIM buffer consisted of 250 mM sucrose, 25 mM KCl, 5 mM MgCl2, 10 mM Tris-HCl (pH 7.4), 1 mM DTT, 0.2 U/uL RNAsin, 25 uM actinomycin D, and 1 cOmplete protease inhibitor mini tablet (Roche cat. #11836170001), dissolved in 10 mL of buffer. NIM-DP buffer consisted of 0.1% triton-X 100 detergent in NIM buffer.

To begin nucleus isolation, tissue samples were removed from -70°C and placed directly on ice to thaw. No more than five samples underwent isolation at any one time in order to minimize sample-handling time before fixation. First, 1 mL of NIM-DP buffer was added to the sample storage microtube, and the sample was transferred to a 2 mL Dounce tube using a 1000 uL pipette tip. Samples were homogenized using a Douce pestle in 1 mL of NIM-DP buffer, first with 10 strokes of the loose pestle, and then 17 strokes of the tight pestle. Two sets of dounce tubes and pestles were used at a time, and each was washed five times with PBS in-between samples. Dounced tissue homogenates were then spun in a fixed-angle benchtop centrifuge at 800 x g for 8 minutes at 4°C. The supernatant was removed, and the nuclear pellet, clearly visible, was resuspended thoroughly in 500 uL of NIM buffer using 30 strokes with a P200 pipette tip. Suspensions were then filtered through a 70 um filter before being spun at 200 x g for 10 min at 4°C in a swinging-bucket centrifuge. The supernatant was removed, and nuclei pellets were resuspended in 187.5 uL of Parse Biosciences’s “Resuspension Buffer” and strained through a 20 um filter to begin Parse’s Evercode Nuclei Fixation protocol (version 1.1, kit chemistry v3). In total, samples/nuclei were suspended in actinomycin D and kept on ice or in the pre-chilled centrifuge for the entirety of the isolation protocol, which took approximately 1 hour before proceeding directly to the Parse fixation protocol.

snRNA-seq barcoding and library construction

snRNA-seq library construction was performed using Parse Biosciences’s Whole Transcriptome v3 kit according to the manufacturer’s protocol. Briefly, nuclei originating from different samples are multiplexed and barcoded according to the SPLiT-pool combinatorial barcoding method²⁸, wherein individual nuclei go through serial rounds of PCR barcoding, pooling, and splitting amongst a 96-well plate containing well-specific oligo barcodes. Split-pooling occurs in four serial barcoding rounds to ensure that each nucleus has a unique combination of barcodes ligated. All samples were loaded and barcoded together to avoid batch effects from separate library preparations. Libraries were sequenced at the UCSD Institute for Genomic Medicine on an Illumina NovaSeq X Plus using paired-end 100-bp reads. Sequencing was targeted at approximately 90,000 reads per nucleus, with a final average read depth of 66,702 reads per nucleus.

Quantification and statistical analysis

Data preprocessing and quality control

Nuclei quality control was performed as outlined in Supplementary Figure 1. Raw nucleus x gene count matrices from demultiplexed samples were generated from the FASTQ files using Parse Bioscience’s Trailmaker software pipeline, version 1.4.1. A “raw” Seurat object created by Trailmaker containing 110,728 nuclei, along with the digital gene-count matrix, associated sample metadata, and DoubletFinder scores⁶⁹ was downloaded and used to begin nuclei QC. Percent of reads mapping to mitochondrial reads was calculated for each nucleus. Nuclei that did not meet the following criteria were removed: transcript count (UMI) >= 1500; gene count >= 1500; and percent-mitochondrial-read count <= 5%. After this filtering step, 108,748 nuclei remained. Allen Brain Cell (ABC) taxonomic labels were assigned using the MapMyCells online tool. Next, single-nucleus raw gene counts were transformed using Seurat’s SCTransform function (v0.4.1), with the vars.to.regress argument set as the percent of reads mapping to mitochondrial genes (percent.mt) and SPLiT-Seq sublibrary. SCTransformation was computed across nuclei from all samples. (Note: Across-sample SCT-normalization, rather than within-sample normalization followed by integration was determined appropriate after normalizing, clustering, and visualizing UMAPs on a pilot sequencing run from this dataset, where we observed no batch effects due to sample or SPLiT-seq sublibrary). SCT-normalized counts were then used to perform principal component analysis (PCA), nearest-neighbor graph construction, initial Leiden cluster assignment, and uniform manifold approximation and projection (UMAP) embedding. This resulted in the segmentation of 108,748 nuclei into 47 Leiden clusters.

At this step, manual curation of 47 detected Leiden clusters was performed by experimenter J.O. cross-referencing plots and tables of UMI count, gene count, percent mitochondrial reads, doublet scores, and putative ABC cell type within each Leiden cluster. Based on these judgements, nuclei in 6 small Leiden clusters were removed because of low counts, high mitochondrial reads, high doublet scores, or the collision of many different putative ABC cell types. Finally, any nuclei with a DoubletFinder score < 0.02 were removed. After these quality-control steps, 93,066 nuclei remained. These again went through final PCA, Leiden clustering, and UMAP embedding to create the UMAP in Figure 1.

ABC taxon assignment stability

Stability of ABC taxonomy assignment for single nuclei was performed by comparing two separate runs of the curated data through the ABC atlas MapMyCells tool. The first run was performed using the native mouse transcriptome. The second run was performed after removing ARGs that had been ascertained in a published study of activity-dependent gene expression². The stability of assignment was then determined by comparing taxon equality within the same barcode and visualized with Sankey flow diagrams (Fig. 1E-H).

Log2(FC) methods

Average log2(FC) expression values in Figures 3 and 6 represent pseudobulk fold-change expression values calculated using DESeq2.

Differential gene expression

To identify differentially expressed genes (DEGs) across activity conditions, we performed pseudobulk DGE analysis independently for each Subclass using DESeq2⁷⁰. For each Subclass, nuclei were subset and grouped by their sample of origin. Pseudobulk counts were aggregated at the sample level using Seurat’s AggregateExpression() function. Genes expressed in fewer than 1% of nuclei were excluded. DESeq2 was then used to construct a DESeqDataSet with sample-level counts and a design formula including activity condition (expression ∼ activity_condition). Only condition contrasts meeting a minimum per-group nucleus count threshold (n = 30 nuclei per condition) were retained for testing. Wald tests were performed, and log2 fold changes were shrunken using the apeglm method. Genes with absolute log2 fold change > 0.585 (corresponding to 1.5-fold change) and fdr < 0.05 were classified as differentially expressed. For analysis discovering DEGs between ZTs, log2 fold change values were not shrunken. This was done for ZT analyses but not activity analyses to minimize type II error rate; there are fewer biological replicates at each ZT (Supplementary Table 1). Results were saved separately for each Subclass and contrast.

Hierarchical clustering

For heatmap plots accompanied by a dendrogram, genes were ordered according to hierarchical clustering of the Euclidean gene distances using the package ComplexHeatmap⁷¹. In the heatmaps of Figure 5, genes were first clustered into 4 k-means clusters, and only after were organized according to hierarchical methods within each k-means cluster.

ZT:Activity interaction testing

To test for statistical interactions between the effect of ZT and the effect of activity on gene expression, we used the linear formula

to test for a significant ZT:Activity interaction term with DESeq2. Pseudobulk gene counts were computed for each subclass separately. log2 fold changes were not shrunken in this analysis in order to increase sensitivity of this test. For a priori gene sets, multiple-hypothesis correction was performed by calculating the FDR among the a priori gene set. For example, the FDR rates reported in Figure 6A-D were corrected for 15 multiple comparisons for the set of 15 IEGs. Only genes with an FDR < 0.05 were labeled as significant DEGs.

ERG vs LRG classification

We used a similar method to a prior single-cell study to classify DEGs as mutually exclusive early-response genes (ERGs) or late-response genes (LRGs)¹⁰. Classification was performed separately for excitatory neurons, inhibitory neurons, and glial Subclasses. First, DEGs for all Subclasses and contrasts (EE30m vs SE, EE6h vs SE) were aggregated. Genes were provisionally labeled as ERG, LRG, or both based on differentially expressed in EE 30 min or EE 6 hr. Genes initially labeled as “both” were reclassified based on their peak fold-change timepoint. In cases where gene classification differed across Subclasses, the modal classification was assigned. Final classifications were saved for downstream analysis (Figure 3).

Tau calculation

We calculated a tau score for each gene, adapted from ref. ³⁶, to discover and quantify highly cell-type specific gene expression. For example, a gene with a tau score of 0 would have promiscuous expression among the cell types examined, whereas a gene with a tau score approaching or greater than 1 would have highly cell-type specific expression. For analyses on non-fold-change expression data among different cell types (e.g. Figure 5), we directly applied the tau formula, in which tau was bounded by (0, 1):

where n is the number of subclasses, and expression_i represents the expression of the gene in subclass i; in these analyses, our tau threshold was set at 0.75 due to the index being bounded by (0, 1). For analyses on fold-change expression where expression values can be negative, the tau score was bounded by (0, ∞); for these analyses, our tau threshold was set higher at 0.8. In practice, tau scores above 0.75 were observed to be highly specific to one or a small number of subclasses.

Defining active cells

In order to determine the number of cells activated within a given condition, it was determined whether cells expressed the IEGs above a certain threshold: The IEG list was defined a priori based on the following 15 canonical IEGs: Arc, Btg2, Dusp1, Dusp5, Egr1, Egr2, Egr3, Egr4, Fos, Fosb, Fosl2, Junb, Npas4, Nr4a1, and Nr4a3. The induction threshold was determined for each gene individually, and was defined as the 90th percentile of expression in SE cells of the CA1 subclass. Thus, for each gene in the gene list, a binary classification was assigned as either induced or not-induced based on whether or not it surpassed the expression threshold (Fig. 4A). Finally, a cell was defined as “active” if it induced at least 20% of the genes (i.e., 3 or more) in the activity-regulated-gene list.

Data availability

We are not able to upload data files to the GEO database at time of submission due to the US Federal Government shutdown, which has limited access to the GEO FTP server. Fastq files for sequenced reads and final Seurat objects will be publicly deposited in the GEO database. Metadata for biological replicates can be downloaded from the GEO database or https://github.com/Bloodgood-Lab/Olmstead-2025-preprint. We used the Trailmaker pipeline running split-pipe v1.4.1 through Parse Biosciences' web portal for creation of the cell-count matrix and unprocessed Seurat object, which will be publicly deposited in the GEO database. All the code and package versions to reproduce this manuscript are publicly available at https://github.com/Bloodgood-Lab/Olmstead-2025-preprint.

Acknowledgements

We thank members of the Bloodgood lab and Dhanajay Bambah-Mukku for continuous discussion and feedback about experimental results and data analysis. We are indebted to Eric Halgren for providing computing resources. We thank Elizabeth Chamiec-Case and Lauren Valdez for their assistance with experimental troubleshooting. We kindly acknowledge the detailed experimental and analytic support provided by Brian Davis, Sam You, and Daniel Diaz of Parse Biosciences. Initial analysis of these data was published in a dissertation by J.A.O.²⁷ This work was supported by funding awarded to BLB by the National Institute of Neurological Disorders and Stroke (R01 NS111162). This publication includes data generated at the UC San Diego IGM Genomics Center utilizing an Illumina NovaSeq X Plus that was purchased with funding from a National Institutes of Health SIG grant (#S10 OD026929).

The authors declare no competing interests.

Additional information

Author contributions

• Jack A Olmstead: Conceptualization, Investigation, Methodology, Data curation, Formal analysis, Software, Visualization, Writing – original draft, Writing – review and editing

• Lauren E King: Investigation, Writing – review and editing

• Brenda L Bloodgood: Conceptualization, Supervision, Project administration, Funding acquisition, Writing – original draft, Writing – review and editing

Funding

HHS | NIH | National Institute of Neurological Disorders and Stroke (NINDS) (R0 NS111162)

• Brenda L Bloodgood

HHS | NIH | National Institute of General Medical Sciences (NIGMS) (T32 GM154642)

• Jack A Olmstead