Introduction

The positive effects of exercise in maintaining both physical and cognitive health have been widely established during the past few decades. In addition to alleviating the risk of metabolic diseases [1,2,3], cardiovascular disease [1, 4, 5], and several types of cancers [6], the benefits of physical exercise on brain function have also been extensively studied and reviewed [7,8,9,10]. Numerous studies in model organisms have shown that exercise promotes significant structural and functional changes in the hippocampus, which is a region of the brain that is essential for learning and memory processes. These changes indicate an overall positive impact of exercise on hippocampal function, including protective effects against aging and neurological disorders [11,12,13,14], increased adult neurogenesis [7, 15], enhanced synaptic plasticity, and improved cognitive performance [16]. Multiple studies have also demonstrated these beneficial effects in humans, indicating that aerobic exercise leads to an increase in the size of the human hippocampus, as well as improvements in memory performance [17,18,19].

While the molecular mechanisms underlying these effects are not fully understood, exercise was shown to increase the production of several neurotrophic factors [20,21,22] and other signaling molecules, such as neurotransmitters and cytokines [23,24,25,26]. Exercise was also found to induce changes in gene expression measured via quantitative polymerase chain reaction (qPCR), microarray, or RNA sequencing. For example, exercise has been found to upregulate the expression of genes that are involved in cell proliferation and differentiation, as well as genes that are involved in enhancing neuronal and synaptic plasticity [16, 27,28,29,30]. Many of these studies have focused on rodent models and have used running wheels as a behavioral paradigm to test the effects of physical exercise.

Recent studies have demonstrated that single nuclei RNA sequencing (snRNA-seq), which is a powerful next-generation sequencing (NGS) technique and has several advantages over traditional RNA sequencing methods (such as enabling transcriptomic analysis of tissues that are difficult to dissociate or have low cell yields), enables detection of previously unknown cell types and transcriptional subtypes [31,32,33]. Furthermore, it can be leveraged to better understand adult neurogenesis at the single-cell level in both mice and humans [31, 34]. Although a number of previous studies have looked into the transcriptomic and epigenetic changes induced by exercise on the brain [11, 35,36,37,38,39,40,41], to our best knowledge, none have investigated hippocampal cell-type specific transcriptomic changes following 4 weeks of voluntary exercise in mice without employing additional behavioral tests. In this study, we used running wheels as a form of voluntary exercise and aimed to further elucidate exercise-related mechanisms in the hippocampus using snRNA-seq.

Our results indicate that exercise selectively affects the proportions of cells in specific clusters of neuronal cell-types and suggest that the transcription factor Prdm16 may play a key role in orchestrating the maturation of newborn dentate gyrus neurons upon exercise. In addition, our data hint to an important role in crosstalk between NF-κB, Wnt/β-catenin, Notch, and retinoic acid (RA) signaling pathways in a specific excitatory neuron cluster belonging to the Cornu Ammonis (CA) region, which could collectively regulate and enhance synaptic plasticity in these neurons upon exercise. In conclusion, our findings provide an important novel dataset and help to further elucidate the molecular processes associated with improved neuronal plasticity in response to aerobic exercise.

Material and Methods

Animal Cohort

The experiments were performed according to the ethical guidelines of the national animal protection law and were approved by the responsible Ethical Committee of the State of Lower Saxony, Germany. Male C57BL/6 J mice were purchased from Janvier and housed individually in an animal facility with a 12-h light–dark cycle at constant temperature (23 °C) with ad libitum access to food and water. Animals were subjected to voluntary running at 3 months of age. For this, mice were kept in single cages with free access to the running wheel 24 h/day [42]. In the control group, the running wheel was blocked.

Sample Preparation, Single Nuclei Isolation, Library Preparation, and Sequencing

For RNA sequencing analysis, 4 runners and 4 control mice at the age of 3 months were assessed. Sample preparation, isolation of nuclei, and library preparation for sequencing were performed according to previously published protocols [43]. After 4 weeks in cages with or without running wheels, mice were sacrificed via cervical dislocation, the whole brain was removed from the skull, and the preparation of the hippocampus was performed in a petri dish filled with ice-cold PBS buffer under a dissection microscope (Leica M60, Leica Microsystems). Please note that both hippocampal hemispheres were collected from each mouse which resembles one sample. The tissue was immediately frozen on liquid nitrogen and stored at − 80 °C until further processing. Nuclei from frozen mouse brain tissues were isolated according to the previously published protocol with certain modifications [44]. Briefly, frozen tissues were Dounce homogenized in 500 µl EZ prep lysis buffer (Sigma NUC101) supplemented with 1:200 RNAse inhibitor (Promega, N2615) for 45 times in a 1.5 ml Eppendorf tube using micro pestles. The volume was increased with lysis buffer up to 2000 µl and incubated on ice for 5 min. Lysates were centrifuged for 5 min at 500 × g at 4 °C and supernatants were discarded. The pellet was resuspended into 2000 µl lysis buffer and incubated for 5 min on ice. After 5 min of centrifugation (500 × g at 4 °C), the resulting pellet was resuspended into 1500 µl nuclei suspension buffer (NSB, 0.5% BSA, 1:100 Roche protease inhibitor, 1:200 RNAse inhibitor in 1 × PBS) and centrifuged again for 5 min (500 × g at 4 °C). The pellet was finally resuspended into 500 µl NSB and stained with 7AAD (Invitrogen, Cat: 00–6993-50). Single nuclei were sorted using BD FACS Aria III sorter. Sorted nuclei were counted in Countess II FL Automated Counter. Approximately 4000 single nuclei per sample were used for GEM generation, barcoding, and cDNA libraries according to 10X Chromium Next GEM Single Cell 3ʹ Reagent v3.1 protocol. Pooled libraries were sequenced in Illumina NextSeq 550 in order to achieve > 50,000 reads/nuclei.

Preprocessing and Exploratory Analysis of snRNA Sequencing Data

The raw sequencing data generated from the snRNA-seq libraries were demultiplexed with cellranger mkfastq using CellRanger software v.4.0.0 (10XGenomics). Reads were aligned to the mm10 genome (NCBI:GCA_000001635.8) (GRCm38.p6) to obtain gene counts. The cellranger count pipeline was used with default options to generate a gene-count matrix by mapping reads to the pre-mRNA reference to account for unspliced nuclear transcripts.

All following downstream analyses were performed in R (version 4.1.0). The feature-barcode matrices generated for all 8 samples by CellRanger were converted to sample-wise Seurat objects using the Run10X() and CreateSeuratObject() functions from the Seurat package (version 4.0.2) [45]. For each sample, the object was filtered to remove nuclei where less than 200 genes and/or 500 reads and/or greater than 50,000 reads were detected. Genes that were expressed in fewer than 10 nuclei were filtered out. Moreover, any nuclei expressing more than 0.1% mitochondrial gene counts (“percent.mt”) were also removed. Normalization was done using the SCTransform() function from the Seurat and sctransform (version 0.3.2) packages [46], with the arguments method = “glmGamPoi” and vars.to.regress = “percent.mt.” PCA embeddings were calculated using the RunPCA() function, and RunHarmony() function from the harmony package (version 0.1.0) [47] (Korsunsky et al., 2019) was used to remove confounding effects due to sequencing batches. Leiden clustering [48] was performed using the FindNeighbors() and FindClusters() functions in Seurat with 50 principal components and a resolution parameter of 0.2. Uniform Manifold Approximation and Projection (UMAP) [49] was used for the visualization of clusters with the following arguments as parameters: umap.method = “umap-learn,” metric = “correlation,” n.neighbors = 30, n.components = 2, min.dist = 0.3, spread = 0.5, reduction = “harmony,” dims = 1:50. For each sample object, the DoubletFinder package (version 2.0.3) [50] was used to detect and remove any potential doublets. All individual sample objects were then merged into one Seurat object, which was re-normalized, integrated, and re-clustered using the SCTransform, RunPCA(), RunHarmony(), FindNeighbors(), and FindClusters() functions, with the same parameter values as before. UMAP was again used for the visualization of clusters in the final merged object. Cell-type annotation of clusters was done using the expression of selected marker genes (as listed in Supplementary Table 1. Three small clusters (22, 23, and 24) identified after the initial clustering (Supplementary Fig. 1A) were removed from the dataset before further analysis, since they did not show specific expression of any cell-type markers and subsequently could not be assigned to any known cell-type.

Cell-Type Abundance and Prioritization Analyses

The differences in the proportion of exercise and control cells in each cluster were calculated using the scProportionTest R package (https://github.com/rpolicastro/scProportionTest) (version 0.0.0.9), which calculates the p-value of the magnitude difference of abundance for each cluster using a permutation test (no. of permutations = 1000) and generates a confidence interval for the same using bootstrapping. Clusters with an observed fold difference magnitude of at least 1.5 with an FDR threshold of 0.05 were considered to have significant changes in cell-type proportions between exercise and control conditions.

To identify clusters (cell-types) in which the response to perturbation because of the given stimulus, i.e., exercise, is prioritized, we used the Augur R package (version 1.0.3) [51] [52].

Augur is a method that quantifies the responsiveness of each inferred cell type to a given treatment; exercise in this case. In a nutshell, the method trained 21 independent classifiers (one classifier per inferred cell type in our dataset) using as input features the gene expression data of the cells belonging to exercise and control groups. Augur ranked the discriminative ability of each of those 21 classifiers by using the mean cross-validated AUC as a criterion; with the general rule that the higher the AUC, the higher the priority.

Differential Gene Expression and Regulon Activity Analyses

The top gene markers for each cluster, as well as genes differentially expressed between exercise and control cells in each cluster, were identified using the FindMarkers() function in Seurat with the following parameters: logfc.threshold = 0.1, test.use = “MAST,” latent.vars = c(“batch,” “percent.mt”), assay = “RNA.” We adopted the single-cell regulatory network inference and clustering (SCENIC) workflow [53, 54] to construct gene regulatory networks from our data and evaluate the activity of these gene regulatory networks in cells from different cell-type clusters, using the SCENIC R package (version 1.3.1). Briefly, this analysis was performed as follows: with the gene expression matrix as input, gene sets that were co-expressed with potential regulators or transcription factors (TFs) were identified using the GENIE3 algorithm (implemented using the GENIE3 R package (version 1.16.0} [55]. Using the RcisTarget R package (version 1.14.0) [53], cis-regulatory motif analysis was performed to further exclude false positives and indirect targets of TFs and get direct-binding targets. With this step, only those co-expression modules that showed significant motif enrichment of a potential regulator (TF) were selected, and these final modules of TFs and their putative downstream target genes were referred to as “regulons.” Finally, the AUCell algorithm [53] was used to quantify the activity of these regulons in cells by assigning an enrichment or activity score to each regulon. This was calculated as an “area under the recovery curve” (AUC) across the ranking (based on expression values) of all genes expressed in a cell. The top 5 regulons specifically active in a selected cluster (compared to other clusters) were identified using the calcRSS() and plotRSS_oneSet() functions in the AUCell R package (version 1.16.0). The differentially active regulons across exercise and control cells within a cluster were determined using the FindMarkers() function in Seurat by setting the argument assay = “AUC,” where the “AUC” assay slot in the Seurat object contained the regulon AUC score matrix generated by AUCell.

For functional annotation and enrichment analysis of selected differentially expressed genes or differentially active regulons, an overrepresentation analysis was carried out using one of the following methods: clusterProfiler R package (version 4.0.0) [56], WebGestalt R package (version 0.4.4) [57], and the ShinyGO web application [58]. Gene regulatory relationship networks were constructed using the output of the SCENIC workflow and the Cytoscape software (version 3.9.1) [59].

RNA Velocity Analysis and PAGA Velocity Graph Generation

To understand cellular dynamics at a deeper level, we used the scVelo Python toolkit (version 0.2.4) [60] with default parameters on our dataset, which utilizes splicing kinetics to calculate RNA velocity estimates for single cells. Using these velocity estimates, we used the PAGA graph abstraction method [61] within the scVelo workflow to infer connectivities and directed transitions between selected cell-type clusters in our data. A PAGA graph was then visualized by associating a node with each cluster of interest and connecting these nodes based on the statistical measure of connectivity or transition between them.

Results

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
figure 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 

Full size image

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
figure 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

Full size image

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
figure 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)

Full size image

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 [112], is upregulated in ExN1 and ExN6. Zbtb16, upregulated in oligodendrocytes and ExN7, is involved in neural progenitor cell proliferation and neuronal differentiation during development [113]. Fkbp5 is a glucocorticoid receptor (GR) binding protein, which negatively regulates GR, and is interestingly upregulated upon exercise only in oligodendrocytes and microglia. It has been associated with mental disorders as well as in regulating neuronal synaptic plasticity [114]. Among the downregulated genes, Ddx5 is an RNA helicase that has been shown to regulate neurogenesis and neuronal differentiation, likely by inhibiting reprogramming to pluripotency [115, 116]. Tcf4, downregulated in ExN1 and InN1, has been linked to neurodevelopment [117], regulation of dendritic spine morphology [118], and regulation of neuronal migration through Bmp7 [119]. Altered Tcf4 function has also been linked to schizophrenia and mice in which Tcf4 is overexpressed develop cognitive impairment, which is attenuated by exercise [120].

Fig. 4
figure 4

Differentially expressed genes in excitatory neurons suggest regulation of synaptic plasticity and neuron differentiation to promote increased neurogenesis upon exercise. A (top panel) and B (top panel) UpSet plots depicting A commonly upregulated and B commonly downregulated genes upon exercise, between different clusters. Each row corresponds to a cluster, and bar charts on the right show the size of the set of genes up/downregulated upon exercise in that cluster. Each column corresponds to a possible intersection (commonly up/downregulated gene(s)): the filled-in cells show which cluster set is part of an intersection. Gene names are labeled for intersections with < 2 common genes among clusters. A (bottom panel) and B (bottom panel) Dot plots showing GO biological process terms enriched among A genes significantly upregulated and B genes significantly downregulated upon exercise, in clusters with significantly enriched terms. C Heatmap depicting the average regulon activity scores (RAS) for regulons with significant differential activity between cells from exercise and control conditions, in specific clusters. The number of genes comprising each regulon is indicated within parentheses. Regulons within red boxes indicate those upregulated upon exercise, and the ones within blue boxes indicate the ones downregulated upon exercise. D Cell-type specific regulons for ExN6 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 ExN6 are labeled on the plot, with the number of genes comprising each regulon indicated within parentheses. E Violin plots showing the average module score (expression) of the 3 regulons specific to ExN6 that also show upregulated activity upon exercise (Hes1, Klf8, Nr2f2). F 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). G Network plot with TFs and their respective regulon genes, specific to or showing differential activity in ExN6 upon exercise, indicating gene regulatory interactions in ExN6 cluster (genes/TFs that are significantly deregulated upon exercise or are specifically overexpressed in ExN6 are represented bigger in size, while selected genes belonging to these regulons which are deregulated in the same direction as the regulon TF are colored in red or blue)

Full size image

We further identified regulons that significantly differ in their activity across the exercise and control cells in certain clusters (Fig. 4C). Functional enrichment analysis of these regulon genes suggests that these TFs regulate a variety of biological processes in selected cells upon exercise (Supplementary Fig. 7 and Supplementary Table 21). Since ExN6 clearly exhibits a larger number of regulated genes as well as regulons with differential activity upon exercise, we sought to characterize these neurons by looking at the specific transcription factors active in this cluster, using the SCENIC workflow. The top 5 ExN6-specific regulons included Hes1, Klf8, Nr2f2, Fosl2, and Sox5 (Fig. 4D and Supplementary Table 22), and the top 3 regulons (Hes1, Klf8, and Nr2f2 regulons) also showed significant upregulation upon exercise in the differential regulon activity analysis as shown in Fig. 4C, further corroborated by increased overall expression of these regulon genes in the ExN6 exercise cells (Fig. 4E and Supplementary Fig. 7). The functional enrichment analysis of the combined list of genes comprising the top 5 regulons suggested regulation of timing of neuron differentiation (mediated by Hes1 and Sox5) as the term with the highest enrichment (Fig. 4F and Supplementary Table 23). Additionally, Hes1 is an important effector gene of the Notch signaling pathway, which has been linked to learning and memory by regulating synaptic plasticity in mature neurons [121, 122]. Interestingly, Hes1 was also induced in the early (0–1 h) time window upon acute exercise in skeletal muscle tissue [123], and another study showed that the inactivation of Hes1 in excitatory neurons resulted in abnormal fear and anxiety behaviors [124]. Klf8 is a TF playing critical roles in various biological processes such as the regulation of cell cycle progression. Elevated protein levels of Klf8 were shown to induce the activation of the Wnt/β-catenin signaling pathway, by promoting the expression and stabilization of β-catenin which further interacts with and inhibits NF-κB to regulate synaptic plasticity and cognition [125]. Moreover, Wnt/β-catenin signaling has been strongly implicated in the regulation of synaptic assembly, neuroprotection, neurotransmission, and synaptic plasticity, while disruption or deregulation of Wnt signaling has been linked to several neurodegenerative diseases such as AD [125] (Oliva, 2013). Nr2f2 (or COUP-TFII) belongs to the COUP-TF transcription factor family which is involved in regulating neuronal differentiation and migration of neural progenitors [126], and was additionally observed to be a member of the Klf8 regulon from the ExN6 SCENIC analysis. Interestingly, the regulon for retinoic acid (RA) receptor Rarb showed decreased activity in exercise cells, along with the regulons for Lhx6, Foxp2, and Pbx3 (as seen in Fig. 4C). Rarb is a direct transcriptional target of RA signaling, and this suggests its possible role in regulating neuronal plasticity upon exercise, which could be linked to the already-established roles of RA signaling in neuronal differentiation [127].

These findings suggest a pattern of crosstalk between NF-κB, Wnt/β-catenin, Notch, RA, and other signaling pathways, which could collectively regulate and enhance synaptic plasticity specifically in ExN6 neurons upon exercise. These relationships are visualized and highlighted in the ExN6 gene regulatory network (Fig. 4G; to explore the network in high magnification please view Supplementary Fig. 8a).

Discussion

In this study, we aimed to investigate the effects of physical exercise, in the form of voluntary wheel running, on the mouse hippocampal transcriptome at the single-cell level. Voluntary wheel running is an established approach in rodent models to test the effects of aerobic exercise on the brain [7, 15, 128, 129]. We decided to investigate the hippocampal transcriptomic changes at the 4-week time point after exercise since this experimental duration has commonly been employed in similar studies and has been shown to be sufficient in enhancing neuronal plasticity and would moreover allow us to capture specific states of adult neurogenesis [35, 40].

Through our analysis, we first observed that 4 weeks of exercise training leads to changes in the proportions of specific cell-types, most significant and pronounced in two clusters: a reduction in cells in the inhibitory cluster InN2 and an increase in cells in the excitatory cluster ExN11. Further investigation into the gene expression patterns in InN2 indicated that these cells, almost all originating from the control (sedentary) samples, showed high expression of the TF Prdm16. Prdm16 has already been studied in the context of cell-fate switch regulation in brown adipose tissue, but interestingly, many recent studies have implicated its potential role in neural development and differentiation, particularly in stem cell homeostasis and maintenance and positional specification of neurons [66,67,68,69,70, 130]. For example, one study showed that Prdm16 in mouse medial ganglionic eminence (MGE) progenitors plays an essential role in controlling the timing of maturation of forebrain GABAergic interneurons, with Prdm16 expression promoting the proliferation of interneurons and repressing maturation [70]. The expression of other top gene markers for InN2 indicated genes such as CaCCs Ano1 and Ano2 regulating ion transport, neuronal excitability, and plasticity, and zinc-finger protein encoding genes Zfhx3, Zic1, Zic4 and Zfp521 involved in neuronal and stem-cell differentiation. We further looked into the transcriptional regulation of InN2 by identifying highly active regulons in this cluster, with the results again implicating Prdm16 and its target genes, which showed enrichment of GO terms involved in neuronal maturation and adult neurogenesis, such as rhythmic process (Ankfn1, Hs3st2, Rorb, Zfhx3) and regulation of BMP signaling (Prdm16, Gpc3, Zfp423). This was substantiated by the highly specific expression of the Prdm16 gene in InN2, suggesting that the observed loss of exercise cells in InN2 is linked to the high expression of Prdm16 in cells from the control samples. Moreover, Dlx1, another highly active InN2 regulon that plays a role in regulating the functional longevity of GABAergic interneurons in the hippocampus, regulates Prdm16 and drives the expression of Prdm16-regulated Meis2 gene. Interestingly, a recent study also indicated that high Prdm16 expression is essential in promoting the disappearance of radial glia and the ending of cortical neurogenesis in the postnatal mouse brain [131]. Collectively, these findings, supported by the existing literature, allowed us to hypothesize that InN2 cells characterized by high Prdm16 expression could represent a specific developmental stage during adult neurogenesis captured only in cells from sedentary mice, possibly as a result of faster maturation of these Prdm16-high neurons upon exercise.

It is important to note that this interpretation is speculative. It would be necessary in future experiments to validate our data via additional approaches including RNAsope or spatial transcriptomics. Especially, the visualization of Prdm16 in hippocampal cells would be of interest. Moreover, our data remains descriptive, and future perturbation experiments would be essential to understand the neurobiological meanings of our observations.

Nevertheless, our interpretation is further supported when we investigated the NeuN-positive cluster ExN11 that showed an increase in the number of cells after exercise. Through RNA velocity measurements and PAGA graph generation, we observed a potential cell-type transition (in the context of trajectory inference) from the ExN11 cluster to the mature DG cluster ExN1, suggesting that ExN11 may represent newborn and maturating excitatory neurons that eventually integrate with the mature granule cells in the DG. ExN11 also seemed to express both neuronal and glial gene markers, which is consistent with studies that characterized radial glia-like precursor cells with astrocytic as well as neural stem cell-like properties during neurogenesis [91, 132].

Upon exploring the gene regulatory landscape for ExN11, we further identified specifically active regulons such as Meis1, Cebpd, Jun, Zfhx2, and Foxo1, all of whom are involved in neurogenesis, neuronal as well as glial differentiation, morphogenesis, and development.

These data may suggest that ExN11 neurons could represent specific states of adult-born and maturing neurons upon exercise. This view would be in agreement with previous data showing that adult neurogenesis in the dentate gyrus can be reactivated through stimuli that alter the levels of specific blood born factors that affect adult neuronal stem cells [133,134,135,136]. However, care has to be taken when interpreting these data that remain of course descriptive at this stage. Without further mechanistic studies, this interpretation of our data remains speculative. In addition, it should also be mentioned that the ExN11 cluster is also characterized by the expression of myelin-associated genes such as Mbp. While this may be attributed to the fact that these cells represent maturing but not finally differentiated newborn neurons, it was also shown that physical exercise can lead to increased myelination as a result of oligodendrocyte proliferation [137]. Although we do not see evidence of increased oligodendrocyte precursors or mature oligodendrocytes upon exercise, it would be premature to exclude the possibility that some of the changes observed indicate altered myelination. Future research should further explore this topic.

In summary, these data may indicate that the downregulation of Prdm16 in maturing new born hippocampal neurons could be an important regulatory step by which exercise enhances the process of adult neurogenesis. It could therefore be interesting to test strategies that reduce Prdm16 function in the hippocampus as a therapeutic approach to improve cognitive function.

Finally, we also looked at the broad gene expression patterns in the dataset that differed between exercise and sedentary conditions in specific cell-types. Although most clusters only showed very mild gene expression changes between the two conditions, ExN6 (belonging to CA excitatory neurons) exhibited a markedly high number of regulated genes as well as regulons with differential activity upon exercise, prompting us to explore it further. The most differentially active as well as specific regulon in ExN6 was Hes1, an important effector gene of the Notch signaling pathway, which itself regulates learning and memory by regulating synaptic plasticity in mature neurons [121, 122]. GO term analysis indicated a role of Hes1 in regulating the timing of neuron differentiation, mediated along with another specific ExN6 regulon Sox5. Multiple studies have shown that oscillating Hes1 expression regulates the maintenance and proliferative capacity of neuronal progenitor cells, acting as a repressor of the commitment of the latter to a neuronal fate, which in turn influences the number of neurons produced, the timing of neurogenesis as well as overall brain morphogenesis [138, 139]. Interestingly, Hes1 was also induced in the early (0–1 h) window upon acute exercise in skeletal muscle tissue [123], further indicating the possibility of an oscillating pattern of Notch signaling activation upon exercise. Our findings suggest that Hes1 upregulation (and subsequent Notch activation) in ExN6 at the 4-week exercise time point could mark a late maturation state of adult neurogenesis in adjacent DG neurons, potentially through lateral inhibition of neuronal fate commitment of newborn cells. This is in line with our previous hypothesis that Prdm16 upregulation in InN2 neurons at this time could also be associated with an important maturation step in adult neurogenesis in the DG and that ExN11 neurons represent a maturing state of excitatory cells that are past the neurogenesis phase. Analysis of other top gene markers and regulons of ExN6 (Klf8, Nr2f2, Fosl2, and Sox5) suggested the involvement of pathways regulating neuronal and synaptic plasticity, such as Notch signaling, retinoic acid (RA) signaling, and Wnt/β-catenin signaling, which is further linked to inhibition of NF-κB to regulate synaptic plasticity and cognition. Together, these findings point to a unique pattern of crosstalk between these various signaling pathways, suggesting that they act collectively to regulate and enhance synaptic plasticity in ExN6 neurons upon exercise.

This study of course comes with its own limitations, such as technical drop-outs in the snRNA-seq dataset that limit gene expression analyses. We also look solely at a fixed snapshot of changes occurring at the 4-week time-point after exercise, and it would be beneficial to consider a longer range and/or multiple time-points for temporal changes that occur with these experimental paradigms, especially when it comes to gene expression, which would give a better insight into the dynamic processes we highlight with our study. Additionally, the different neuronal clusters/subtypes and their interactions could be more deeply characterized through functional studies in the future. Another issue to be mentioned is the fact that we detected comparatively few astrocytes [140]. This might be due to the specific nuclei isolation protocols [141, 142] suggesting the necessity to address the role of such cells in future single-cell/nucleus RNAseq experiments. In addition, an epigenome analysis at the single-cell level, for example via ATAC-seq, could elucidate to what extend epigenetic processes, that have been implicated with exercise-improved learning behavior [10, 11], may contribute to the observed changes in gene-expression. Finally, our study is currently limited to male mice but we have initiated similar experiments in female mice as part of the European JPND project EPI-3E.

Despite the fact that our findings are currently limited to interpretations from bioinformatic analyses and a single snRNA-seq dataset, we believe that they encourage further research focusing on functional interpretation and validation of these results and probe important questions in the field of aerobic exercise and its impact on brain health.

In conclusion, our study provides a useful single-cell sequencing resource for exploring exercise-induced transcriptomic changes in the adult mouse hippocampus. Some of the observed changes allow for a deeper insight into the molecular processes that control adult hippocampal neurogenesis. For example, we provide evidence that a cluster of neurons characterized by high Prdm16 expression may represent a specific state in the maturation of adult-born hippocampal neurons, and our data suggest that a decrease of Prdm16 levels is critical to eventually promote the maturation of newborn neurons in the adult brain. We also identified a cluster of excitatory CA neurons exhibiting pronounced gene expression changes leading to regulation and enhancement of synaptic plasticity in the hippocampus upon exercise. Since adult hippocampal neurogenesis is impaired during aging and neurodegenerative diseases such as Alzheimer’s diseases, targeting Prdm16 could be a novel approach to enhance plasticity in the diseased brain.