snRNA-seq Analysis of the Mouse Hippocampus Reveals Changes in Cell-Type Abundance Upon Exercise

To understand how physical exercise affects the hippocampus, we designed an experiment where 8 wild-type (WT) mice were divided into 2 equally sized groups—runners that had free and voluntary access to running wheels in their cages (the “exercise” group), and sedentary mice that had similar wheels in their cages that were however blocked (the “control” group) (Supplementary Table 1). This voluntary exercise paradigm lasted 4 weeks, after which we isolated whole hippocampal nuclei from all mice for snRNA-seq analysis (Fig. 1A). A total of 22817 nuclei across the 8 samples were recovered after quality control, doublet removal, and filtering (see methods). After accounting for batch effects, clustering was performed and nuclei were grouped into 21 unique clusters. Subsequently, these clusters were assigned to major hippocampal cell types including excitatory neurons (ExN), inhibitory neurons (InN) as well as microglia, astrocytes, oligodendrocytes, and oligodendrocyte precursor cells based on the expression of cell-type specific marker genes (Fig. 1B, C; Supplementary Fig. 1 and Supplementary Table 2). Excitatory neurons made up the majority of cells, followed by oligodendrocytes and inhibitory neurons (Fig. 1D). Additionally, we also identified excitatory neuron clusters that belonged to either the CA regions (ExN2-10 and ExN12-13) or the dentate gyrus (ExN1 and ExN11) within the hippocampus, using marker genes for these regions [62] (Fig. 1C, E and Supplementary Table 2). Additionally, distinct marker genes for each of these clusters were determined computationally (see methods) (Supplementary Table 3), and gene ontology (GO) analyses of the top markers specific to each of the neuronal clusters suggested different biological functions for these neurons within the hippocampus (Supplementary Fig. 2 and Supplementary Table 4–6).

Fig. 1

Single-nucleus RNA sequencing analysis of the whole hippocampus in exercising vs. sedentary mice reveals changes in the abundance of specific cell-types. A Experimental design: The cohort of WT mice belonging to the exercising experimental group, or “runners” (n = 4), had free access to running wheels in their cages for a duration of 4 weeks. Mice belonging to the control or “sedentary” group (n = 4) were similarly housed but the running wheels were blocked. The 4-week-long voluntary exercise paradigm was followed by the isolation of hippocampal nuclei for single-nucleus RNA sequencing. B UMAP plots showing clusters of nuclei or “cells,” colored by experimental group (left panel) and cell-type (right panel). (ExN1-13, excitatory neurons; InN1-4, inhibitory neurons; ODC, oligodendrocytes; OPC, oligodendrocyte precursor cells; AST, astrocytes; MGL, microglia) C Violin plots showing average module score (expression) of marker genes specific to the different cell types/hippocampal regions, after cell-type specific annotation of clusters (ExN, excitatory neurons; InN, inhibitory neurons; ODC, oligodendrocytes; OPC, oligodendrocyte precursor cells; AST, astrocytes; MGL, microglia; CA, Cornu Ammonis; DG, Dentate Gyrus). D Bar graph indicating the proportions of broad cell-types observed in the dataset (ExN, excitatory neurons; InN, inhibitory neurons; ODC, oligodendrocytes; OPC, oligodendrocyte precursor cells; AST, astrocytes; MGL, microglia). E UMAP plot with cell clusters colored by average module score (expression) of marker genes specific to the dentate gyrus (DG), highlighting the two DG excitatory neuron clusters. F Analysis of differences in cell-type proportions between cells from exercise and control samples, using permutation testing. The x-axis denotes the fold difference in cell-type proportions (log2 scale). Points marked in red indicate clusters with significantly different proportions of cells between the two groups. Horizontal lines around the points indicate the confidence interval for the magnitude of difference for a specific cluster, calculated via bootstrapping. G UMAP plot with cells colored by Augur cell-type prioritization upon perturbation resulting from exercise, measured using the area under the curve (AUC) scores 

Next, we applied permutation testing to analyze if exercise would lead to detectable differences in the proportions of cells originating in each cluster. Out of all clusters, InN2, ExN10, and ExN11 showed significant changes in cell-type abundance (Fig. 1F). Using Augur [52], we also performed cell-type prioritization analysis to detect which cell types are most perturbed in response to exercise. To measure the perturbation levels, we employed the area under the receiver operating characteristic curve (AUC). An AUC value of 0.5 indicates that cells from the exercise condition in a cluster have no significant difference in perturbation compared to cells from the control condition (random chance). Conversely, a value of 1.0 signifies that every cell from the exercise condition exhibits higher perturbation compared to the control condition. Augur performs subsampling in a way that changes in relative abundances of cell-types across conditions do not confound the analysis. Hence, if certain cell-types are found to be strongly perturbed due to exercise, it could imply that differences in cell-type proportions appear due to these perturbations and not the other way around. We observed that clusters InN2 and ExN11 were the only two clusters with an AUC value > 0.6, indicating changes in cell-type proportions within these clusters upon exercising (Fig. 1G). More specifically, upon exercise, fewer cells were detected for cluster InN2 while more cells were detected within cluster ExN11. In particular, the increased proportion of cells within cluster ExN11 may reflect increased adult neurogenesis, which has been repeatedly described upon exercise (Kempermann, 1997) (Van Praag, 1999) [63]. Thus, while the changes observed for cluster ExN11 may indicate the increased number of excitatory neurons upon neurogenesis, the changes seen within cluster InN2 are more difficult to explain. Therefore, we decided to first analyze cluster InN2 in greater detail.

Exercise Induces Changes in an Inhibitory Neuron Population Expressing Prdm16

We identified four distinct inhibitory neuron clusters in the dataset based on the overall expression of marker genes Gad1 and Gad2. Although we could not attribute the expression of canonical inhibitory neuron subtype markers to these clusters, each of these 4 clusters was instead characterized by high specific expression of the following genes: Sox6 (InN1), Prdm16 (InN2), Cnr1 (InN3), and Egfr (InN4) (Supplementary Fig. 3). As indicated in Fig. 1E, we observed a selective reduction in abundance (log2 fold difference =  − 2.69) of inhibitory neurons in cluster InN2 in the exercise condition, as compared to the controls (Fig. 2A). Further functional analysis of the top marker genes for InN2 cells distinguishing them from the other inhibitory clusters (i.e. differentially expressed (upregulated) in InN2 as compared to InN1, InN3, and InN4) broadly implicated genes involved in transmembrane transport of ions as well as dendritic spine development (Fig. 2B and Supplementary Table 7). Among the top 20 of these markers, a few were more distinctly expressed in InN2 as compared to the rest. These included Prdm16, Ano1, Ano2, Zfhx3, Zic1, Zic4, Zfp521, and Ankfn1 (Supplementary Fig. 4A). Prdm16 is a transcription factor that belongs to the PRDM family of transcriptional regulators. It has been widely characterized as an important cell-fate switch regulator in brown adipose tissue [64, 65], but recent studies have also identified its role in neural stem cell homeostasis (including maintenance of the neural stem cell pool), as well as the development, differentiation, and positional specification of neurons [66,67,68,69,70]. Ano1 and Ano2 are calcium-activated chloride channels (CaCCs) of the anoctamin protein family. CaCCs play crucial roles in regulating the excitability of smooth muscle cells and also some types of neurons [71]. Specifically, Ano1 has been shown to contribute to the process maturation of radial glial cells during cortex development [72], while Ano2 modulates action potential waveforms in hippocampal neurons, along with the integration of excitatory synaptic potentials [73]. Moreover, Ano2 has been implicated in the calcium-dependent regulation of synaptic weight in GABAergic inhibition in the cerebellum, hence regulating ionic plasticity [74]. Zfhx3 is a large transcription factor with 23 zinc fingers and 4 homeodomains [75], which has been shown to promote neuronal differentiation during neurogenesis in development and primary cultures [76, 77]. Zic1 and Zic4 are also transcription factors containing zinc finger domains, which are involved in nervous system development [78, 79], specifically by regulating neuron differentiation and maintaining neural precursor cells in an undifferentiated state [79]. Zfp521 encodes another zinc-finger protein that regulates many genes involved in neural differentiation. Studies have shown that Zfp521 is essential and sufficient for driving the intrinsic neural differentiation of mouse embryonic stem cells [80] and is also sufficient for converting adult mouse brain-derived astrocytes into induced neural stem cells [81]. Little is known about Ankfn1 but it has been genetically linked to cannabis dependence [82].

Fig. 2

Selective loss of InN2 inhibitory neuron cluster upon exercise indicates a role of the transcription factor Prdm16. A UMAP plot with cells colored by the experimental group, showing the loss of InN2 inhibitory neurons upon exercise. B Heatmap plots of functional annotation for specific genes among the top 50 markers of the InN2 cluster that had significantly enriched GO (biological process, molecular function, and cellular component) terms. C Cell-type specific regulons for InN2 cluster identified using the SCENIC workflow. The y-axis denotes the regulon specificity score (RSS) (with high RSS values indicating high cell-type/cluster specificity, and vice versa). The x-axis denotes the rank of each regulon within the selected cluster, based on the RSS. The top 5 ranked regulons for InN2 are labeled on the plot, with the number of genes comprising each regulon indicated within parentheses. (Regulons ending with “extended” also include motifs linked to the transcription factor by lower confidence annotations). D Dot plot showing significant GO biological process terms enriched among the collective list of genes and transcription factors (TFs) making up the top 5 regulons, as indicated in the RSS-Rank plot in (C). E Violin plots depicting the normalized expression of Prdm16 (top panel), and the average module score (expression) for the genes included in the Prdm16 regulon (bottom panel) in all clusters. F Violin plot showing the expression of Prdm16 in the InN2 cluster, split between exercise and control cells. G Network plot showing the 17 genes comprising the Prdm16 regulon. H Dot plot showing significant GO biological process terms enriched among the list of 17 Prdm16 regulon genes

To further understand what makes the InN2 inhibitory neurons distinct, we used the SCENIC workflow [53] to identify regulons (sets of transcription factors and their putative targets) that were specifically active in the InN2 cluster (Fig. 2C and Supplementary Table 8). Enrichment analysis of the collective list of genes comprising the top 5 InN2-specific regulons (Prdm16, Rfx2, Gtf2ird1, Dlx1 and Zfp941) indicated significant enrichment of gene ontology (GO) terms for regulation of synaptic signaling and transmembrane ion transport, membrane depolarization and response to acetylcholine, and interestingly, regulation of learning and memory (Fig. 2D and Supplementary Table 9). Among these top 5 regulons, the transcription factor Prdm16 and its 17 target genes comprised the most specific active regulon in this cluster. This was substantiated by the highly specific expression in InN2 of both Prdm16 and the combined expression of the 17 regulon genes (Fig. 2E). The observed reduction of exercise cells in InN2 was subsequently indicated in the stark upregulation of Prdm16 in the control cells (Fig. 2F). Five candidates from the Prdm16 regulon (Prdm16, Ankfn1, Meis2, Myo5b, Zfhx3, Zfp521) are also shared with the top 20 InN2 marker genes (differentially expressed in InN2 as compared to InN1, InN3, and InN4) (Supplementary Fig. 4), most of which have been implicated in regulating neuronal differentiation (Liu) [69, 70, 76, 77, 80, 83]. Additionally, functional analysis specifically on the Prdm16 regulon genes suggests enrichment of GO terms such as rhythmic process (Ankfn1, Hs3st2, Rorb, Zfhx3) and regulation of bone morphogenetic protein (BMP) signaling (Prdm16, Gpc3, Zfp423), which is known to modulate neuronal maturation and adult neurogenesis [84] (Fig. 2G, H and Supplementary Table 10).

Among the other top regulons, the gene candidates within the Dlx1 regulon indicated that Dlx1 itself positively regulates Prdm16. Dlx1 is a TF that, along with its related gene Dlx2, is involved in regulating neuronal migration (of newborn neurons, away from the proliferative zone) as well as the functional longevity of GABAergic interneurons in the hippocampus [85]. Moreover, a recent study indicated that the expression of the Meis2 gene driven by Dlx1/2 promoted the fate determination of striatal neurons in mice [86]. Since Meis2 is regulated by Prdm16 (as seen in the Prdm16 regulon), and it is also one of the top markers of the InN2 cluster which is comprised mostly of control cells, these findings collectively suggest that exercise leads to repression of Prdm16 in specific hippocampal neurons (which are characterized by high Prdm16 expression). Although we identified these neurons initially as cluster InN2, it is possible that these cells represent a specific developmental stage during adult neurogenesis that is only captured in control mice in our dataset, since exercise is known to increase proliferation but also maturation of adult newborn neurons [87,88,89].

Exercise Affects a Subpopulation of Excitatory Neurons Linked to Enhanced Adult Neurogenesis

Following the cell-type abundance and perturbation analyses (Fig. 1E and F), we also sought to take a deeper look into the ExN11 cluster belonging to the dentate gyrus. These cells showed an increase in abundance (log2 fold difference = 0.79) in the exercise condition, as compared to the controls. This increase can also be seen in a small fraction of cells from ExN11 that seem to integrate into the mature granule cell cluster ExN1 of the dentate gyrus (as indicated in Fig. 3D by the expression of GC markers Prox1 [62, 90, 91] and Calb1 [92], which shows an overall low but relatively higher expression in ExN1), with the number of these integrating cells increasing more than twofold after exercise (Fig. 3A).

Fig. 3

Increased abundance of neurons in the ExN11 cluster suggests increased neurogenesis upon exercise. A (left panel) UMAP plots highlighting the ExN11 cluster, split by the experimental conditions (exercise and control). B Dotplot showing significant GO biological process terms enriched among the top 50 genes differentially expressed between ExN11 and ExN1 clusters, and among the top 50 genes differentially expressed between the ExN11 cluster and the ODC and OPC clusters. C Partition-based graph abstraction (PAGA) graph for InN2, ExN11, and ExN1 clusters, with velocity-directed edges constructed from RNA velocity measurements. Edges denote either connectivities (dashed) or transitions (solid/arrows). D Box-plots showing a significant decrease in the proportion of developing cells (ExN11/InN2) and an increase in the proportion of mature granule cells (ExN1) upon exercise. Each dot indicates one sample (n = 4/4 for exercise and control samples). E Violin plots showing normalized expression of selected marker genes for radial glia/immature neurons/exercise-mediated neurogenesis. F Cell-type specific regulons for ExN11 cluster identified using the SCENIC workflow. The y-axis denotes the regulon specificity score (RSS) (with high RSS values indicating high cell-type/cluster specificity, and vice versa). The x-axis denotes the rank of each regulon within the selected cluster, based on the RSS. The top 5 ranked regulons for ExN11 are labeled on the plot, with the number of genes comprising each regulon indicated within parentheses. G Dot plot showing significant GO biological process terms enriched among the collective list of genes and transcription factors (TFs) making up the top 5 regulons, as indicated in the RSS-Rank plot above. H Combined gene-regulatory network plot with TFs specific to InN2 and ExN11 clusters, and their respective regulon genes (larger nodes represent genes/TFs that connect two or more regulon networks)

In order to understand the underlying gene expression patterns of ExN11 cells, we looked at the top computationally detected marker genes for this cluster. Most of these markers were not specific to ExN11, and in fact, were also highly expressed in either ExN1 neurons or oligodendrocyte lineage cells (ODC and OPC clusters). This becomes more evident after plotting the combined expression of the top 50 genes differentially expressed between ExN11 and ExN1 (Supplementary Fig. 5A and Supplementary Table 11), and the top 50 genes differentially expressed between ExN11 and the oligodendrocyte lineage clusters (ODC and OPC) (Supplementary Fig. 5B and Supplementary Table 12). The functional enrichment analysis of these two lists of differentially expressed markers indicates enrichment of GO terms related to the ensheathment of neurons, regulation of synaptic plasticity, neurogenesis, and neuronal migration (Fig. 3B and Supplementary Table 13), suggesting that increased neurogenesis upon exercise results in more excitatory neurons in ExN11.

Next, we generated a partition-based graph abstraction (PAGA) graph with velocity-directed edges constructed from RNA velocity measurements [60] [61], after subsetting the dataset to analyze our clusters of interest: InN2, ExN11, and ExN1 (Fig. 3C). The graph indicates a potential connectivity between these three clusters, and a transition from ExN11 to ExN1, which supports our findings so far suggesting that ExN11 represents newborn excitatory neurons that eventually integrate into the dentate gyrus. Moreover, we see that the proportion of cells coming from ExN11 and InN2 combined, taken among the total cells from the three clusters of interest (ExN1, ExN11, and InN2), decreases in exercise, while the proportion of mature granule cells (ExN1) increases in the exercise condition as compared to controls (Fig. 3D). This suggests that exercise leads to faster maturation of developing or immature neurons (InN2 and ExN11) to mature granule cells (ExN1).

We also looked at the expression of markers corresponding to different stages of neurogenesis in this cluster, in order to assign these cells to a specific stage of adult neurogenesis [93] (Fig. 3E and Supplementary Fig.5C). However, we observed either no expression (Neurod, Tbr2, Tubb3) or very low expression (Gfap, Nestin, Pax6) of many canonical markers, possibly due to dropout effects in the sequencing data. As compared to other neurons, ExN11 cells show higher expression of GLAST (Slc1a3), which is a marker for radial glia cells, but at the same time, we also observed low expression of Dcx and TUC-4 (Dpysl3) and high expression of Prox1 which could suggest that these are transiently amplifying (or transit amplifying) progenitor cells [94] (Torii, 1999). This latter cell-type is known to show transient expression of markers: initially expressing glial or stem cell markers like Gfap and Nestin/Sox2, and later expressing granule cell markers such as Prox1. ExN11 cells also express stathmin (Stmn1) at relatively higher levels. Stathmin is a microtubule destabilizing protein that is considered to be an immature neuron marker since it has been implicated in controlling the transition from dividing neuronal precursors to postmitotic neurons in the subgranular zone of the DG during adult neurogenesis [95]. Bdnf (brain-derived neurotrophic factor), which is known to be involved in exercise-mediated neurogenesis in the hippocampus [96,97,98], showed very low expression throughout the dataset, which made it difficult to reliably quantify any overall changes in its expression levels in the exercise condition as compared to the sedentary controls. However, it did have relatively higher expression levels in ExN11, although we found no significant differences between the exercise and control conditions in this cluster.

We further analyzed the gene regulatory landscape for ExN11 using the SCENIC workflow to identify regulons that are specifically active in this cluster and found the Meis1 regulon to be the most specific out of the top 5 regulons (Fig. 3F and Supplementary Table 14). Meis1 is a transcription factor that has been shown to play a role in promoting neuronal differentiation and possibly also neurogenesis [99,100,101]. Other top transcription factors include Cebpd (CCAAT enhancer binding protein delta) which has been implicated in hippocampal neurogenesis and regulation of learning and memory [102, 103], Jun which is an immediate early gene and transcriptional regulator involved in cell proliferation [104], Zfhx2 which is involved in neuronal differentiation [105], and Foxo1 which regulates the long-term maintenance of adult neural stem cells and belongs to the FoxO family that plays an important role in the maintenance of autophagic flux and neuronal morphogenesis in adult neurogenesis [106, 107]. A functional enrichment analysis of the collective list of genes comprising these top 5 ExN11-specific regulons indicated significant enrichment of gene ontology (GO) terms for gas transport (Car2, Car14—carbonic anhydrases that play a role in neuronal excitability and signaling [108]), regulation of neural development, neurogenesis, and also glial cell differentiation, which is in line with our previous observation of glial markers in this cluster (Fig. 3G and Supplementary Table 15). We also observed that some of the TFs from the top regulons of ExN11 and InN2 clusters shared downstream target genes and also regulated each other, such as Rfx2 regulating Jun as well as Prdm16. These common gene-regulatory relationships between our two clusters of interest are highlighted in Fig. 3H (to explore the network in high magnification please view Supplementary Fig. 6). Building on our previous hypothesis, these results further suggest that the observed loss of Prdm16-expressing neurons upon exercise might represent a molecular snap shot of enhanced adult neurogenesis and especially the maturation of new born neurons. In this scenario, ExN11 neurons could be a consequence of their differentiation into and possible switch to an excitatory neuron subtype.

Differential Gene Expression and Gene Regulatory Network Analyses Reveal Broad Transcriptional Changes Linked to Synaptic Plasticity

To get a better understanding of the broad transcriptomic changes mediated by exercise, we looked at differential gene expression between cells from the exercise and control conditions in all other clusters apart from InN2 and ExN12 (a complete list of differentially expressed genes is given in Supplementary Table 16). Most clusters revealed either no significant changes in gene expression or very few deregulated genes across the two conditions, with the exception of ExN1, ExN2, ExN6, InN1, and ODC (Fig. 4A (top) and B (top)). The top significant GO (biological process) terms enriched for these deregulated genes indicate broad processes involved in specific clusters upon exercise (Fig. 4A (bottom) and B (bottom), Supplementary Table 17–18). In particular, ExN6 shows a distinctly large number of both up- and downregulated genes, with the former being enriched for dendrite morphogenesis and synaptic signaling pathways and the latter for postsynapse organization and synapse assembly processes (Supplementary Table 19–20). UpSet plots in Fig. 4A and B identified certain genes that are commonly up- or downregulated upon exercise across multiple clusters. Vps13a, for example, is upregulated in exercise in 8 out of the 21 clusters and also downregulated in 2 excitatory neuron clusters, suggesting a role of chorein (the protein encoded by Vps13a) which is a powerful regulator of cytoskeletal architecture and cell survival in many cell types [109]. Moreover, Vps13a mutations lead to chorea-acanthocytosis, a rare disorder that also affects the brain and reduced chorein levels have been linked to Alzheimer’s disease [110]. Plxna4, which has been shown to regulate synaptic plasticity [111] and dendrite morphogenesis in the hippocampus [