The ability to slow or reverse biological ageing would have major implications for mitigating disease risk and maintaining vitality1. Although an increasing number of interventions show promise for rejuvenation2, their effectiveness on disparate cell types across the body and the molecular pathways susceptible to rejuvenation remain largely unexplored. Here we performed single-cell RNA sequencing on 20 organs to reveal cell-type-specific responses to young and aged blood in heterochronic parabiosis. Adipose mesenchymal stromal cells, haematopoietic stem cells and hepatocytes are among those cell types that are especially responsive. On the pathway level, young blood invokes new gene sets in addition to reversing established ageing patterns, with the global rescue of genes encoding electron transport chain subunits pinpointing a prominent role of mitochondrial function in parabiosis-mediated rejuvenation. We observed an almost universal loss of gene expression with age that is largely mimicked by parabiosis: aged blood reduces global gene expression, and young blood restores it in select cell types. Together, these data lay the groundwork for a systemic understanding of the interplay between blood-borne factors and cellular integrity.
This is a preview of subscription content
Subscription info for Chinese customers
We have a dedicated website for our Chinese customers. Please go to naturechina.com to subscribe to this journal.
Get time limited or full article access on ReadCube.
All prices are NET prices.
Code used for downstream analyses of the raw count matrices is available at GitHub:https://github.com/rpalovics/parabiosis/releases/tag/0.1.0 (https://doi.org/10.5281/zenodo.5762659).
The dataset can be explored interactively at https://ccb-web.cs.uni-saarland.de/parabiosis/. Gene counts and metadata are available on Figshare: https://figshare.com/projects/Molecular_hallmarks_of_heterochronic_parabiosis_at_single_cell_resolution/127628. Raw sequencing data are available in the Gene Expression Omnibus under accession code GSE132042 and via the following public AWS S3 bucket: https://registry.opendata.aws/tabula-muris-senis/. The versions of the noted databases used are available at the following URLs: STRING database, https://version-11-0.string-db.org/; MitoCarta 3.0, https://www.broadinstitute.org/files/shared/metabolism/mitocarta/mouse.mitocarta3.0.html; GENCODE vM19, https://www.gencodegenes.org/mouse/release_M19.html.
López-Otín, C., Blasco, M. A., Partridge, L., Serrano, M. & Kroemer, G. The hallmarks of aging. Cell 153, 1194–1217 (2013).
Mahmoudi, S., Xu, L. & Brunet, A. Turning back time with emerging rejuvenation strategies. Nat. Cell Biol. 21, 32–43 (2019).
Schaum, N. et al. Single-cell transcriptomics of 20 mouse organs creates a Tabula Muris. Nature 562, 367–372 (2018).
Schaum, N. et al. Ageing hallmarks exhibit organ-specific temporal signatures. Nature 583, 596–602 (2020).
Almanzar, N. et al. A single-cell transcriptomic atlas characterizes ageing tissues in the mouse. Nature 583, 590–595 (2020).
Castellano, J. M., Kirby, E. D. & Wyss-Coray, T. Blood-borne revitalization of the aged brain. JAMA Neurol. 72, 1191–1194 (2015).
Villeda, S. A. et al. The ageing systemic milieu negatively regulates neurogenesis and cognitive function. Nature 477, 90–96 (2011).
Katsimpardi, L. et al. Vascular and neurogenic rejuvenation of the aging mouse brain by young systemic factors. Science 344, 630–634 (2014).
Smith, L. K. et al. β2-microglobulin is a systemic pro-aging factor that impairs cognitive function and neurogenesis. Nat. Med. 21, 932–937 (2015).
Khrimian, L. et al. Gpr158 mediates osteocalcin’s regulation of cognition. J. Exp. Med. 214, 2859–2873 (2017).
Conboy, I. M. et al. Rejuvenation of aged progenitor cells by exposure to a young systemic environment. Nature 433, 760–764 (2005).
Ma, S. et al. Caloric restriction reprograms the single-cell transcriptional landscape of Rattus norvegicus aging. Cell 180, 984–1001 (2020).
Das, M. M. et al. Young bone marrow transplantation preserves learning and memory in old mice. Commun. Biol. 2, 73 (2019).
Baht, G. S. et al. Exposure to a youthful circulation rejuvenates bone repair through modulation of β-catenin. Nat. Commun. 6, 7131 (2015).
Kovina, M. V., Zuev, V. A., Kagarlitskiy, G. O. & Khodarovich, Y. M. Effect on lifespan of high yield non-myeloablating transplantation of bone marrow from young to old mice. Front. Genet. https://doi.org/10.3389/fgene.2013.00144 (2013).
Wang, C.-H. et al. Bone marrow rejuvenation accelerates re-endothelialization and attenuates intimal hyperplasia after vascular injury in aging mice. Circ. J. 77, 3045–3053 (2013).
Smith, L. K. et al. The aged hematopoietic system promotes hippocampal-dependent cognitive decline. Aging Cell 19, e13192 (2020).
Lee, S. S. et al. A systematic RNAi screen identifies a critical role for mitochondria in C. elegans longevity. Nat. Genet. 33, 40–48 (2003).
Zhang, M. J., Pisco, A. O., Darmanis, S. & Zou, J. Mouse Aging Cell Atlas analysis reveals global and cell type specific aging signatures revision 1. Elife 10, e62293 (2021).
Kular, J. K., Basu, S. & Sharma, R. I. The extracellular matrix: structure, composition, age-related differences, tools for analysis and applications for tissue engineering. J. Tissue Eng. 5, https://doi.org/10.1177/2041731414557112 (2014).
Pinti, M. et al. Aging of the immune system: focus on inflammation and vaccination. Eur. J. Immunol. 46, 2286–2301 (2016).
Villeda, S. A. et al. Young blood reverses age-related impairments in cognitive function and synaptic plasticity in mice. Nat. Med. 20, 659–663 (2014).
Picelli, S. et al. Full-length RNA-seq from single cells using Smart-seq2. Nat. Protoc. 9, 171–181 (2014).
Darmanis, S. et al. A survey of human brain transcriptome diversity at the single cell level. Proc. Natl Acad. Sci. USA 112, 7285–7290 (2015).
Picelli, S. et al. Tn5 transposase and tagmentation procedures for massively scaled sequencing projects. Genome Res. 24, 2033–2040 (2014).
Hennig, B. P. et al. Large-scale low-cost NGS library preparation using a robust Tn5 purification and tagmentation protocol. G3 (Bethesda) 8, 79–89 (2018).
Luecken, M. D. & Theis, F. J. Current best practices in single‐cell RNA‐seq analysis: a tutorial. Mol. Syst. Biol. 15, e8746 (2019).
Wolf, F. A., Angerer, P. & Theis, F. J. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 19, 15 (2018).
Polański, K. et al. BBKNN: fast batch alignment of single cell transcriptomes. Bioinformatics 36, 964–965 (2020).
Traag, V. A., Waltman, L. & Van Eck, N. J. From Louvain to Leiden: guaranteeing well-connected communities. Sci. Rep. 9, 5233 (2019).
McInnes, L., Healy, J. & Melville, J. UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction. Preprint at https://arxiv.org/abs/1802.03426 (2018) .
Korsunsky, I. et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat. Methods 16, 1289–1296 (2019).
Pedregosa, F. et al. Scikit-learn: machine Learning in Python. J. Machine Learn. Res. 12, 2825–2830 (2011).
Mann, H. B. & Whitney, D. R. On a test of whether one of two random variables is stochastically larger than the other. Ann. Math. Stat. 18, 50–60 (1947).
Benjamini, Y. & Hochberg, Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. B 57, 289–300 (1995).
Enge, M. et al. Single-cell analysis of human pancreas reveals transcriptional signatures of aging and somatic mutation patterns. Cell 171, 321–330 (2017).
Kimmel, J. C. et al. Murine single-cell RNA-seq reveals cell-identity- and tissue-specific trajectories of aging. Genome Res. 29, 2088–2103 (2019).
Gerstner, N. et al. GeneTrail 3: advanced high-throughput enrichment analysis. Nucleic Acids Res. 48, W515–W520 (2020).
Carbon, S. et al. The Gene Ontology Resource: 20 years and still GOing strong. Nucleic Acids Res. 47, D330–D338 (2019).
Rath, S. et al. MitoCarta3.0: an updated mitochondrial proteome now with sub-organelle localization and pathway annotations. Nucleic Acids Res. 49, D1541–D1547 (2021).
Gu, Z., Eils, R. & Schlesner, M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 32, 2847–2849 (2016).
Wickham, H. ggplot2. Wiley Interdiscip. Rev. Comput. Stat. 3, 180–185 (2011).
Hagberg, A., Swart, P. & Schult, D. Exploring Network Structure, Dynamics, and Function Using Networkx (Los Alamos National Laboratory, 2008); https://www.osti.gov/servlets/purl/960616.
Szklarczyk, D. et al. STRING v11: protein–protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 47, D607–D613 (2019).
Ewels, P. A. et al. The nf-core framework for community-curated bioinformatics pipelines. Nat. Biotechnol. 38, 276–278 (2020).
Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013).
Li, B. & Dewey, C. N. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 12, 323 (2011).
Soneson, C., Love, M. I. & Robinson, M. D. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences [version 1; peer review: 2 approved]. F1000Research 4, 1521 (2015).
Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).
Chen, E. Y. et al. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinform. 14, 128 (2013).
We thank the members of the laboratory of T.W.-C. and the Chan Zuckerberg Biohub for feedback and support. Financial support for library preparation, sequencing and AWS time was provided by the Chan Zuckerberg Biohub. Further financial support was provided by the Department of Veterans Affairs (BX004599 to T.W.-C.), the National Institute on Aging (R01-AG045034, AG064928, AG072255, and DP1-AG053015 to T.W.-C.), the NOMIS Foundation (T.W.-C.), The Glenn Foundation for Medical Research (T.W.-C.), Nan Fung Life Sciences (T.W.-C.), and the Wu Tsai Neurosciences Institute and Bertarelli Foundation (T.W.-C.). This work was supported by the National Institute of Aging and the National Institutes of Health under award number P30AG059307.
The authors declare no competing interests.
Peer review information
Nature thanks Johan Auwerx, Yousin Suh, Fan Zhang and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Extended data figures and tables
a, Outline of computational analyses. Single-cell count data are processed per tissue, see Methods ‘Quality control’-’Cell type annotation’. Differential gene expression is then conducted per cell type and comparison (AGE, ACC, REJ) within each tissue, see Methods ‘Differential gene expression’. All of the next panels present data after quality control. b, Number of cells per tissue and replicate. Replicates are colored by their condition. c, Number of replicates per tissue. Replicates are colored by their condition. d, Total number of cells per tissue. e, Fraction of cells within each condition per tissue. f-i, For each experimental condition within each tissue: total read counts (f), the percent of reads mapped to ribosomal genes (g), mitochondrial genes (h), and ERCC spike-ins (i) plotted against the mean number of genes expressed. j, Average LISI scores of mouse replicates calculated over the batch corrected tissue specific UMAP embeddings plotted against the mean LISI scores of tissue specific UMAP embeddings calculated from neighborhood graphs without batch correction. k, Mean entropy batch mixing of mouse replicates calculated over the tissue specific batch-corrected neighborhood graph plotted against the mean entropy batch mixing calculated from neighborhood graphs without batch correction. l, Result of saturation analyses shown per condition (Y, A, IY, HY, IA, HA). Downsampling was carried out per condition within each tissue separately. Results indicate the number of detected genes as the function of the downsampled total counts.
All panels present data after quality control. a, Number of cells per cell type and replicate. Replicates are colored by their condition. b, Number of replicates per cell type. Replicates are colored by their condition. c, Total number of cells per cell type. d, Fraction of cells within each condition per cell type. e-h, For each experimental condition within each cell type, total read counts (e), the percent of reads mapped to ribosomal genes (f), mitochondrial genes (g), and ERCC spike-ins (h) plotted against the mean number of genes expressed.
a, Number of DEGs plotted against the total number of cells within the control and treatment groups. Each dot represents a DGE comparison of a cell type. b-d, Cumulative distributions of the calculated effect size (b), -log10(adj. p-value) (c) and log2 fold change values (d). Distributions are shown separately for ACC, REJ and AGE DGE. Vertical lines indicate the cutoffs applied throughout the study. e, Summary of ACC DGE results. Each cell type that has at least 50 cells in IY and HY is studied in the context of ACC and hence shown. From top to bottom: control and treatment sample sizes indicated separately for AGE and ACC, the number of genes differentially expressed in AGE and ACC, overlaps between AGE and ACC. Overlaps are normalized by the number of DEGs in the union of ACC and AGE DEGs. f, Summary of REJ DGE results. Each cell type that has at least 50 cells in IA and HA is studied in the context of REJ and hence shown. From top to bottom: control and treatment sample sizes indicated separately for AGE and REJ, the number of genes differentially expressed in AGE and REJ, overlaps between AGE and REJ. Overlaps are normalized by the number of DEGs in the union of REJ and AGE DEGs. g, Percent of DEGs that change in the same direction with AGE and ACC are plotted against the total number of DEGs within AGE and ACC for each comparison. Percentages are based on the union of DEGs as defined in (e) h, Percent of DEGs that change in the opposite direction with AGE and REJ are plotted against the total number of DEGs within AGE and REJ for each comparison. Percentages are based on the union of DEGs as defined in (f). i, Fraction of DEGs changing in the same direction with AGE and ACC plotted against the fraction of DEGs changing in the opposite direction with AGE and ACC. Each dot represents a cell type of the study. Colored area indicates where more DEGs change in the same direction than in the opposite direction. j, Fraction of DEGs changing in the same direction with AGE and REJ plotted against the fraction of DEGs changing in the opposite direction with AGE and REJ. Each dot represents a cell type of the study. Colored area indicates where more DEGs change in the opposite direction than in the same direction.
a, Violin plots showing the number of differentially expressed genes as the function of the number of replicates per comparison. The number of replicates are defined as the minimum number of replicates within the control and treatment groups. Results are shown separately for TMS (left) and parabiosis (right). b, Comparison of differential gene expression results with and without subsampling in case of each cell type specific comparison. Spearman correlation values indicate (dis)similarities between: p-values derived from the original and subsampled data (left), effect sizes calculated from the original and subsampled data (middle), and effect sizes calculated on the original data and p-values derived from the subsampled datasets. c, Number of DEGs identified at different p-value and effect size cutoffs per comparison in AGE (left), ACC (middle), and REJ (right).
a, Top list of the 50 most frequent DEGs identified for ACC and REJ. Results are shown separately for up and downregulation. Columns with darker bars indicate top lists where only changes consistent with AGE are shown. These include genes changing in the same direction with ACC and AGE, as well as genes changing in the opposite direction with REJ and AGE. b, DGE results for marrow HSCs for ACC (left) and for REJ (right). from top to bottom: volcano plots (top) show top DEGs. Comparisons of log2-fold changes (middle) show changes with parabiosis on the x-axis and with normal ageing on the y-axis. DEGs with adj. p-value < 0.05, eff. size > 0.6 are shown. Areas where ACC and AGE change in the same direction as well as where REJ and AGE change in the opposite direction are highlighted. Top pathways (GO Biological Process) with highest ‘Combined scores’ defined as in Enrichr50 are shown at the bottom. c, Most enriched pathways (GO Biological Process) among the 100 most frequent DEGs shared across multiple cell types. Results shown for ACC and ACC-AGE same direction (top), and REJ and REJ-AGE opposite direction (bottom). Combined scores are defined as in Enrichr50. d, Gene expression violin plots for liver hepatocytes, GAT MSCs and marrow HSCs of select genes encoding proteins of the electron transport chain. Significance values show the adj. p-values of the Wilcoxon–Mann–Whitney test (two-sided) based differential gene expression, see Methods: ‘Differential gene expression’.
a, Log2-fold changes with AGE, ACC and REJ of genes associated with the 5 OXPHOS complexes. Changes with adj. p-val.<0.05 and eff. size>0.6 are shown. Each column corresponds to one complex and the three separate colors distinguish between AGE, ACC and REJ. b, Spearman correlation of gene expression values with age in case of genes associated with the 5 OXPHOS complexes in the Tabula Muris Senis bulk dataset. Data has been analyzed as in2. Correlation values with adj. p-value < 0.05 are shown. c, Log2-fold changes with ACC and REJ of genes associated with the 5 OXPHOS complexes in the bulk parabiosis dataset.
a–c, Mean number of genes expressed within each cell type, x and y axes indicate Y and A (a), IY and HY (b) and IA and HA (c), each dot represents a cell type. d–f, Cell-cell variability within each cell type in Y and A (d), IY and HY (e) and IA and HA (f), each dot represents a cell type. g–i, Overdispersion within each cell type in Y and A (g), IY and HY (h) and IA and HA (i), each dot represents a cell type.
a, Top 10 most differently affected pathways over all ACC tissues and cell types (top, largest difference at the top) and top 10 most differently affected pathways over all AGE (bottom, largest difference at the bottom) tissues and cell types. b, Same as for a, comparing REJ and AGE pathways. c, Heatmap showing the top 30 most strongly affected pathways in AGE, ACC and REJ in hematopoietic stem cells (HSCs) of the marrow. Pathways related to mitochondria are highlighted in green on the left. Each entry of the heatmap shows the significance level and the number of genes associated with the pathway. d, Same as for c, showing the top 30 most strongly affected pathways for beta cells of the pancreas.
a–c, AGE (a), ACC (b), and REJ (c) DGE based cosine similarity matrices of the cell types studied, see Methods section ‘Ageing and rejuvenation similarity analysis’. All matrices are clustered with complete link hierarchical clustering. d, Force-directed network visualization of the STRING links between DEGs common to MSCs from GAT, MAT, SCAT, bladder, limb muscle, and diaphragm. e, Force-directed network visualization of the STRING links between DEGs common to MSCs (GAT and MAT), hepatocytes, basal and epidermal cells (skin), and HSCs and macrophages (marrow). All links with >0.9 STRING confidence score (scale from 0-1) are queried and shown. f, g Most enriched pathways (GO Biological Process) among the nodes of the networks shown in (d, e), combined scores are defined as in Enrichr50.
Statistics related to quality control, batch correction and cell type annotation for the data used after quality control throughout the study. Results shown separately, for each tissue analysed.
AGE, ACC and REJ DGE results shown for each cell type separately. Blank plots indicate insufficient data (that is, fewer than 50 cells per control or treatment group of the specific comparison).
Heatmaps showing the top 30 most strongly affected pathways in AGE, ACC and REJ for each tissue and cell type. Pathways related to mitochondria are highlighted in green. Heatmaps show the significance level and the number of genes associated with the corresponding pathway.
This file contains the legends for Supplementary Tables 1–8.
This file contains Supplementary Tables 1–8; see separate PDF for table legends.
About this article
Cite this article
Pálovics, R., Keller, A., Schaum, N. et al. Molecular hallmarks of heterochronic parabiosis at single-cell resolution. Nature 603, 309–314 (2022). https://doi.org/10.1038/s41586-022-04461-2