A single-cell transcriptomic atlas of exercise-induced anti-inflammatory and geroprotective effects across the body

Abstract

Exercise benefits the whole organism, yet, how tissues across the body orchestrally respond to exercise remains enigmatic. Here, in young and old mice, with or without exercise, and exposed to infectious injury, we characterized the phenotypic and molecular adaptations to a 12-month exercise across 14 tissues/organs at single-cell resolution. Overall, exercise protects tissues from infectious injury, although more effectively in young animals, and benefits aged individuals in terms of inflammaging suppression and tissue rejuvenation, with structural improvement in the central nervous system and systemic vasculature being the most prominent. In vascular endothelial cells, we found that readjusting the rhythmic machinery via the core circadian clock protein BMAL1 delayed senescence and facilitated recovery from infectious damage, recapitulating the beneficial effects of exercise. Our study underscores the effect of exercise in reconstituting the youthful circadian clock network and provides a foundation for further investigating the interplay between exercise, aging, and immune challenges across the whole organism.

Keywords: BMAL1; aging; circadian clock; exercise; inflammation; single-cell RNA sequencing.

Graphical abstract

Figure 1

Construction of single-cell transcriptomic atlas of exercise-induced effects across multiple tissues in mice (A) Schematic diagram of mouse exercise and multi-dimensional analysis of its systematic effects. Y-Ctrl, young control; Y-Ex, young exercise; O-Ctrl, old control; O-Ex, old exercise. (B) Boxplot showing the average running distance (km per day) of mice in the Y-Ex and O-Ex groups every 3 months. n = 22–24 mice per group. (C) Body weights (top) and physical score (bottom) of mice from different groups as indicated. n = 22–24 mice per group. (D) Bar plots showing the mean time on a rotarod (top left), grip strength of four limbs (top right), distance on the treadmill (bottom left), and the percentage of correct alternations by Y-maze test (bottom right) of mice from indicated groups. n = 15–24 mice per group. (E) ELISA analysis of the concentration of IL-1β (top) and the concertation ratio of AST to ALT (bottom) in the serum of indicated groups. n = 5–12 mice per group. (F) t-SNE plots showing different cell types across the five tissues based on data from snRNA-seq (left) or nine tissues by scRNA-seq (right). Oligo, oligodendrocyte; Ast, astrocyte; ExN, excitatory neuron; InN, inhibitory neuron; Granule, granule cell; MLI, molecular layer interneuron; Epe, ependymal cell; OPC, oligodendrocyte precursor cell; Mic, microglia; Men, meningeal cell; Tendon, tendon fibroblast; Fast IIX, type IIX fast-twitch fiber; Fast IIA, type IIA fast-twitch fiber; Fast IIB, type IIB fast-twitch fiber; NMJ_post, postsynaptic muscle fiber of neuromuscular junction; Adi, adipocyte; MTJ, myotendinous junction; MuSC, muscle stem cell; Fib, fibroblast; Car, cardiocyte; Per, pericyte; EC, endothelial cell; SMC, smooth muscle cell; TC, T cell; M1, type I macrophage; M2, type II macrophage; AT1, type I alveolar epithelial cell; AT2, type II alveolar epithelial cell; Epi, epithelial cell; AM, alveolar macrophage; CD-IC, collecting duct intercalated cell; CD-PC, collecting duct principal cell; CD-Trans, collecting duct transitional cell; DCT, distal convoluted tubule; LOH, ascending loop of Henle; DLOH, descending loop of Henle; PT-S1, segment1 of proximal tubule; PT-S2, segment2 of proximal tubule; PT-S3, segment3 of proximal tubule; Hep, hepatocyte; Cho, cholangiocyte; Kup, Kupffer cell; NKT, natural killer T cell; SPG, spermatogonia; RS, round spermatid; ES, elongated spermatid; SPC, spermatocyte; Neu, neutrophil; Mono, monocyte; Bas, basophil; LP-B, late pro-B cell; Mega, megakaryocyte; mDC, myeloid dendritic cell; pDC, plasmacytoid dendritic cell; CD4_Naive, CD4⁺ naive T cell; CD4_Mem, CD4⁺ memory T cell; CD8_Naive, CD8⁺ naive T cell; CD8_Mem, CD8⁺ memory T cell; CD8_CTL, CD8⁺ cytotoxic T cell; BC, B cell; Pla, plasmocyte. The quantification data in (C)–(E) are shown as the means ± SEM, and two-tailed Student's t test p values are indicated.

Figure 2

Exercise reprograms transcriptional profiles across tissues in an age-dependent manner (A) Venn diagram showing the YE-specific DEGs, OE-specific DEGs, as well as common DEGs overlapping between YE DEGs (Y-Ex/Y-Ctrl) and OE DEGs (O-Ex/O-ctrl) based on data from scRNA-seq and snRNA-seq. The DEG numbers were calculated by summing the cell type-specific DEGs from the 14 tissues or organs, yielding the total number of upregulated and downregulated genes. (B) Bar plots showing the DEG number (left) or DEG relative ratio (corresponding DEG number to total number of DEGs including YE-specific DEGs, common DEGs, and OE-specific DEGs, right) of indicated DEGs across different tissues based on data from scRNA-seq and snRNA-seq. (C) Bubble plot showing the difference in the DEG ratio of YE DEG number or OE DEG number to total DEG number in different cell types of the 14 tissues based on data from scRNA-seq and snRNA-seq. Cell types that have a DEG ratio difference of more than 60% are named YE-specific cell types or OE-specific cell types. Cell types that have a DEG ratio difference lower than 20% are named Common cell types. (D) Network plots showing the enriched GO terms/pathways of upregulated DEGs during exercise in the indicated groups (left) or indicated tissues (right) based on data from scRNA-seq and snRNA-seq. (E) WGA staining of skeletal muscle tissues from the indicated groups. The relative cross-sectional areas are quantified as fold changes. Scale bars, 100 μm. n = 7–12 mice per group. (F) Electron microscopy analysis of the number of intermyofibrillar mitochondria in skeletal muscle (top) and the percentage of total interfibrillar mitochondrial area in heart tissues (bottom) from the indicated groups. Scale bars, 1 μm. n = 4 mice per group and 10 images were taken from each mouse. (G) Immunofluorescence analysis of ChAT-positive motor neurons (MNs) in the motor cortex and spinal cord from the indicated groups. Scale bars, 50 μm. n = 7–12 mice per group. (H) Immunofluorescence analysis of neurofilament (NF)-positive terminal buttons in skeletal muscle tissues from the indicated groups. The number of terminal buttons per 0.6 mm² is measured across different groups. Scale bars, 50 μm. n = 7–12 mice per group. (I) Immunochemistry analysis of NeuN-positive neurons in the brain cortex from the indicated groups. The thickness of the cortex from the indicated regions was measured across different groups. Scale bars, 50 μm. n = 7–12 mice per group. 1# and 2# cortexes are the motor and sensory cortex above the corpus callosum (CC) and dentate gyrus (DG) regions, respectively. The quantification data in (F)–(I) are shown as the means ± SEM, and two-tailed Student's t test p values are indicated.

Figure 3

Exercise safeguards young mice against acute inflammatory damage (A) Schematic diagram of mice treated with vehicle or LPS after a 12-month exercise period and subsequent analysis of five tissues as indicated. (B) Venn diagrams showing the overlapped genes between Y-LPS (Y-LPS/Y-Ctrl) and YE-LPS (Y-Ex-LPS/Y-LPS) (top), as well as between O-LPS (O-LPS/Y-Ctrl) and OE-LPS (O-Ex-LPS/O-LPS) (bottom), based on data from scRNA-seq and snRNA-seq. YE-rev-LPS-DEGs or OE-rev-LPS-DEGs were defined as subsets of overlapped DEGs that were changed in the opposite direction in young control and young exercise samples, or old control and old exercise samples, respectively. (C) Bar plots showing the DEG number (top) and DEG ratio (bottom) of LPS-specific DEGs and Ex-rev-LPS DEGs to total LPS DEGs (Y-LPS and O-LPS) across different tissues in both young and aged mice based on data from scRNA-seq and snRNA-seq (left) or bulk RNA-seq (right). LPS DEGs excluding EX-rev-LPS DEGs are defined as LPS-specific DEGs. (D) Bar plots showing the DEG number (top) and DEG ratio (bottom) of LPS DEGs and Ex-rev-LPS DEGs to total LPS DEGs in different cell types across five tissues in both young and aged mice based on data from scRNA-seq and snRNA-seq. (E) Heatmaps showing the expression profile of genes in response to LPS across different tissues in indicated groups based on data from scRNA-seq and snRNA-seq. (F) Network plots showing enriched GO terms/pathways for the Rev-LPS DEGs downregulated by exercise in the indicated comparisons (top) or indicated tissues (bottom). (G) Immunofluorescence analysis of CD45-positive cells (immune cells) in liver (left) and lung (right) tissues from the indicated groups. The percentages of CD45-positive cells were quantified. Scale bars, 50 μm. n = 7–12 mice per group. (H) Immunofluorescence analysis of neutrophils (anti-granulocytes antibody, HIS48) in liver (left) and lung (right) tissues from the indicated groups. The percentage of neutrophils was quantified. Scale bars, 50 μm. n = 7–12 mice per group. (I) Immunochemistry analysis of F4/80-positive cells in liver (left) and lung (right) tissues from the indicated groups. The percentages of F4/80-positive cells were quantified. Scale bars, 50 μm. n = 7–12 mice per group. (J) TUNEL staining in liver tissues from the indicated groups. The percentages of TUNEL-positive cells were quantified. Scale bars, 50 μm. n = 7–12 mice per group. (K) Immunochemistry analysis of IL-1β-positive cells in liver tissues from the indicated groups. The percentages of IL-1β-positive cells were quantified. Scale bars, 50 μm. n = 7–12 mice per group. The quantification data in (G)–(K) are shown as the means ± SEM, and two-tailed Student's t test p values are indicated.

Figure 4

Exercise reverses aging-related gene expression across tissues (A) Venn diagrams showing the overlapped DEGs between Aging DEGs (O-Ctrl/Y-Ctrl) and OE DEGs (O-Ex/O-Ctrl) based on data from scRNA-seq and snRNA-seq. Pro-aging DEGs were the subset of OE DEGs changed in the same direction as Aging DEGs, and Rev-aging DEGs were the subset of OE DEGs changed in the opposite direction as Aging DEGs. The DEG numbers were calculated by summing the cell type-specific DEGs from the 14 tissues or organs, yielding the total number of upregulated and downregulated genes. (B) Bar plots showing the DEG ratio (left) and DEG number (right) of Rev-aging, Pro-aging DEGs, and Aging-specific DEGs across different tissues based on data from scRNA-seq and snRNA-seq. The DEG ratio was calculated by comparing Rev-aging, Pro-aging DEGs, and Aging-specific DEGs with total Aging DEGs. Aging DEGs excluding Rev-aging and Pro-aging DEGs are defined as Aging-specific DEGs. (C) Bar plots showing the DEG number (top) and DEG ratio (bottom) of Rev-aging, Pro-aging DEGs, and Aging-specific DEGs to total Aging DEGs across different cell types based on data from scRNA-seq and snRNA-seq. (D) Dot plot showing the ratios of rescued DEGs (comparing upregulated or downregulated Rev-aging DEGs with total aging-DEGs) across all cell types of different tissues. The size of dots indicates the number of Rev-aging DEGs. (E) Bar plots showing the DEG ratio of Rev-aging DEGs, Pro-aging DEGs, and Aging-specific DEGs to total Aging DEGs in endothelial cell and seven T cell subpopulations as indicated across different tissues. (F) Heatmaps showing the enriched GO terms/pathways for upregulated (top) and downregulated (down) Rev-aging DEGs by exercise across different tissues in aged mice based on data from scRNA-seq and snRNA-seq. (G) Heatmaps showing the changes in ligand-receptor interactions between different cell types across all the tissues in the Aging (O-Ctrl/Y-Ctrl) and OE (O-Ex/O-Ctrl) comparison groups based on data from scRNA-seq and snRNA-seq. (H) Bubble plots showing the enriched GO terms/pathways for exercise-erased (left) or exercise-rescued (right) ligand-receptor interactions across all the cells in the indicated tissues.

Figure 5

Exercise alleviates a panel of aging-associated phenotypes across tissues (A) Immunofluorescence analysis of IBA1-labeled microglial cells in the cortex and hippocampal DG region, cerebellum, spinal gray matter (GM), and white matter (WM) from the indicated groups. Scale bars, 50 μm. The percentages of IBA1-positive microglia were quantified. n = 7–12 mice per group. (B) Immunostaining of GFAP in the hippocampus, cerebellum, spinal gray matter (GM), and white matter (WM) from the indicated groups. Scale bars, 50 μm. The percentages of GFAP-positive cells were quantified. n = 7–12 mice per group. (C) Immunostaining of CD45 in lung, liver, heart, kidney, skeletal muscle, and small intestine tissues from the indicated groups. Scale bars, 50 μm. The percentages of CD45-positive immune cells were quantified. n = 5–12 mice per group. (D) Immunofluorescence analysis of neutrophils in liver, lung, and kidney tissues from the indicated groups. Scale bars, 50 μm. The percentages of neutrophils were quantified. n = 5–12 mice per group. (E) Immunochemistry analysis of IL-1β-positive areas in liver and kidney tissues from the indicated groups. The percentages of IL-1β-positive cells were quantified. Scale bars, 50 μm. n = 7–12 mice per group. The quantification data in (A)–(E) are shown as the means ± SEM, and two-tailed Student's t test p values are indicated.

Figure 6

Exercise resets the circadian clock machinery across aged tissues (A) Network plots showing the upregulated aging TFs (left) and downregulated aging TFs (right). (B) Heatmaps showing the expression profiles of overlapped Rev-aging DEGs based on the sc/snRNA-seq and bulk RNA-seq data across different cell types from different tissues. The color key from blue to yellow indicates gene expression levels from low to high. (C) Bubble plots showing the cell type frequency and tissue frequency (more than 5 tissues and 20 cell types) of upregulated (left) and downregulated (right) rescued DEGs. (D) Heatmaps showing the expression profiles of circadian rhythm-related Rev-aging DEGs for different cell types across different tissues. The color key from blue to red indicates gene expression levels from low to high.

Figure 7

BMAL1 plays a protective role in endothelial cells (A) Bar plot showing the expression changes of Arntl (Bmal1) in different cell types from all the tissues analyzed in the indicated groups. (B) Bar plots showing the rhythm index of different cell types from all the tissues analyzed in the indicated groups. (C) Transverse sections of heart tissues from the indicated groups of mice were subjected to RNA in situ hybridization (RNA-ISH) with Bmal1 riboprobes. Left, representative ISH images; right, the relative expression level of Bmal1 in mouse hearts was quantified. Scale bars, 50 μm. n = 7 mice per group. (D) Western blotting of BMAL1 protein levels in CAECs at the indicated passages. n = 3 biological repeats per group. Young, passage 4; Senescent, passage 11. RS, replicative senescence. (E) Western blotting of BMAL1 protein levels in CAECs transduced with lentivirus with non-targeting or Bmal1-targeting sgRNAs in the CRISPR-Cas9-mediated knockout system. (F) SA-β-gal analysis of CAECs transduced with lentivirus with non-targeting or Bmal1-targeting sgRNAs in the CRISPR-Cas9-mediated knockout system. Scale bars, 100 μm. n = 4 biological repeats per group. (G) Transcript levels of indicated genes in CAECs transduced with lentivirus with non-targeting or Bmal1-targeting sgRNAs in the CRISPR-Cas9-mediated knockout system. n = 3 biological repeats per group. (H) Monocyte adhesion assay in CAECs transduced with lentivirus with non-targeting or Bmal1-targeting sgRNAs in the CRISPR-Cas9-mediated knockout system. Scale bars, 50 μm. n = 3 biological repeats per group. (I) Western blotting of BMAL1 protein levels in CAECs transduced with lentivirus-expressing Flag-tagged GAL4 or BMAL1. n = 3 biological repeats per group. (J) SA-β-gal analysis of CAECs transduced with lentivirus-expressing Flag-tagged GAL4 or BMAL1. The percentage of SA-β-gal-positive cells was quantified. Scale bars, 100 μm. n = 3 biological repeats per group. (K) Transcript levels of the indicated genes associated with endothelial cell dysfunction in CAECs transduced with lentivirus-expressing Flag-tagged GAL4 or BMAL1. n = 3 biological repeats per group. (L) ELISA analysis of the secretion of IL-6 in CAECs transduced with lentivirus-expressing Flag-tagged GAL4 or BMAL1. The absorption value is normalized to the cell number. n = 3 biological repeats per group. (M) Immunofluorescence staining of Ki67 in CAECs transduced with lentivirus-expressing Flag-tagged GAL4 or BMAL1. Scale bars, 50 μm. n = 4 biological repeats per group. (N) Assessment of tube formation capability in CAECs transduced with lentivirus-expressing Flag-tagged GAL4 or BMAL1. The relative branch points (right upper) and total length of tubes (right lower) were quantified as fold changes. Scale bar, 100 μm. n = 3 biological repeats per group. (O) Wound scratch assay showing the migration ability of CAECs transduced with lentivirus-expressing Flag-tagged GAL4 or BMAL1. Red lines represent scar boundaries. The relative migration ability was quantified as the fold change. Scale bar, 100 μm. n = 3 biological repeats per group. (P) Immunochemistry analysis of CD31-positive (left) and VCAM-1-positive (right) cells in the heart using adjacent tissue sections from the indicated groups. The percentages of CD31- or VCAM-1-positive areas were quantified. Scale bars, 50 μm. n = 5–12 mice per group. The quantification data in (C)–(P) are shown as the means ± SEM, and two-tailed Student's t test p values are indicated.