Introduction

Physical exercise is a potent modulator of central nervous system plasticity, promoting structural and functional adaptations across diverse neuronal networks. While extensive research has characterized exercise-induced changes in the brain, particularly within hippocampal circuits (1–5), the spinal cord represents a remarkably understudied yet critical component of the exercise-neuroplasticity axis. The spinal cord integrates descending motor commands with ascending sensory information, positioning it to coordinate exercise-related neural adaptations (6).

The therapeutic potential of exercise-induced spinal plasticity has been predominantly studied in injury contexts. After spinal cord injury, exercise promotes neurotrophin expression, enhances synaptic transmission, and facilitates oligodendrocyte proliferation (7–13). Voluntary wheel running is particularly effective compared to forced paradigms, promoting enhanced neurotransmitter modulation and axonal innervation (12, 14–17). These studies implicate pathways such as brain-derived neurotrophic factor (BDNF) signaling, synaptic vesicle cycling, and glial-mediated responses in these processes (7–9, 11, 15, 17, 18). However, a critical knowledge gap exists regarding exercise effects in the healthy spinal cord. Understanding these foundational adaptations is essential for various reasons, e.g. to understand the consequences of increased and repeated (hyper-) activation of motor systems in exercise (often over hours and to failure), thereby revealing neuroprotective processes that may be leveraged in age- or disease-related degeneration.

A multicellular involvement in the response of the spinal cord to repeated bouts of exercise could be expected from prior data in the spinal cord as well as insights gained in the brain. Glial cells demonstrate exquisite sensitivity to activity-dependent stimuli (4, 10, 12, 16, 19–22). Microglia orchestrate neuroplasticity responses and coordinate cellular interactions during repair, with exercise promoting anti-inflammatory responses while suppressing pro-inflammatory activation, and playing additional roles, e.g. in mediating neuropathic pain (23–27). Astrocytes in the brain respond to physical activity through metabolic and structural adaptations that support neuronal function (28–33), but their specific responses in the healthy spinal cord remain largely unexplored (34). Oligodendrocytes are also highly exercise-sensitive, with voluntary running enhancing precursor proliferation through Wnt signaling and increasing nestin-positive progenitors in spinal white and gray matter (21). Motor skill learning specifically requires rapid oligodendrocyte production, and exercise promotes peroxisome proliferator-activated receptor γ coactivator 1α (PGC1α)-mediated myelin synthesis, suggesting these cells facilitate adaptive myelination in response to increased neural activity (12, 13).

Recent advances in single-nucleus RNA sequencing (snRNA-seq) have revolutionized our understanding of cellular heterogeneity within spinal cord circuits, revealing extensive glial diversity in the human and mouse spinal cords (35–44). These approaches overcome technical barriers of adult tissue dissociation and enable unbiased resolution of exercise-responsive populations. When integrated with proteomics, they provide comprehensive insight into both cell-type-specific transcriptional dynamics and tissue-wide protein-level adaptations. Two recent studies applied snRNA-seq to examine how aerobic exercise influences the aging spinal cord, showing that training can partially restore youthful cellular and regulatory features (45, 46). While these datasets advance our broader understanding of exercise–aging interactions, several features of exercise-induced adaptation in the healthy adult spinal cord remain insufficiently defined. In particular, neither study captured acute temporal dynamics following an exercise bout, assessed inter-individual variability due to pooled sequencing designs, or examined lumbar segments that directly innervate hindlimb musculature. Additionally, both relied on fresh-tissue dissociation, which is susceptible to stress-induced transcriptional artifacts that newer fixed-tissue approaches help to mitigate.In this study, we use multi-omic profiling to characterize how endurance training and acute exercise remodel the healthy mouse spinal cord. We combine voluntary wheel running with bulk proteomics and fixed-tissue snRNA-seq (10x Genomics Flex) to profile lumbar spinal cord segments at baseline and at 6 h and 24 h after an acute bout of exercise in both sedentary and trained mice. This design allows us to distinguish transient from long-lasting adaptations, quantify cell-type–specific transcriptional programs and signaling pathways with independent biological replicates, and focus particularly on glial plasticity and neuron–glia crosstalk. Our work therefore provides a temporally and anatomically resolved, multi-omic framework for exercise-induced adaptations in the young, healthy spinal cord.

Results

Endurance Training Remodels Spinal Cord Composition and Glial Gene Expression

To investigate how prolonged physical activity impacts the spinal cord, we subjected adult male C57BL/6 mice to six weeks of voluntary wheel running (TRAIN), with sedentary littermates serving as controls (SED) (Figure 1A). Whole-body composition and exercise capacity were assessed before and after the intervention. Mice with access to running wheels voluntarily ran between 4 and 8 km/day, averaging approximately 6 km/day over the training period (Figure 1B). While total body mass increased in both groups (Figure S1A), trained mice preserved lean mass, whereas sedentary mice showed significant decreases in relative lean mass combined with an increase in fat mass (Figure S1B–C). Treadmill-based endurance testing demonstrated an approximately threefold increase in maximal exercise capacity in trained mice (Figure 1C), validating the effectiveness of our aerobic training protocol.

Figure 1.

To determine whether the endurance training paradigm induces molecular remodeling in the spinal cord at rest, i.e. untainted by the last acute exercise bout, we performed quantitative bulk proteomics on lumbar spinal cord lysates. 159 proteins were differentially expressed between trained and sedentary mice (p < 0.05, |log₂FC| > 0.5); including 97 upregulated and 62 downregulated proteins (Figure 1D). Gene Ontology (GO) enrichment analysis revealed upregulated proteins to be linked to mitochondrial metabolism, vesicle-mediated transport, autophagy, and stress response pathways, whereas downregulated proteins were enriched for chromatin organization and immune-related functions (Figure 1E). Top upregulated proteins—including TIMM22, NDUFA11, UBQLN4, MAP1LC3B, DYNLL2, TUBB4A, NRF1, ETV3, and FABP4 —reflect regulation of neurodevelopmental signaling, chromatin remodeling, cytoskeletal and metabolic adaptation, and CNS homeostasis (Figure 1F).

To resolve the cellular origin of these proteomic changes, we performed high-throughput single-nucleus RNA sequencing (snRNA-seq) using the 10x Genomics Fixed RNA Profiling (Flex) protocol (Figure 1A, S1D-E). This method enabled nuclei isolation from frozen adult central nervous system (CNS) tissue, providing improved recovery of intact transcriptional profiles, with enhanced sensitivity for low-abundance transcripts and greater uniformity across nuclei. Using canonical gene markers, we annotated 11 major cell-types, including neurons, astrocytes, oligodendrocyte lineage cells, microglia, vascular cells, and ependymal cells, with finer substructure within multiple glial populations (Figure 1G). Estimates of cell-type abundance remained broadly stable between groups (Figure 1H), indicating that exercise training does not alter gross cellular composition in a persistent manner. Nevertheless, transcriptional reprogramming was evident, with glial cells—particularly oligodendrocytes— showing the most prominent differential gene expression in response to training (Figure 1I). These data suggest glial populations play a central role in mediating spinal adaptations to chronic physical activity.

Glia are the Primary Transcriptomic Responders to Training

To identify the key cellular contributors to spinal cord remodeling post-endurance training, we focused on glial populations showing the most pronounced transcriptional changes: oligodendrocytes, oligodendrocyte precursor cells (OPCs), and the Astrocyte-1 subcluster (Figure 2A–C). Among these, oligodendrocytes had the greatest number of differentially expressed genes (DEGs). Gene set enrichment analysis (GSEA) indicated that shared biological processes—such as synaptic modulation, axon development, and intercellular signaling— were activated across these glial populations (Figure 2D). However, the specific genes driving these enrichments were largely cell-type specific, with only minimal DEG overlap (Figure S2A), except for the shared genes between oligodendrocytes and OPCs, consistent with their developmental relationship.

Figure 2.

In response to training, oligodendrocytes upregulated genes associated with vesicle trafficking and synaptic function (Slc17a7, Syt2, Snap25, Vgf), ion transport and membrane potential (Atp1a3, Slc4a4, Panx1), and axon-glia interactions (Ret, Dkk3, Sparcl1, Hapln4) (Figure 2E–F, S2B–C). Increased expression of metabolic and structural genes (Acsl6, Mlc1, Reln, Slc24a4) suggests enhanced activity-dependent support of axons. Trained OPCs were enriched for genes related to neurofilament organization, intracellular signaling, and synaptic modulation (Figure 2B), with upregulation of synaptic, metabolic, and structural genes (Figure 2G), consistent with emerging roles for OPCs in myelin plasticity and neuronal responsiveness (47, 48).

Astrocyte-1 displayed enriched programs related to gap junction regulation, cell–cell communication, and axonogenesis (Figure 2C), with significant upregulation of cytoskeletal and signaling regulators (Vim, Rtn1, Map1b, Hes5, Trf, Sparcl1) in trained compared to untrained spinal cord (Figure 2H). Although several DEGs (Gja1, Syp, Vgf, Nefh) were expressed in multiple glial subtypes, their expression trajectories differed, suggesting functionally distinct roles. These findings demonstrate that while glial subtypes converge on certain biological processes, their transcriptomic responses are highly cell-type specific and aligned with their specialized roles in the spinal cord.

Expanding the analysis to other cell-types, we found that microglia-2 displayed training-induced expression changes in genes associated with synaptic regulation and secretory function (Dnm1, Snap25, Stxbp1, Sparcl1, Nrxn1), as well as enriched GSEA terms for gliogenesis and synaptic modulation (Figure S2D, H). These data support a shift away from classical immune activation toward a more synaptically engaged, neuroprotective microglial phenotype.

By contrast, neuronal responses were more heterogeneous and generally of smaller magnitude. Excitatory neurons (ExN-1) showed modest upregulation of ion transport and homeostasis genes (Slc1a2, Atp1a2, Slc4a4, Gja1) and downregulation of immune-responsive genes (Nr4a1, Ly6a), suggesting reduced inflammatory signaling after training (Figure S2E, I). Inhibitory neuron subclusters responded variably: Inh-1 displayed upregulation of genes related to neuronal growth and cytoskeletal organization, while Inh-5 was enriched for genes in MAPK signaling, membrane transport, and neurotransmitter uptake (Figure S2F–G). Genes upregulated included Agt, Apoe, Gja1, Uts2, Mt1, Mt3, and various solute carrier transporters (Figure S2J–K), reflecting increased inhibitory tone and metabolic support.

Endurance Training Reshapes Intercellular Signaling in the Spinal Cord

To determine if transcriptional changes were accompanied by altered intercellular communication, we used CellChat (49) to infer ligand–receptor (L–R) interactions among spinal cord cell-types. Endurance training remodeled multiple signaling networks, suppressing stress- and inflammation-associated pathways such as neurotensin (NTS), leukemia inhibitory factor receptor (LIFR), and carcinoembryonic antigen-related cell adhesion molecule (CEACAM) signaling, while enhancing activity in pathways including Notch, non-canonical WNT (ncWNT), and Galectin (Figure S2L). Notably, the predicted communication patterns did not simply mirror the magnitude of differential expression within each cell-type: neurons (inhibitory, cholinergic, excitatory) and supportive cells (fibroblasts, astrocytes, vasculature) emerged as major senders, while astrocytes, neurons, and OPCs were the predominant receivers (Figure 1I, S2L–M). Thus, even populations with relatively few DEGs contributed disproportionately to signaling networks, highlighting that exercise-induced plasticity involves a reorganization of intercellular communication beyond what is captured by cell-intrinsic transcriptional changes alone.

Further analysis of cell-type-specific L–R usage across key signaling pathways (Figure 2I–J, S3) confirmed that several communication axes were enhanced by training. WNT, PTN, and NECTIN signaling from cholinergic neurons to glial and vascular targets were mainly observed only in trained animals (Figure S2I), suggesting recruitment of additional structural and trophic signaling cues. Likewise, NRXN-DAG1 interactions between cholinergic neurons and astrocytes were strengthened by exercise training, potentially enhancing synaptic formation, stability, and communication (Figure S3A). NRXN signaling from astrocytes to neurons, OPCs, and other astrocytes also increased (Figure S3B), indicating broader post-training engagement of astrocyte-mediated synaptic and glial networks.

Laminin signaling derived from astrocytes and oligodendrocytes increased toward neuronal, vascular, and fibroblast targets with exercise (Figure S3C–D), consistent with remodeling of extracellular matrix and tissue architecture. Microglia-2 in trained mice exhibited reduced JAM-mediated signaling and amplified Semaphorin (SEMA4) and FGF signaling to several targets (Figure 2J, S3E), marking a possible transition from immune surveillance toward growth- and guidance-related signaling. Together, these findings show that physical activity drives marked, cell-type-specific remodeling of spinal cord intercellular communication networks.

Acute Exercise Induces Distinct Temporal and Cell-Specific Molecular Responses Shaped by Training

Given the extensive transcriptional and proteomic remodeling evoked by endurance training at rest, we next sought to determine how the spinal cord temporally responds to a single acute bout of exhaustive exercise—and whether this response is modulated by prior training. Sedentary (SED) and trained (TRAIN) mice underwent a single session of treadmill-based running to exhaustion before lumbar spinal cords were collected at 6-hours (6h) and 24-hours (24h) post-exercise (Figure 3A).

Figure 3.

Bulk proteomics revealed a dynamic molecular response in sedentary mice, with modest changes at 6h and robust upregulation of 710 proteins by 24h (Figure 3B–D). In contrast, trained mice exhibited blunted proteomic responses at both time points, with fewer differentially expressed proteins and smaller fold changes (Figure 3D–E), implying that training buffers the acute molecular response and fosters more robust homeostatic control. Gene ontology enrichment analysis supported this interpretation. In sedentary mice, acute exercise strongly activated immune-related, metabolic, and mitochondrial pathways—particularly at 24h (Figure S4A–B). In trained animals, these same pathways were suppressed (Figure S4C–D), consistent with a shift from reactive to regulated responses.

To resolve the cellular contributors to these bulk changes, we performed single-nucleus RNA sequencing (snRNA-seq) on all six experimental conditions. Integration with prior datasets identified 11 major cell-types with stable clustering across samples (Figure 3F). Overall cell-type proportions remained largely unchanged, but there were subtle, dynamic shifts. Notably, Microglia-1 abundance decreased 24h post-exercise in sedentary animals, with a milder trend in trained mice (Figure 3H–I). Spatial overlays of nuclei suggest this reduction stems from a shift toward Microglia-2, potentially representing altered microglial activation (Figure 3G). Trained mice displayed a modest increase in oligodendrocytes at 24h and both groups showed a minor increase in inhibitory neuron subcluster 1 (Inh-1) after exercise (Figure 3G–I, S4F–G). These changes were absent at baseline, indicating they likely represent transient, exercise-induced transcriptomic shifts rather than chronic adaptive effects.

These data show that acute exhaustive exercise triggers rapid, extensive proteomic and cellular remodeling of the spinal cord—particularly in untrained animals. Prior training buffers both the magnitude and nature of these responses, suggesting that long-term activity impacts spinal cord adaptation.

Training State Differentially Tunes Glial Transcriptional Programs after Acute Exercise

To address how prior training modulates glial responses to acute exercise, we characterized transcriptional dynamics in spinal cord cell-types 6 and 24h after a single bout of exhaustive exercise. Four glial populations—oligodendrocytes, astrocyte-1, OPCs, and microglia-2— demonstrated extensive DEGs at both timepoints (Figure 4A–C, S5A). At 6h post-exercise, DEGs were predominantly upregulated, with trained mice showing more DEGs. By 24h, DEGs expanded in both groups with a mix of up- and downregulation, representing a biphasic response. Additional glial subtypes, including astrocyte-2 and microglia-1, were also activated with time-dependent patterns, notably with microglia-1 mainly responding at 24h (Figure S5B–C).

Figure 4.

Overlap analysis confirmed minimal DEG sharing between glial populations (Figure S5D–G), highlighting strong cell-specificity. While oligodendrocytes dominated the baseline transcriptional response to chronic exercise, the acute response strongly engaged astrocyte-1 and OPCs besides oligodendrocytes, particularly in trained mice. This re-weighting suggests that prior training tunes multiple glial types for transient activation following physiological challenge, reflecting a more distributed and responsive glial network.

GSEA demonstrated distinctive functional programs among these glia. Oligodendrocytes were enriched for metabolic and synaptic terms (Figure 4D), with transient upregulation of metabolic genes (Etnppl, Pla2g3, Plaat5, Itpr2, Cth) at 6h in both groups, with stronger induction in trained animals, suggesting enhanced metabolic readiness (Figure 4E). Synaptic genes (Syt1, Syt2, Dnm1, Gria2, and Glra1) were elevated at baseline in trained animals, remained stable at 6h, and declined at 24h, while sedentary mice showed lower baseline expression and delayed induction at 24h (Figure 4F). These divergent trajectories imply that training pre-establishes a synaptic support state, which undergoes rapid resolution post-exercise, whereas sedentary conditions have a delayed or dampened transcriptional response.

Astrocyte-1 regulated metabolic, synaptic, and signaling genes in a training-dependent manner (Figure 4G), including robust baseline or 6h induction and then suppression at 24h of lipid metabolic (Fabp7, Pnpla2, Srebf1) and cholesterol biosynthesis (Hmgcr, Cyp51, Mvd) genes (Figure 4H). Mitochondrial genes (Cox6c, Ndufa4) were induced at modestly higher levels in trained animals at baseline and 6h (Figure 4H). Genes for synaptic function and cytoskeletal remodeling (Dnm1, Map1a/b, Rtn1, Tubb3, Rab11a, Syn1, Syp, Sez6l2), and immediate early/neurotrophic genes (Vgf, Trf) were upregulated at baseline in trained mice, then normalized or suppressed by 24h (Figure 4H–I). In sedentary animals, these genes were unchanged or weakly induced at 6h, suggesting that training primes astrocytes for prompt activation and rapid resolution post-exercise.

OPCs exhibited fewer but distinctive changes in metabolic and signaling genes, most notably in trained animals (Figure S5H). At baseline, trained OPCs expressed higher levels of anabolic and cholesterol biosynthesis genes (Fdps, Msmo1, Hmgcs1, Sqle, Dhcr7/24), which were downregulated at 6h in TRAIN and at 24h in SED (Figure S5J). Other genes, including those related to glycolysis, redox, chromatin remodeling, and Wnt signaling (Aldoa, Phka1, Gpd1, Kdm4a, Tcf7l2), peaked at 6h or 24h in a training-dependent manner. These results indicate that OPCs coordinate anabolic and plasticity programs more efficiently after training.

Microglia-2 exhibited global suppression of classical homeostatic (P2ry12, Csf1r, Hexb) and immune (Tyrobp, C1qc, Ctss, Lamp2) genes at 24h, regardless of training state (Figure S5I, S5K). A transient upregulation of remodeling- and plasticity-associated genes (Adora1, Tmsb4x, Daam2) at 6h also suggests subtle priming toward a neuroprotective state in response to acute activity.

Training Modulates the Timing and Nature of Intercellular Communication after Acute Exercise

Applying CellChat L–R analysis across all spinal cord cell-types at baseline, 6h, and 24h post– acute exercise revealed that communication number and strength peaked at 6h and diminished sharply by 24h in both groups (Figure 5A–B), paralleling biphasic transcriptional regulation in glia (Figure 4A–C, S5A-C). Trained mice displayed slightly higher overall inferred interactions at baseline and 6h and lower inferred interactions at 24h.

Figure 5.

Pathway-level analysis identified distinct training-dependent remodeling of intercellular communication patterns (Figure 5C–D). In trained animals, positive remodeling was characterized by higher baseline or early activity of acetylcholine, PDGF, FLRT, SLIT, BMP, FGF, PTPR, and SEMA4 signaling, with several of these pathways (e.g., BMP, FGF, PDGF) showing early resolution by 24h. Conversely, sedentary animals preferentially engaged NTS, LIFR, CXCL, JAM, and sustained FGF/PTPR signaling, often peaking at 6h and persisting through 24h, indicative of a prolonged inflammatory/trophic state.

Cholinergic neuron–derived signaling showed marked training-dependent remodeling across multiple pathways (Figure 5E–G, S6). Guidance cues such as Semaphorin and Slit were strongly enhanced in trained animals: Sema4 signaling to astrocytes and vasculature increased particularly at 6h for astrocyte-1, at baseline and 6h for astrocyte-2, and at all timepoints for vascular targets, while Slit signaling rose from cholinergic neurons to both astrocytes and themselves (Figure 5E-F). In both cases, the elevated signaling strength reflected increased expression of receptors in the target populations (astrocytes and cholinergic neurons), rather than increased ligand expression in the sending neurons (Figure S6A–C). This shift suggests that glial cells, often the primary receivers, act as downstream transcriptional responders to neuronal cues. Within neuronal circuits, acetylcholine signaling dynamics further diverged: in trained mice, self-signaling among cholinergic neurons peaked sharply at 6h and resolved by 24h, whereas in sedentary mice it persisted through 24h. By contrast, cholinergic-to-inhibitory signaling was stronger in SED at 6h, while excitatory neuron targeting favored baseline activity in TRAIN (Figure 5G, S6D). Additional trophic pathways also showed accelerated resolution with training. FGF signaling was elevated in trained animals at baseline (toward cholinergic neurons and OPCs) and at 6h (toward excitatory neurons), with corresponding increases in Fgfr1/3 expression in OPCs (Figure S6E–F). By contrast, sedentary mice showed prolonged FGF activity in multiple, persisting at 24h post-exercise (Figure S6E-G). Similarly, PTPR signaling was enhanced at baseline and 6h in trained mice but diminished by 24h, whereas sedentary animals showed delayed and sustained activity at 24h (Figure S6H). Modest receptor upregulation (Ptprd, Ptprf, Ptprs) accompanied these changes (Figure S6F–G). The fine-tuned regulation of PTPR signaling in TRAIN likely reflects balanced OPC proliferation and maturation supporting efficient myelin adaptation.

Together, these results demonstrate that cholinergic neuron-driven signaling to glial cells is one of the most dynamically altered features of the trained spinal cord—both acutely and chronically—and represents a principal upstream driver of adaptive glial remodeling. Receptor-level tuning in glia likely enables the trained spinal cord to remodel rapidly and resolve signaling states more efficiently than in sedentary conditions.

Glial Populations Fine-Tune Communication Networks in Response to Exercise

Although cholinergic neurons emerged as central hubs for intercellular signaling after exercise (Figure 5, S6), glial populations—astrocyte-1, oligodendrocytes, and microglia-2—also undergo substantial, training-mediated communication remodeling, reflecting a nuanced regulatory role in shaping the spinal cord environment (Figure 6, S7).

Figure 6.

Astrocytes showed pronounced differences in overall signaling between training states. JAM signaling was stronger in sedentary animals, particularly at baseline and 24h, where astrocyte-1 and astrocyte-2 sent enhanced autocrine signals, to microglia-2, and to OPCs (Figure 6A). This was mirrored by higher expression of JAM receptors (F11r, Jam2, Jam3) in sedentary astrocytes at these timepoints (Figure 6B–D, S7A). By contrast, semaphorin signaling from astrocytes to other astrocytes and to vascular cells was increased at baseline and 6h post-exercise in trained animals (Figure S7C), mirroring patterns seen in cholinergic neurons (Figure 5E).

Other glial pathways also showed training-specific regulation. PDGF signaling, originating from astrocyte-1, microglia-2, and oligodendrocytes, was enriched in trained mice at baseline and 6h, with higher ligand–receptor expression in astrocytes and OPCs (Figure 6B, 6E–F, S7A–B). These interactions are associated with growth and differentiation programs, suggesting that training primes glia for adaptive remodeling. In contrast, FGF signaling from oligodendrocytes was elevated in sedentary animals at baseline and 6h and persisted at 24h (Figure S7D), whereas trained animals exhibited lower levels at 24h. The broad elevation of oligodendrocyte-derived FGF signaling across multiple spinal cord cell-types in sedentary animals suggests a reliance on prolonged trophic support and maintenance cues, whereas its earlier resolution in trained animals may reflect a shift toward more efficient, differentiation-driven remodeling.

Microglia-2 demonstrated divergent signaling programs, with SED animals showing broadly elevated JAM-mediated inflammatory signaling towards all cell-types at baseline that was resolved following exercise (Figure 6D, G). In contrast, trained microglia showed stronger engagement in growth- and guidance-related pathways—including SEMA4 and FGF—already at baseline; these signals were then transiently downregulated at 6h and 24h post-exercise, respectively. Notably, in SED animals, SEMA4 and FGF signaling emerged only at 6h post-exercise. This pattern indicates that even a single bout of exercise biases microglia away from inflammatory outputs and toward trophic and structural support.

In sum, glial populations act as fine-tuners of the spinal cord environment. While cholinergic neurons appear to provide primary “drive,” glia regulate the timing, precision, and qualitative balance of communication pathways. This dual mechanism allows exercise to orchestrate coordinated, cell-type–specific remodeling—where upstream neuronal cues initiate broad adaptation and glial cells enforce an adaptive balance between inflammation, structure, and trophic support.

Discussion

Physical activity exerts profound benefits on neural health, yet how the spinal cord, a critical hub for sensorimotor integration, adapts to exercise under healthy conditions has remained largely unexplored. Previous work has focused primarily on exercise in the context of injury, neurodegeneration, aging, or systemic disease, where exercise often mitigates pathology or enhances recovery(7, 9, 10, 13, 14, 16, 17, 19, 50, 51). As examples, two recent studies examined exercise effects in the aging spinal cord using snRNA-seq, demonstrating that aerobic exercise can partially restore youthful cellular signatures and attenuate age- associated microenvironmental decline (Du et al. 2025; Li et al. 2023).

By contrast, the present study provides a multi-omic characterization of how the healthy spinal cord remodels in response to an acute endurance exercise bout in a training state- dependent manner, as well as to chronic training at rest. Using bulk proteomics, fixed tissue snRNA-sequencing, and cell–cell communication inference in the lumbar spinal cord, we uncover transcriptional and signaling programs that engage predominantly in glial populations and are shaped by time post-exercise and training status. Our results highlight glial plasticity and cholinergic neuron–glia crosstalk as central mechanisms through which exercise remodels the spinal cord microenvironment.

Glia emerge as the primary transcriptional responders

Across analyses, oligodendrocytes, astrocytes, OPCs, and microglia displayed the most extensive transcriptional remodeling, whereas neurons exhibited more modest and heterogeneous responses, placing glia under the spotlight for mediating exercise-induced adaptations (Figure 7A). Functionally, the transcriptional programs engaged across these populations were strikingly convergent—enrichment for metabolic coordination, lipid turnover, synaptic signaling support, and structural remodeling. However, the specific gene programs were largely cell-type specific, underscoring tailored glial contributions (Figure 7A). These findings align with mounting evidence that glia are not merely passive supporters, but active modulators of circuit excitability, synaptic remodeling, and motor performance (52–55).

Figure 7.

Training state dictates the temporal dynamics of glial responses

By profiling both steady-state (baseline) adaptations after training at rest and acute responses at 6h and 24h following a single bout of exercise, we distinguished long-lasting remodeling from transient shifts. A clear biphasic pattern emerged: metabolic and signaling genes were induced at 6h and largely resolved by 24h. Crucially, trained animals exhibited earlier induction and faster resolution, whereas sedentary animals showed delayed and more persistent activation. For example, synaptic genes in oligodendrocytes and astrocytes were already elevated at baseline in trained mice, remained stable at 6h, and declined at 24h, while in sedentary animals, they were induced only at 6h for astrocytes and at 24h for oligodendrocytes. These temporal distinctions between trained and sedentary animals are summarized in Figure 7B, which illustrates representative gene programs across key functional categories and highlights how prior training shifts the timing, amplitude, and persistence of transcriptional responses. The training state-dependent re-weighting implies that prior exercise primes glia for rapid but transient responses, thereby reducing the persistence of injury-like or inflammatory phenotypes that characterize sedentary animals. Such dynamics likely underlie the enhanced resilience and adaptive capacity observed in trained spinal cords.

Neuron–glia communication is a central driver of adaptation

One of the more curious findings was that cholinergic neurons, despite showing few DEGs, emerged as dominant hubs of intercellular signaling. CellChat analysis revealed training- dependent remodeling of cholinergic neuron–glia pathways including Semaphorin, Slit, PDGF, FGF, and PTPR. In trained animals, these signals were engaged at baseline or within 6h and resolved by 24h, whereas in sedentary animals they were absent or persisted at 24h. These pathways converge on processes of axon guidance, synaptic adhesion, vascular remodeling, and myelin adaptation—precisely the themes transcriptionally enriched in glial populations (Figure 7A). The overlap strongly suggests that cholinergic neurons coordinate glial transcriptional programs through paracrine signaling, positioning cholinergic neuron–derived cues as upstream regulators of glial plasticity. This observation parallels recent work implicating motor neurons in shaping oligodendrocyte and astrocyte states via activity- dependent signaling (55–58).

Microglial remodeling complements astrocyte and oligodendrocyte adaptation

Microglia exhibited a distinct profile characterized by diminished inflammatory JAM signaling at baseline in trained animals alongside early engagement of growth- and guidance-related pathways such as SEMA4 and FGF. This was accompanied by transient induction of remodeling-associated genes such as Adora1, Tmsb4x, and Daam2. Functionally, this suggests a shift away from inflammatory activation and toward a more protective, growth-supportive phenotype. This pattern resonates with recent findings that exercise promotes anti- inflammatory microglial states in the brain (59–61) and extends this principle to the spinal cord.

Mechanistic implications and integration with existing literature

Our findings extend previous work showing that exercise enhances neurogenesis, angiogenesis, and myelin remodeling in the brain (1, 32, 57, 62–64) and injured spinal cord (7, 9–14, 16, 21, 50) by demonstrating that similar principles apply in the healthy spinal cord. The convergence of metabolic and synaptic programs across glial types supports their role as activity sensors and circuit stabilizers and is summarized in the cell-type–specific schematic in Figure 7A. The central contribution of cholinergic neurons highlights a model in which these neurons orchestrate adaptive plasticity indirectly, by modulating glial activity rather than relying on broad neuronal transcriptional changes. At the same time, our work complements recent aging-focused studies that applied snRNA-seq to examine how aerobic exercise influences spinal cord cell states (45, 46). These studies showed that exercise can partially restore youthful transcriptional trajectories in aged animals, but their designs did not address several key features of exercise adaptation in the healthy adult spinal cord. Specifically, neither study captured acute temporal dynamics after an acute exercise bout, assessed inter- individual variability due to pooled sequencing replicates, nor examined lumbar spinal segments directly involved in hindlimb locomotion. Both studies also relied on fresh-tissue dissociation, which can introduce stress-related transcriptional signatures that are minimized by fixed-tissue approaches. By focusing on healthy adult mice, incorporating independent biological replicates, profiling lumbar segments, and resolving transcriptional and signaling programs at 6 h and 24 h after an acute bout of exercise in addition to baseline adaptations, our design captures aspects of exercise-induced plasticity that previous atlases could not detect, e.g. temporal waves of gene expression, inter-individual variability, and region-specific remodeling in lumbar segments innervating hind limbs. These distinctions are reflected in the temporal model presented in Figure 7B, which illustrates how prior training alters the amplitude and resolution of transcriptional responses across major glial populations. Together, these observations position our study as a complementary contribution to the emerging literature, defining the cellular and temporal logic of exercise-induced spinal cord adaptation in a physiological, non-aging context.

Limitations and future directions

Several limitations should be noted. First, cholinergic neurons were underrepresented in the snRNA-seq dataset, limiting our ability to detect subtle transcriptional changes. Enriching for cholinergic neurons using ChAT-Cre lines, as applied in prior work (41, 44, 65), will be important for future studies. Second, our timepoints (6h, 24h) may have missed very early (0– 4h) or late (>24h) neuronal responses to an acute bout of exercise. Given the timing of peak firing and neuromuscular junction remodeling, expanding the temporal resolution will be critical. Third, although bulk proteomics validated broad proteome-wide changes, cell-type– resolved proteomics or spatial transcriptomics would strengthen the mechanistic links between transcription, translation, and local circuit remodeling. Then, our study was only performed in male mice, using specific endurance exercise and training paradigms. Sex- dependent differences, as well as the effect of other training protocols, for example resistance exercise, will have to be assessed. Finally, while our study documented extensive molecular adaptation, functional validations and assays of spinal circuit excitability or NMJ remodeling will be essential to establish causal links between glial plasticity and motor performance.

Conclusions

In summary, endurance exercise remodels the healthy spinal cord through coordinated transcriptional plasticity of glia and neuron–glia communication. Training primes oligodendrocytes, astrocytes, and OPCs for rapid metabolic and synaptic responses, accelerates resolution of transcriptional activation, and dampens inflammatory microglial states (Figure 7). Beyond transcription, exercise reshapes communication networks, with cholinergic neurons acting as major signaling hubs that likely drive changes in the transcriptional response of glia. Our findings identify coordinated neuron–glia interactions as fundamental mediators of exercise-induced spinal cord adaptation and provide a framework for understanding how physical activity promotes resilience in neural circuits and benefits motor system function in a healthy state. These findings are important on several levels: first, exercise, in particular when performed to exhaustion, represents a physiological state of prolonged (hyper-) activation of motor systems, most likely contributing to fatigue and therefore exercise cessation (66). Insights on the acute perturbations in untrained and trained mice as well as the persistent chronic adaptations therefore help to understand the molecular and cellular underpinnings of spinal cord plasticity. Second, the pathways and interactions resulting in increased resilience, reduced inflammation, and faster return to baseline in the trained spinal cord, together with findings from other reports describing benefits of physical activity in different contexts of spinal cord pathologies, will help to identify and potentially leverage mechanisms and targets for therapeutic approaches.

Materials and methods

Animals

Male C57BL/6J mice (from Janvier Labs) were maintained in the animal facility of the Biozentrum, University of Basel under standard housing conditions (22 °C; 12-h light/dark cycle) with unrestricted access to water and standard chow. Experimental protocols were approved by the Cantonal Veterinary Office Basel-Stadt and complied with Swiss animal welfare regulations. Animals were habituated to the housing environment for at least one week prior to the start of experiments. At 11 weeks of age, mice were randomly divided into sedentary controls or a training group. The training group was given access to running wheels equipped with sensors to continuously log distance, while sedentary mice were provided with plastic housing structures for enrichment. Wheel running was voluntary and monitored for 6 weeks. Body composition was assessed at baseline and 4 weeks after the start of the running wheel period using an EchoMRI-100 system in awake, briefly restrained animals. Body weight was monitored weekly until the end of the experimental duration.

Exhaustive treadmill test

Endurance capacity was assessed on a motorized treadmill set to a 5° incline at both the beginning and end of the 6-week training phase. The final test was used both to measure adaptation and to provide an acute bout of exercise. Running wheels were removed 48 h before treadmill testing to limit immediate effects of voluntary activity. Mice were acclimatized to the treadmill over two days: day 1 consisted of 5 min each at 0, 5, and 8 m/min, and day 2 of 5 min each at 5, 8, and 12 m/min. On day 3, mice performed the maximal test protocol, starting with 5 min at 5 m/min, followed by 5 min at 8 m/min, and then incremental increases of 2 m/min every 3 min up to 32 m/min. Beyond that, speed was increased by 2 m/min every 10 min until exhaustion. Exhaustion was defined as the inability to resume running despite gentle prodding and brief, mild electrical stimulation (0.5 mA, 200 ms pulses at 1 Hz). Animals were sacrificed either 6 h or 12 h following exhaustion. Mice assigned to the baseline sedentary and trained groups were not exposed to the treadmill to avoid introducing acute transcriptional changes.

Proteomics sample preparation and analysis

Lumbar spinal cord sections of 0.5 cm were finely minced, suspended in lysis buffer (5% SDS, 10mM TCEP, 0.1 M TEAB) and lysed by sonication using a PIXUL Multi-Sample Sonicator (Active Motif) with Pulse set to 50, PRF to 1, Process Time to 10 min and Burst Rate to 20 Hz. Lysates were incubated for 10 min at 95°C, alkylated in 20 mM iodoacetamide for 30 min at 25°C and proteins digested using S-Trap™ micro spin columns (Protifi) according to the manufacturer’s instructions. Shortly, 12 % phosphoric acid was added to each sample (final concentration of phosphoric acid 1.2%) followed by the addition of S-trap buffer (90% methanol, 100 mM TEAB pH 7.1) at a ratio of 6:1. Samples were mixed by vortexing and loaded onto S-trap columns by centrifugation at 4000 g for 1 min followed by three washes with S-trap buffer. Digestion buffer (50 mM TEAB pH 8.0) containing sequencing-grade modified trypsin (1/25, w/w; Promega, Madison, Wisconsin) was added to the S-trap column and incubate for 1h at 47 °C. Peptides were eluted by the consecutive addition and collection by centrifugation at 4000 g for 1 min of 40 ul digestion buffer, 40 ul of 0.2% formic acid and finally 35 ul 50% acetonitrile, 0.2% formic acid. Samples were dried under vacuum and stored at −20 °C until further use.

Dried peptides were resuspended in 0.1% aqueous formic acid, 0.02% DDM (n-Dodecyl-B-D- maltoside) and subjected to LC–MS/MS analysis using a timsTOF Ultra Mass Spectrometer (Bruker) equipped with a CaptiveSpray nano-electrospray ion source (Bruker) and fitted with a Vanquish Neo (Thermo Fisher Scientific). Peptides were resolved using a RP-HPLC column (100 um × 30 cm) packed in-house with C18 resin (ReproSil Saphir 100 C18, 1.5 um resin; Dr. Maisch GmbH) at a flow rate of 0.4 ul/min and column heater set to 60°C. The following gradient was used for peptide separation: from 2% B to 25% B over 25 min to 35% B over 5 min to 95% B over 0.5 min followed by 1 min at 95% B to 2% B over 0.1 min followed by 5 min at 2% B. Buffer A was 0.1% formic acid in water and buffer B was 80% acetonitrile, 0.1% formic acid in water.

The mass spectrometer was operated in dia-PASEF mode with a cycle time estimate of 0.95 s. MS1 and MS2 scans were acquired over a mass range from 100 to 1700 m/z. A method with 8 dia-PASEF scans separated into 3 ion mobility windows per scan covering a 400-1000 m/z range with 25 Da windows and an ion mobility range from 0.64 to 1.37 Vs cm² was used. Accumulation and ramp time were set to 100 ms, capillary voltage was set to 1600V, dry gas was set to 3 l/min and dry temperature was set to 200 °C. The collision energy was ramped linearly as a function of ion mobility from 59 eV at 1/K0 = 1.6 V s cm^-2 to 20 eV at 1/K0 = 0.6 V s cm^-2.

The acquired files were searched using the Spectronaut (Biognosys v19.0) directDIA workflow against a Mus musculus database (consisting of 17085 protein sequences downloaded from Uniprot on 20220222) and 393 commonly observed contaminants. Quantitative fragment ion data (F.Area) was exported from Spectronaut and analyzed using the MSstats R package v.4.13.0. Data was normalised using the default normalisation option “equalizedMedians”, imputed using “AFT model-based imputation” and p-values and q-values for pairwise comparisons were calculated using the limma package.

Nuclei isolation and FACS

Nuclei were isolated from frozen lumbar spinal cord using a protocol adapted from previous published studies(65, 67). Tissue was homogenized in 1 mL of low-sucrose lysis buffer (320 mM sucrose, 10 mM HEPES pH 8.0, 5 mM CaCl₂, 3 mM magnesium acetate, 0.1 mM EDTA, 1 mM DTT, 0.1% NP-40) with a glass Dounce homogenizer (20 strokes with pestle A followed by 20 strokes with pestle B). Homogenates were incubated on ice for 10 min, filtered through a 40 μm cell strainer, and rinsed with 500 μL lysis buffer. The filtrate was collected in low-binding tubes, briefly spun to recover residual volume, and centrifuged at 500 × g for 5 min at 4 °C. Pellets were resuspended in low-sucrose buffer, spun a second time, and washed in 500 μL of wash buffer (PBS with 1% BSA and 0.2 U/μL SUPERaseIn). For enrichment of intact nuclei, 900 μL of 1.5 M sucrose solution was mixed with 500 μL of the nuclei suspension and layered over 500 μL of 1.5 M sucrose buffer in a low-binding tube. Gradients were centrifuged at 3,000 × g for 40 min at 4 °C (swinging-bucket rotor). Pellets were resuspended in wash buffer and centrifuged for 10 min at 500 × g to remove residual sucrose. The resulting nuclei pellet was fixed for 1 h at room temperature using the 10X Genomics Fixation Kit (#1000497) and quenched according to the manufacturer’s protocol. Nuclei were stained with DAPI and sorted on a BD Aria Fusion cytometer to enrich for DAPI-positive events. A small fraction of sorted nuclei was reserved for imaging quality control, while the remainder was cryopreserved following the 10X Genomics Flex protocol (CG000478 | Rev C) and stored at −80 °C until library preparation.

Single-nucleus library preparation

Single-nucleus RNA-seq was performed with the Chromium Next GEM Fixed RNA kit (10X Genomics; User Guide Rev D, CG000522). Nuclei were counted using Acridine Orange/Propidium Iodide staining on a LUNA-FL fluorescence counter. After probe hybridization, samples were pooled to ensure equal representation (pool-and-wash workflow). Each run included 16 multiplexed samples, with a minimum of 80’000 nuclei loaded per sample. Library quality was confirmed using a Bioanalyzer, and sequencing was carried out on an Illumina NovaSeq6000 platform targeting approximately 15’000 reads per nucleus.

Single-nucleus RNA-seq data processing and analysis Preprocessing, quality control, data integration, and clustering

Feature-barcode matrices were generated with the Cell Ranger pipeline (10X Genomics, v7.1.0). Ambient RNA contamination was removed using CellBender (v0.3.0)(68). Count matrices were imported into Seurat v5 (69) for downstream analysis. Nuclei with fewer than 500 or more than 8,000 detected genes, >25,000 UMI counts, or mitochondrial read content above 1% were excluded. Putative doublets were identified using DoubletFinder, and clusters enriched for doublet scores were removed. Data were normalized with SCTransform, regressing out sequencing depth and mitochondrial percentage. Integration across replicates and experimental conditions was performed using Seurat’s IntegrateLayers with batch correction via Harmony. Dimensionality reduction was conducted with PCA followed by UMAP, and clusters were identified using the Louvain algorithm at a resolution of 0.5.

Cell-types were annotated based on canonical marker expression for excitatory neurons (Slc17a6, Slc17a8), inhibitory neurons (Gad1, Gad2, Slc6a5), cholinergic neurons (Chat, Slc18a3, Slc5a7), astrocytes (Gfap, Aldh1l1, Aqp4), oligodendrocytes (Mbp, Plp1, Mag), oligodendrocyte precursor cells (Pdgfra, Cspg4), microglia (P2ry12, Tmem119), vascular cells (Pecam1, Flt1), fibroblasts (Pdgfrb, Vtn), and pericytes (Rgs5). Contaminant clusters were excluded, and closely related clusters were merged into major cell classes.

Pseudobulk differential expression analysis

For each biological replicate, raw UMI counts were aggregated by cell-type and condition to construct pseudobulk count matrices. DESeq2 (70) was used for differential expression testing across training-state and time course comparisons:

• Training effect: TRAIN vs SED

• Acute sedentary responses: SED6H vs SED, SED24H vs SED

• Acute trained responses: TRAIN6H vs TRAIN, TRAIN24H vs TRAIN.

Lowly expressed genes (row sum < 10 across all samples) were excluded. DESeq2 models included condition and time as factors (∼ condition + time + condition:time). Genes were considered significantly regulated if they had an adjusted p-value < 0.01, absolute log₂ fold change ≥ 0.5, and baseMean ≥ 10. Log₂ fold-change shrinkage was applied using the “ashr” method to stabilize estimates.

Gene set enrichment analysis

Ranked DEG lists were analyzed with clusterProfiler using MSigDB gene sets (Hallmark, Reactome, WikiPathways, and GO: Biological Process). Enrichment was performed separately for up- and downregulated genes, and results were visualized as dot plots, heatmaps, and UpSet plots to highlight pathway overlap across cell-types and conditions. Complementary ontology-based enrichment was conducted using the Metascape platform (metascape.org)(71) with default parameters. Where appropriate, custom background gene lists were supplied to account for cell-type–specific expression repertoires. Metascape clusters enriched terms based on shared gene membership and hierarchical clustering, and representative terms from each cluster were selected for visualization.

Line plot visualization of expression dynamics

To visualize temporal expression changes of selected genes, pseudobulk counts were aggregated per replicate, condition, and timepoint (PRE, 6 h, 24 h). Expression values were converted to counts per million (CPM), log-transformed, and summarized as the mean ± standard error across replicates. Line plots were generated with ggplot2, displaying mean trajectories per gene and condition with shaded ribbons representing standard error. Genes were shown in faceted panels with free y-axis scaling, and consistent color schemes and formatting were applied to produce publication-quality figures.

Cell–cell communication analysis

Intercellular communication analysis was performed using the CellChat R package (version 2.1.2)(49), following the standard comparative workflow. Processed Seurat objects were subset by condition and time point (SED, SED6H, SED24H, TRAIN, TRAIN6H, TRAIN24H) and converted into CellChat objects using the createCellChat function, grouping cells by curated major clusters. The CellChat mouse ligand–receptor database was applied. For each subset, we identified overexpressed genes and interactions, inferred communication probabilities, filtered out low-support interactions (<10 cells per partner), and aggregated pathway-level networks.

To enable temporal and training-state comparisons, CellChat objects were merged within sedentary and trained groups and across matched time points. Comparative analyses included quantification of the number and strength of interactions, ranking of pathways by contribution, and visualization of signaling roles using bubble plots, chord diagrams, and role heatmaps. Pathway similarity networks were computed using both functional and structural similarity, followed by UMAP embedding and network clustering to identify groups of related signaling pathways. Network centrality analysis was applied to determine dominant sending and receiving populations across conditions. Direct pairwise comparisons between matched sedentary and trained time points (e.g., SED24H vs TRAIN24H) were also performed to highlight training-dependent differences in signaling patterns.

Statistics

Differential gene expression was performed using DESeq2 in a pseudobulk framework, with Wald tests for pairwise comparisons and Benjamini–Hochberg correction for multiple testing. Genes with adjusted p-value < 0.01, baseMean ≥ 10, and absolute log₂ fold change ≥ 0.5 were considered significantly regulated. For interaction terms, DESeq2’s likelihood ratio test (LRT) was used to assess the effect of training state on time-dependent responses.

For non-sequencing data (e.g., body composition, behavioral readouts), comparisons between two groups were performed using unpaired Student’s t-tests when normality was confirmed (Shapiro–Wilk test) or Mann–Whitney U tests otherwise. Multi-factorial comparisons (e.g., endurance performance, cell-type proportions, regulon activity) were analyzed using two-way ANOVA followed by Tukey’s post hoc test. Unless otherwise specified, n = 4 mice per group.

Supplementary figures

Supplementary Figure S1.

Supplementary Figure S2.

Supplementary Figure S3.

Supplementary Figure S4.

Supplementary Figure S5.

Supplementary Figure S6.

Supplementary Figure S7.

Data availability

Proteomic data have been deposited to the ProteomeXchange Consortium (https://www.proteomexchange.org/) via the MassIVE partner repository (https://massive.ucsd.edu/) with MassIVE data set identifier MSV000099207 and ProteomeXchange identifier PXD068447. The transcriptomics data have been deposited to the Gene Expression Omnibus (GEO) and are available under the accession number GSE313063.

Acknowledgements

We thank the FACS Core Facility of the Biozentrum, especially Stella Stefanova, for assistance with flow sorting and all animal caretakers from the animal facility at the Biozentrum for their help. We thank the Quantitative Genomics Facility of the University of Basel, in particular Phillippe Demougin, for preparing cDNA libraries. Calculations were performed at sciCORE, the scientific computing center at University of Basel, with support by the SIB - Swiss Institute of Bioinformatics.

Additional information

Author Contributions

Conceptualization: S.M. and C.H. Methodology: S.M., S.D., D.R, and C.H. Investigation: S.M., S.D., G.S., and S.A.S. Formal analysis and visualization: S.M. and D.R. Writing: S.M. and C.H. Funding acquisition: C.H. All authors have read and approved the final version of this manuscript and agree to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. All persons designated as authors qualify for authorship, and all those who qualify for authorship are listed.

Funding

This study was supported by the Swiss National Science Foundation (Grant no. CRSII5_209252), the Biozentrum and the University of Basel.

Declaration of generative AI and AI-assisted technologies in the writing process

During the preparation of this manuscript, the authors used ChatGPT (o4, 5) and Perplexity to assist with grammar checking, sentence structure optimization, and improving clarity. These tools were used to refine the language and presentation of ideas. However, the intellectual content, data interpretation, and conclusions remain entirely the authors’ own. The final manuscript was thoroughly reviewed and edited by the authors, who take full responsibility for its content.

Funding

Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (SNF) (CRSII5_209252)

• Christoph Handschin