Skip to main content

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

  • Article
  • Published:

Heterogeneity in old fibroblasts is linked to variability in reprogramming and wound healing

Abstract

Age-associated chronic inflammation (inflammageing) is a central hallmark of ageing1, but its influence on specific cells remains largely unknown. Fibroblasts are present in most tissues and contribute to wound healing2,3. They are also the most widely used cell type for reprogramming to induced pluripotent stem (iPS) cells, a process that has implications for regenerative medicine and rejuvenation strategies4. Here we show that fibroblast cultures from old mice secrete inflammatory cytokines and exhibit increased variability in the efficiency of iPS cell reprogramming between mice. Variability between individuals is emerging as a feature of old age5,6,7,8, but the underlying mechanisms remain unknown. To identify drivers of this variability, we performed multi-omics profiling of fibroblast cultures from young and old mice that have different reprogramming efficiencies. This approach revealed that fibroblast cultures from old mice contain ‘activated fibroblasts’ that secrete inflammatory cytokines, and that the proportion of activated fibroblasts in a culture correlates with the reprogramming efficiency of that culture. Experiments in which conditioned medium was swapped between cultures showed that extrinsic factors secreted by activated fibroblasts underlie part of the variability between mice in reprogramming efficiency, and we have identified inflammatory cytokines, including TNF, as key contributors. Notably, old mice also exhibited variability in wound healing rate in vivo. Single-cell RNA-sequencing analysis identified distinct subpopulations of fibroblasts with different cytokine expression and signalling in the wounds of old mice with slow versus fast healing rates. Hence, a shift in fibroblast composition, and the ratio of inflammatory cytokines that they secrete, may drive the variability between mice in reprogramming in vitro and influence wound healing rate in vivo. This variability may reflect distinct stochastic ageing trajectories between individuals, and could help in developing personalized strategies to improve iPS cell generation and wound healing in elderly individuals.

This is a preview of subscription content, access via your institution

Access options

Buy this article

Prices may be subject to local taxes which are calculated during checkout

Fig. 1: Primary fibroblasts from old mice secrete inflammatory cytokines and show increased variability in reprogramming efficiency between mice.
Fig. 2: Old fibroblast cultures exhibit a signature of an inflammatory activated state, which is associated with variability in reprogramming efficiency.
Fig. 3: Age-associated increase in activated fibroblasts and the cytokines that they secrete drive part of the variability in reprogramming between mice.
Fig. 4: Wound healing rate is variable between old mice and correlates with fibroblast subpopulations with distinct cytokine signatures.

Similar content being viewed by others

Data availability

All raw sequencing reads for population RNA-seq, ChIP–seq and single-cell RNA-seq data can be found under BioProject PRJNA316110. The command and configuration files, in addition to a list of all versioned dependencies present in the running environment, are available on the Github repository for this paper (https://github.com/brunetlab/Mahmoudi_et_al_2018) (except for the code for the processing of metabolomics data, which is available upon request).

References

  1. Franceschi, C. & Campisi, J. Chronic inflammation (inflammaging) and its potential contribution to age-associated diseases. J. Gerontol. A Biol. Sci. Med. Sci. 69, S4–S9 (2014).

    Article  PubMed  Google Scholar 

  2. Eming, S. A., Martin, P. & Tomic-Canic, M. Wound repair and regeneration: mechanisms, signaling, and translation. Sci. Transl. Med. 6, 265sr6 (2014).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Lynch, M. D. & Watt, F. M. Fibroblast heterogeneity: implications for human disease. J. Clin. Invest. 128, 26–35 (2018).

    Article  PubMed  PubMed Central  Google Scholar 

  4. Ocampo, A., Reddy, P. & Belmonte, J. C. I. Anti-aging strategies based on cellular reprogramming. Trends Mol. Med. 22, 725–738 (2016).

    Article  PubMed  Google Scholar 

  5. Chaleckis, R., Murakami, I., Takada, J., Kondoh, H. & Yanagida, M. Individual variability in human blood metabolites identifies age-related differences. Proc. Natl Acad. Sci. USA 113, 4252–4259 (2016).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  6. Ong, M. L. & Holbrook, J. D. Novel region discovery method for Infinium 450K DNA methylation data reveals changes associated with aging in muscle and neuronal pathways. Aging Cell 13, 142–155 (2014).

    Article  CAS  PubMed  Google Scholar 

  7. Li, R. et al. Linking inter-individual variability in functional brain connectivity to cognitive ability in elderly individuals. Front. Aging Neurosci. 9, 385 (2017).

    Article  PubMed  PubMed Central  Google Scholar 

  8. Cheung, P. et al. Single-cell chromatin modification profiling reveals increased epigenetic variations with aging. Cell 173, 1385–1397 (2018).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Lapasset, L. et al. Rejuvenating senescent and centenarian human cells by reprogramming through the pluripotent state. Genes Dev. 25, 2248–2253 (2011).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Li, H. et al. The Ink4/Arf locus is a barrier for iPS cell reprogramming. Nature 460, 1136–1139 (2009).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  11. Ravaioli, F., Bacalini, M. G., Franceschi, C. & Garagnani, P. Age-related epigenetic derangement upon reprogramming and differentiation of cells from the elderly. Genes 9, 39 (2018).

    Article  CAS  PubMed Central  Google Scholar 

  12. Mosteiro, L. et al. Tissue damage and senescence provide critical signals for cellular reprogramming in vivo. Science 354, aaf4445 (2016).

    Article  CAS  PubMed  Google Scholar 

  13. Somers, A. et al. Generation of transgene-free lung disease-specific human induced pluripotent stem cells using a single excisable lentiviral stem cell cassette. Stem Cells 28, 1728–1740 (2010).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Guerrero-Juarez, C. F. et al. Single-cell analysis reveals fibroblast heterogeneity and myeloid-derived adipocyte progenitors in murine skin wounds. Nat. Commun. 10, 650 (2019).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  15. Shook, B. A. et al. Myofibroblast proliferation and heterogeneity are supported by macrophages during skin repair. Science 362, eaar2971 (2018).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Frenk, S. & Houseley, J. Gene expression hallmarks of cellular ageing. Biogerontology 19, 547–566 (2018).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Bahar, R. et al. Increased cell-to-cell variation in gene expression in ageing mouse heart. Nature 441, 1011–1014 (2006).

    Article  ADS  CAS  PubMed  Google Scholar 

  18. Wiley, C. D. et al. Analysis of individual cells identifies cell-to-cell variability following induction of cellular senescence. Aging Cell 16, 1043–1050 (2017).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Salzer, M. C. et al. Identity noise and adipogenic traits characterize dermal fibroblast aging. Cell 175, 1575–1590 (2018).

    Article  CAS  PubMed  Google Scholar 

  20. Croft, A. P. et al. Distinct fibroblast subsets drive inflammation and damage in arthritis. Nature 570, 246–251 (2019).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  21. Brady, J. J. et al. Early role for IL-6 signalling during generation of induced pluripotent stem cells revealed by heterokaryon RNA-seq. Nat. Cell Biol. 15, 1244–1252 (2013).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Keyes, B. E. et al. Impaired epidermal to dendritic T cell signaling slows wound repair in aged skin. Cell 167, 1323–1338 (2016).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Nishiguchi, M. A., Spencer, C. A., Leung, D. H. & Leung, T. H. Aging suppresses skin-derived circulating SDF1 to promote full-thickness tissue regeneration. Cell Rep. 24, 3383–3392 (2018).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Coppé, J. P. et al. Senescence-associated secretory phenotypes reveal cell-nonautonomous functions of oncogenic RAS and the p53 tumor suppressor. PLoS Biol. 6, e301 (2008).

    Article  CAS  PubMed Central  Google Scholar 

  25. Ramírez-Solis, R. et al. Genomic DNA microextraction: a method to screen numerous samples. Anal. Biochem. 201, 331–335 (1992).

    Article  PubMed  Google Scholar 

  26. Lo Sardo, V. et al. Influence of donor age on induced pluripotent stem cells. Nat. Biotechnol. 35, 69–74 (2017).

    Article  CAS  PubMed  Google Scholar 

  27. Kareta, M. S. et al. Inhibition of pluripotency networks by the Rb tumor suppressor restricts reprogramming and tumorigenesis. Cell Stem Cell 16, 39–50 (2015).

    Article  CAS  PubMed  Google Scholar 

  28. Vierbuchen, T. et al. Direct conversion of fibroblasts to functional neurons by defined factors. Nature 463, 1035–1041 (2010).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  29. Ackermann, M. & Strimmer, K. A general modular framework for gene set enrichment analysis. BMC Bioinformatics 10, 47 (2009).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Heinz, S. et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol. Cell 38, 576–589 (2010).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Langmead, B., Trapnell, C., Pop, M. & Salzberg, S. L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10, R25 (2009).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Hashimoto, T. B., Edwards, M. D. & Gifford, D. K. Universal count correction for high-throughput sequencing. PLOS Comput. Biol. 10, e1003494 (2014).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  33. Feng, J., Liu, T., Qin, B., Zhang, Y. & Liu, X. S. Identifying ChIP–seq enrichment using MACS. Nat. Protocols 7, 1728–1740 (2012).

    Article  CAS  PubMed  Google Scholar 

  34. Zhang, Y. et al. Model-based analysis of ChIP–seq (MACS). Genome Biol. 9, R137 (2008).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Ross-Innes, C. S. et al. Differential oestrogen receptor binding is associated with clinical outcome in breast cancer. Nature 481, 389–393 (2012).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  36. Li, Q., Brown, J. B., Huang, H. & Bickel, P. J. Measuring reproducibility of high-throughput experiments. Ann. Appl. Stat. 5, 1752–1779 (2011).

    Article  MathSciNet  MATH  Google Scholar 

  37. 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).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Benayoun, B. A. et al. H3K4me3 breadth is linked to cell identity and transcriptional consistency. Cell 158, 673–688 (2014).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Azuara, V. et al. Chromatin signatures of pluripotent cell lines. Nat. Cell Biol. 8, 532–538 (2006).

    Article  CAS  PubMed  Google Scholar 

  40. Bernstein, B. E. et al. A bivalent chromatin structure marks key developmental genes in embryonic stem cells. Cell 125, 315–326 (2006).

    Article  CAS  PubMed  Google Scholar 

  41. Guan, S., Price, J. C., Prusiner, S. B., Ghaemmaghami, S. & Burlingame, A. L. A data processing pipeline for mammalian proteome dynamics studies using stable isotope metabolic labeling. Mol. Cell. Proteomics 10, M111.010728 (2011).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Smith, C. A., Want, E. J., O’Maille, G., Abagyan, R. & Siuzdak, G. XCMS: processing mass spectrometry data for metabolite profiling using nonlinear peak alignment, matching, and identification. Anal. Chem. 78, 779–787 (2006).

    Article  CAS  PubMed  Google Scholar 

  43. van den Berg, R. A., Hoefsloot, H. C., Westerhuis, J. A., Smilde, A. K. & van der Werf, M. J. Centering, scaling, and transformations: improving the biological information content of metabolomics data. BMC Genomics 7, 142 (2006).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Xia, J., Sinelnikov, I. V., Han, B. & Wishart, D. S. MetaboAnalyst 3.0—making metabolomics more meaningful. Nucleic Acids Res. 43, W251–W257 (2015).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Bronte, V. & Zanovello, P. Regulation of immune responses by l-arginine metabolism. Nat. Rev. Immunol. 5, 641–654 (2005).

    Article  CAS  PubMed  Google Scholar 

  46. Phang, J. M., Liu, W., Hancock, C. & Christian, K. J. The proline regulatory axis and cancer. Front. Oncol. 2, 60 (2012).

    Article  PubMed  PubMed Central  Google Scholar 

  47. Rath, M., Müller, I., Kropf, P., Closs, E. I. & Munder, M. Metabolism via arginase or nitric oxide synthase: two competing arginine pathways in macrophages. Front. Immunol. 5, 532 (2014).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Fan, J. et al. Characterizing transcriptional heterogeneity through pathway and gene set overdispersion analysis. Nat. Methods 13, 241–244 (2016).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. O’Malley, J. et al. High-resolution analysis with novel cell-surface markers identifies routes to iPS cells. Nature 499, 88–91 (2013).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  50. Schaffer, B. E. et al. Identification of AMPK phosphorylation sites reveals a network of proteins involved in cell invasion and facilitates large-scale substrate prediction. Cell Metab. 22, 907–921 (2015).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Lam, M. T., Nauta, A., Meyer, N. P., Wu, J. C. & Longaker, M. T. Effective delivery of stem cells using an extracellular matrix patch results in increased cell survival and proliferation and reduced scarring in skin wound healing. Tissue Eng. Part A 19, 738–747 (2013).

    Article  CAS  PubMed  Google Scholar 

  52. Barret, J. P. et al. Accelerated re-epithelialization of partial-thickness skin wounds by a topical betulin gel: results of a randomized phase III clinical trials program. Burns 43, 1284–1294 (2017).

    Article  PubMed  Google Scholar 

  53. Butler, A., Hoffman, P., Smibert, P., Papalexi, E. & Satija, R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat. Biotechnol. 36, 411–420 (2018).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Stuart, T. et al. Comprehensive integration of single-cell data. Cell 177, 1888–1902 (2019).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Joost, S. et al. Single-cell transcriptomics of traced epidermal and hair follicle stem cells reveals rapid adaptations during wound healing. Cell Rep. 25, 585–597 (2018).

    Article  CAS  PubMed  Google Scholar 

  56. Cheng, J. B. et al. Transcriptional programming of normal and inflamed human epidermis at single-cell resolution. Cell Rep. 25, 871–883 (2018).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Der, E. et al. Tubular cell and keratinocyte single-cell transcriptomics applied to lupus nephritis reveal type I IFN and fibrosis relevant pathways. Nat. Immunol. 20, 915–927 (2019).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Sun, Z. et al. A Bayesian mixture model for clustering droplet-based single-cell transcriptomic data from population studies. Nat. Commun. 10, 1649 (2019).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  59. Mertens, J. et al. Directly reprogrammed human neurons retain aging-associated transcriptomic signatures and reveal age-related nucleocytoplasmic defects. Cell Stem Cell 17, 705–718 (2015).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Luminex. xPonent 3.1 Logistic Curve Fitting: Technical Notes. https://www.luminexcorp.com/zh/download/xponent-3-1-logistic-curve-fitting-white-paper/ (2014).

  61. Carter, A. C., Davis-Dusenbery, B. N., Koszka, K., Ichida, J. K. & Eggan, K. Nanog-independent reprogramming to iPSCs with canonical factors. Stem Cell Rep. 2, 119–126 (2014).

    Article  CAS  Google Scholar 

Download references

Acknowledgements

We thank D. Wagh from Stanford Functional Genomics Facility for help with 10x Genomics single-cell RNA-seq; C. Carswell-Crumpton and M. Weglarz from Stanford Shared FACS Facility for FACS support; L. Liu for help with shRNA experiments; V. Sebastiano and M. Kareta for help with iPS cell generation and quality assessment; V. Sebastiano, K. Andreasson, C. Weyand and T. Wyss-Coray for helpful discussions and input to the manuscript; L. Booth, T. Ruetz, P. P. Singh, X. Zhao, P. Navarro, B. Dulken and other members of the Brunet laboratory for helpful discussions and feedback on the manuscript; M. Buckley, B. Demmitt and P. P. Singh for helping with independent code checking; S. Chen, A. Freund, J. Goudeaux, M. Pesch, A. Roux and J. Reuter for feedback on the initial manuscript. The piggyback OSKM transposon constructs were a gift from K. Kaji. This work used the Genome Sequencing Service Center of the Stanford Center for Genomics and Personalized Medicine Sequencing Center, supported by the grant award NIH S10OD025212. Sorting was performed on instruments in the Stanford Shared FACS Facility obtained using the NIH grant award S10RR025518-01 and purchased by Parker Institute for Cancer Immunotherapy. Work was supported by NIH P01 GM099130 (to A.B., M.P.S. and J.W.), CIRM RB4-06087 (to A.B.), a generous philanthropic gift from M. and T. Barakett, the EMBO post-doctoral fellowship (to S.M.), the WennerGren post-doctoral fellowship (to S.M.), the Sweden-America Foundation post-doctoral fellowship (to S.M.), the Dean’s Fellowship (to E.M.), K99/R00 AG049934 (to B.A.B.), the Stanford Graduate Fellowship (to L.X.), the NSF Graduate Research Fellowship (to L.X.), R00 AG049934 and the Hanson-Thorell family fellowship (to B.A.B.), the Hagey Laboratory for Pediatric Regenerative Medicine, The Gunn/Olivier fund, the Johnson Longaker fund and R01 GM116892 (to M.T.L.).

Author information

Authors and Affiliations

Authors

Contributions

L.X. and A.M contributed equally to this study. S.M., E.M. and A.B. planned the study, with input from M.P.S. and J.W. S.M. and E.M. performed the reprogramming experiments, and analysed and interpreted data. S.M., E.M. and L.X. wrote the manuscript with the help of A.B. S.M. generated, processed and analysed bulk and single-cell RNA-seq datasets, analysed the metabolomics data, and performed most THY1-related and conditioned medium experiments. E.M. generated and propagated transgene-free iPS cell lines. All other studies were done by both E.M and S.M., unless otherwise noted. A.M. and S.M. performed wound healing experiments under the supervision of M.T.L. L.X. helped with reprogramming, FACS, and immunofluorescence experiments. F.J. generated the in vitro single-cell RNA-seq data under the supervision of M.P.S. R.S. generated the ChIP–seq libraries under the supervision of J.W. K.H. helped with statistics and PAGODA analysis. X.L. performed metabolomics experiments and helped with metabolomics data analysis and validation under the supervision of M.P.S. K.D. helped with reprogramming and western blotting experiments. L.P. helped with reprogramming and RT–qPCR experiments. C.E.A. and Y.S. performed the induced neuron reprogramming experiment under the supervision of M.W. B.A.B. helped with analysis of the epigenomic data. A.L.S.C. identified and collected the human samples. All authors discussed the results and commented on the manuscript.

Corresponding author

Correspondence to Anne Brunet.

Ethics declarations

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Peer review information Nature thanks Ellen Pure, Fuchou Tang and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

Extended data figures and tables

Extended Data Fig. 1 Primary old fibroblasts from mouse ear, mouse lungs and human skin secrete high levels of inflammatory cytokines, and the ability of individual cultures from ear fibroblasts to reprogram is stereotypical.

a, Cytokine profiling of plasma from young (3 months, n = 21) and old (29 months, n = 19) mice using Luminex multiplex cytokine assay (2 independent experiments). Box-and-whisker plot of log2-transformed fold change in MFI over median of young fibroblasts. Box plots depict the median and interquartile range, with whiskers indicating minimum and maximum values. *P < 0.05, **P < 0.01, ***P < 0.001, n.s., not significant; two-tailed Wilcoxon rank-sum test with Benjamini–Hochberg correction. Exact P values are in Supplementary Table 1a. b, Blood cell composition analysis of plasma from young (3 months, n = 9) and old (29 months, n = 9) mice using the Hemavet Multispecies Hematology Analyzer (2 independent experiments). Data are cell numbers per μl in whole-blood samples. Each dot represents cells from one mouse. Lines depict median. P values, two-tailed Wilcoxon rank-sum test with Benjamini–Hochberg correction. c, Percentage of fibroblasts (PDGFRα+) and immune cells in young (3 months, n = 8) and old (29 months, n = 8) fibroblast cultures at passage 3 (1 experiment), as determined by FACS using the indicated cell-type-specific surface markers. Data are mean ± s.e.m. Primary splenocytes from an 8-week-old mouse were used as a positive control (right). d, Comparison between cytokine profiles of plasma (red triangles) and conditioned medium from cultured ear fibroblasts at passage 3 (orange squares). Results are the mean log2-transformed concentrations (pg μl−1) of cytokines detected. For exact concentrations, see Supplementary Table 1a, b. n.d., not detected. e, Cytokine profiles of conditioned medium from primary cultures (passage 3) of lung fibroblasts from young (3 months, n = 8) and old (20–24 months, n = 9) mice (2 independent experiments). Box-and-whisker plot of log2-transformed fold change in MFI over median of young fibroblasts. *P < 0.05, **P < 0.01; two-tailed Wilcoxon rank-sum test with Benjamini–Hochberg correction. Box plots as in a. Exact P values are in Supplementary Table 1c. f, Cytokine profiles of conditioned medium from primary fibroblast cultures isolated from punch biopsy of pre-auricular skin of healthy human subjects of different ages. Results are shown as Spearman’s rank correlation coefficient (ρ) between donor age (years) and cytokine levels (MFI) in human fibroblast cultures (n = 8) at passage 3 (1 experiment). Each dot represents cells from one individual. P values, two-sided algorithm AS 89 in R. For multiple hypothesis testing, see Supplementary Table 1d. g, Cytokine profiles of conditioned medium collected from passage 23 cultures of iPS cell lines derived from young (3 months, n = 4) and old (29 months, n = 6) fibroblasts (1 experiment). Box-and-whisker plot of log2-transformed fold change in MFI over median of young iPS cells. Only cytokines that were detected at significantly different levels in young and old fibroblasts are shown (for a complete cytokine list and more details, see Supplementary Table 1e, f). Each dot represents an individual iPS cell line. P values, one-tailed Wilcoxon rank-sum test with Benjamini–Hochberg correction. Exact P values are in Supplementary Table 1e. Box plots as in a. h, Comparison of age-dependent changes in cytokine levels between plasma-incubated (described in a) and conditioned-medium-incubated mouse fibroblasts (described in Fig. 1b), their derived iPS cells (described in g) and from human fibroblasts (described in f) (Supplementary Table 1a, b, d, e) based on cytokines that are significantly different in conditioned medium from fibroblasts (Fig. 1b, bottom). Top (also presented in Fig. 1b), ranked fold change (old/young) in levels of the indicated cytokines in plasma, conditioned medium from mouse primary fibroblasts and iPS cells. Bottom, ranked Spearman ρ correlations for the indicated cytokines in conditioned medium from human primary fibroblasts (see f for individual ρ values). i, j, iPS cell lines derived from young and old mice show typical morphologies of mouse iPS cells, express similar levels of pluripotency markers and can give rise to cell types from all three germ layers upon embryoid body formation. i, Representative immunofluorescence images of iPS cell lines derived from young and old mice at passage 23, stained with the indicated antibodies (1 experiment). j, RT–qPCR on the indicated genes in embryoid bodies differentiated in vitro from iPS cell lines from young (n = 5) and old (n = 8) mice at passage 23 (1 experiment). Expression is presented as expression relative to the housekeeping gene Hprt1. Each bar represents one iPS cell line. k, Cytokine profiles of conditioned medium collected from cultures of young (n = 7) and old (n = 7) ear fibroblasts at passage 33 (1 experiment). Box-and-whisker plot of log2-transformed fold change in MFI over median of young fibroblasts. Only the cytokines that exhibited a significant difference in expression levels in conditioned medium from young and old fibroblasts at passage 3 are shown. P values, one-tailed Wilcoxon rank-sum test with Benjamini–Hochberg correction. For a complete cytokine list and exact P values, see Supplementary Table 2r. Box plots as in a. Note that the experiments in fibroblasts at passage 3 and 33 were conducted independently, and therefore statistical comparisons were restricted to within experiments. However, a direct comparison between the levels of secreted factors at passage 3 to 33 revealed that the levels of most cytokines decrease upon passaging. l, Reprogramming efficiency, assessed by SSEA1 staining, of young (n = 14) and old (n = 24) ear fibroblast cultures at passage 3 (3 independent experiments), as log2-transformed fold change over the median of young mice. Each dot represents a fibroblast culture from one mouse. P value, Fligner–Killeen test to assess differences in variance between age groups. m, Reprogramming efficiency assessed by AP staining of young (3 months, n = 7), middle-aged (12–13 months, n = 7) and old (28–30 months, n = 8) chest fibroblast cultures at passage 3 (2 independent experiments), as log2-transformed fold change over the median of young mice. Dots as in l. P values, Fligner-Killeen test to assess differences in variance between age groups with Benjamini–Hochberg correction. n, Reprogramming efficiency of fibroblast cultures are mainly stereotypical to fibroblast cultures from an individual mouse. Correlation plot depicting the reprogramming efficiency (assessed by AP staining) of fibroblast cultures reprogrammed in two experiments (separated by more than one month), with data from experiment 1 on the x axis and data from experiment 2 on the y axis. Data shown are from young (3 months, n = 14), middle-aged (12 months, n = 6) and old (29 months, n = 18) mice. Dots as in l. P values, two-sided algorithm AS 89 in R. There was a positive correlation (Spearman rank correlation, ρ = 0.63, P = 2.1 × 10−5) between reprogramming efficiencies of fibroblast cultures from the same mouse. o, Correlation plot depicting the ability of young (n = 4) and old (n = 5) fibroblast cultures at passage 3 to reprogram into neurons (iN, assessed by PSA-NCAM) or to iPS cells (iPSC, assessed by AP) (1 experiment). Dots as in l. P values, two-sided algorithm AS 89 in R. There was a significant positive correlation (Spearman rank correlation, ρ = 0.84, P = 0.003) between these two features. For individual experiments in a, b, e, ln, see Supplementary Table 7.

Extended Data Fig. 2 Old fibroblasts exhibit distinct transcriptomic, epigenomic and metabolomics profiles compared to young fibroblasts.

a, b, PCA of H3K4me3 (a) and H3K27me3 (b) peak intensities from young ear fibroblast cultures (3 months, n = 2), old cultures that reprogram well (good old, 29 months, n = 2) and old cultures that reprogram poorly (bad old, 29 months, n = 1) (1 experiment). Principal components (PCs) 1 and 2 are shown. c, PCA of metabolomes of young (n = 8), good old (n = 4) and bad old (n = 4) ear fibroblast cultures (1 experiment). Principal components 1 and 2 are shown. d, PCA of metabolomes of young (n = 8), good old (n = 4) and bad old (n = 4) ear fibroblast cultures (1 experiment). Principal components 1 and 3 are shown. e, Unsupervised hierarchical clustering of transcriptomes (RNA-seq) of young (n = 8), good old (n = 5) and bad old (n = 5) ear fibroblast cultures (3 independent experiments). Hierarchical clustering was performed using correlation-based dissimilarity (Pearson’s) as distance measure and average for linkage analysis. The y axis indicates the similarity between samples. be, Mice were the same ages as in a. fh, Unsupervised hierarchical clustering of H3K4me3 (f) and H3K27me3 (g) peaks and metabolomes (h) described in a, b and c, respectively. Hierarchical clustering was performed as in e. i, PCA of transcriptomes (RNA-seq) of good and bad old ear fibroblast cultures. Principal components 1 and 2 are shown. j, PCA of metabolomes of good old (n = 4) and bad old (n = 4) ear fibroblast cultures. Principal components 2 and 3 are shown. k, Selected GO terms enriched in the old transcriptomes (young, n = 8; old, n = 10), with corresponding FDR-adjusted P values (one-sided Fisher’s exact test with Benjamini–Hochberg correction). For a complete list, see Supplementary Table 2d. l, Heat map showing expression (VST-transformed read counts, scaled row-wise) of selected cytokine genes. The scale for expression fold changes is indicated on the left. m, Gene set enrichment analysis (GSEA) plot depicting the transcriptomes of young (n = 8) and old (n = 10) fibroblasts. Top, genes associated with fibroblast activation (see Supplementary Table 2f) are enriched in the old fibroblasts (P < 1.0 × 10−10, two-sided nominal P value). Bottom, heat map shows expression of fibroblast activation genes (VST-transformed read counts, scaled row-wise). The scale for expression fold changes is indicated on the left. nq, Analysis of the epigenomic data as described in a, b. In line with the transcriptomic data, age-dependent changes in epigenomic landscape also revealed enrichment of pathways involved in fibroblast activation, such as cytokines, extracellular matrix components and contractility-related features. n, Left, heat map shows the H3K4me3 peaks within promoter regions that exhibit a significant difference in intensity with age, assessed by Diffbind. Peak intensity is shown as VST-transformed read counts, scaled row-wise. The scale for peak intensity fold changes is indicated on the left. Right, selected enriched KEGG pathways colour coded according to significance (one-sided Fisher’s exact test with Benjamini–Hochberg correction; black, FDR-adjusted P < 0.05; grey, FDR-adjusted P < 0.15). For a complete list of significant KEGG terms, see Supplementary Table 2h. o, Top Venn diagram depicts the overlap of bivalent domains within promoters regions of young and old fibroblasts. For details, see ‘Chromatin immunoprecipitation followed by sequencing and analysis of the epigenomic landscape’. Middle, pie charts show how the unique bivalent domains in young fibroblasts change in old fibroblasts (left pie chart) and vice versa (right pie chart). Bottom, selected enriched KEGG pathways colour coded according to significance. For a complete list of KEGG terms, see Supplementary Table 2l. p, Top, Venn diagram depicting the overlap of broad H3K4me3 domains within promoters regions of young and old fibroblasts. For details see ‘Chromatin immunoprecipitation followed by sequencing and analysis of the epigenomic landscape’). Bottom, selected enriched KEGG pathways colour coded according to significance. For a complete list of significant KEGG terms, see Supplementary Table 2j. q, The relationship between changes in H3K4me3 peak intensity (as described in a) and gene expression for H3K4me3 peaks that are significantly different between the age groups (as described in Fig. 2b). The analysis was restricted to H3K4me3 peaks that are within promoter regions, defined as the transcriptional start site ±2 kb. The y axis denotes the log2-transformed fold change in gene expression between young and old fibroblasts for the gene, and x axis denotes the log2-transformed fold change of H3K4me3 peak intensity assigned to the gene. Genes of interest are labelled. r–t, Pathway analysis of all putatively identified metabolites that were significantly different between young and old fibroblasts (described in c), as well as the differentially expressed genes (described in Fig. 2b; FDR-adjusted P < 0.05, absolute fold change >1.5), using the MetaboAnalyst tool44. Note that the MetaboAnalyst tool does not provide multiple-hypothesis-corrected P values (one-sided Fisher’s exact test). s, Box plots showing the log2-transformed signal intensities of selected metabolites in the arginine and proline pathway from the metabolic profiling of young and old fibroblasts cultures at passage 3 (as described in c), for which the identity was confirmed using commercially available standards. Box plots depict the median and interquartile range, with whiskers indicating minimum and maximum values. *P < 0.05, **P < 0.01; two-tailed Wilcoxon rank-sum test with q value correction; l-arginine, P = 0.069; l-ornithine, P = 0.690; l-glutamate, P = 0.055; creatine, P = 0.022; creatinine, P = 0.002; putrescine, P = 0.081; spermidine P = 0.094; spermine, P = 0.016. t, Schematic representation of the biological functions of key metabolites and genes in the arginine and proline metabolic pathway, and how they relate to regulation of inflammatory cytokines and extracellular matrix synthesis45,46,47. Abundance of putative metabolites (oval) and gene transcripts (squares) in old fibroblasts is colour coded (red, higher in old; blue, lower in old; grey, not significantly different or not detected). Epigenomic changes are indicated with black asterisks. u, Top, top 3 motifs found in promoters of differentially expressed genes between young and old fibroblasts described in Fig. 2b, using HOMER motif analysis. Bottom, top 10 putative upstream regulators identified by the IPA database that are differentially expressed between young and old fibroblasts (FDR-adjusted P < 0.05, absolute fold change >1.5). Heat map depicts log2-transformed fold change in expression (old/young) calculated using DESeq2. The transcription factor identified across both analyses (EBF2) is in red. *P < 0.05, **P < 0.01, ***P < 0.001; cumulative hypergeometric distribution (Homer motif analysis), activation z&#x2010;score in IPA (upstream regulator analysis). For a complete list of significant motifs and upstream regulators, and exact P values, see Supplementary Table 2n. v, Top, heat map of differentially expressed genes in a regression analysis from young to old healthy human primary fibroblasts59 (n = 13, FDR-adjusted P < 0.05). Expression is shown as VST-transformed read counts, scaled row-wise. The depicted KEGG pathways are FDR-adjusted P < 0.15 (one-sided Fisher’s exact test with Benjamini–Hochberg correction) (Supplementary Table 2o, p). Bottom, VST-transformed expression of EBF2 across human samples as a function of age (years). Each square represents transcripts from a patient. NB, newborn. w, Pathway enrichment analysis of KEGG pathways associated with enhanced (up) or reduced (down) reprogramming, comparing H3K4me3 peak intensities between the top (n = 1) and the bottom (n = 1) old reprogramming cultures in our datasets (see Supplementary Table 2a). All depicted KEGG pathways were significantly enriched (FDR-adjusted P < 0.05). *P < 0.05; two-sided nominal P value with Benjamini–Hochberg correction. For a complete list of KEGG terms and exact P values, see Supplementary Table 3f.

Extended Data Fig. 3 Reprogramming erases features of inflammageing and variability between mice.

a–l, To test whether transcriptomic and metabolomics features of inflammageing could be erased by reprogramming, iPS cell lines from young and old fibroblasts at passage 23 were profiled for their transcriptome (RNA-seq) and metabolome. a, All iPS cell lines (n = 11) cluster with previously established bona fide iPS cell lines (n = 3) and ES cells (n = 4)61. Unsupervised hierarchical clustering based on overall transcriptomes of the indicated cell types. The hierarchical clustering was performed using correlation-based dissimilarity (Pearson’s) as distance measure and average for linkage analysis. The y axis indicates the similarity between samples. b, PCA of whole transcriptomes from RNA-seq data of iPS cell lines derived from young (3 months, n = 5) and old (29 months, n = 6) mice at passage 23. Principal components (PCs) 1 and 2 are shown. c, Unsupervised hierarchical clustering of transcriptomes described in b. Hierarchical clustering was performed as in a. d, Strip plot illustrating the log2-transformed fold expression changes of all genes with age for fibroblasts (left; described in Fig. 2b) and iPS cells (right; described in b). Genes detected as significantly upregulated or downregulated (DESeq2, FDR-adjusted P < 0.05, absolute fold change >1.5) with age, are shown in blue and yellow, respectively (see Supplementary Table 2q). e, PCA of metabolomes of iPS cell lines derived from young (n = 5) and old (n = 8) mice at passage 23. Ages as in b. Untargeted metabolomics profiles were generated using ultra-high performance liquid chromatography–mass spectrometry. Principal components 1 and 2 are shown. f, Unsupervised hierarchical clustering of metabolomics profiles described in e. Hierarchical clustering was performed as in a. g, Strip plot illustrating the log2-transformed fold change in signal intensity of all metabolic features with age for fibroblasts (left; described in Extended Data Fig. 2c) and iPS cells (right; described in e). Metabolic features detected as significantly up or down (using a two-tailed Wilcoxon rank-sum test with q value correction) with age are shown in blue and yellow, respectively. h, i, PCA (h) and unsupervised clustering (i) of iPS cell lines derived from young (n = 5) and old (n = 6) mice at passage 23, based on solely the genes that were significantly differentially expressed between young and old at fibroblast level. Ages as in b. Principal components 1 and 2 are shown. Hierarchical clustering was performed as in a. j–l, RT–qPCR of the indicated genes in fibroblasts cultures at passage 3 (j) and 33 (k), and iPS cell cultures at passage 23 (l). The genes shown represent the three major groups of features associated with fibroblast activation that change with age in fibroblasts. Box-and-whisker plot of log2-transformed fold change in MFI over median of young fibroblasts. Box plots depict the median and interquartile range, with whiskers indicating minimum and maximum values. Data are from young (n = 6) and old (n = 6) fibroblast cultures at passage 3, young (n = 5) and old (n = 6) fibroblast cultures at passage 33, young (n = 6) and old (n = 7) iPS cell cultures at passage 23. Ages as in b. *P < 0.05, **P < 0.01; one-tailed Wilcoxon rank-sum test with Benjamini–Hochberg correction. For j: Ccl7 (also known as Mcp3), P = 0.004; Ccl2 (also known as Mcp1), P = 0.004; Acvr2a (which encodes ACVR2α), P = 0.006; Ccl11 (also known as Eotaxin), P = 0.004; Pak6, P = 0.004; Thsb2, P = 0.004; Actn3, P = 0.006; Col1a1 (which encodes COL1α1), P = 0.004; Acta2 (which encodes αSMA), P = 0.004; Lama2, P = 0.004; Dmd, P = 0.008; F2r, P = 0.008. For k: Ccl7, P = 0.027; Ccl2, P = 0.027; Acvr2a, P = 0.027; Ccl11, P = 0.049; Pak6, P = 0.027; Thsb2, P = 0.026; Actn3, P = 0.026; Col1a1, P = 0.035; Acta2, P = 0.027; Lama2, P = 0.027; Dmd, P = 0.027; F2r, P = 0.229. For l: Ccl7, P = 0.800; Ccl2, P = 0.800; Acvr2a, P = 0.800; Ccl11, P = 1.000; Pak6, P = 0.800; Thsb2, P = 1.000; Actn3, P = 0.800; Col1a1, P = 1.000; Acta2, P = 1.000; Lama2, P = 1.000; Dmd, P = 1.000; F2r, P = 1.000. Note that the experiments were conducted independently in fibroblasts at passage 3, 33 and iPS cells, and therefore the statistical comparisons indicated were restricted to each independent experiment. However, a comparison between the expression of secreted factors at passage 3 to 33 shows that expression of Ccl11, but not Ccl2 and Ccl7, significantly decreases upon passaging.

Extended Data Fig. 4 Correlation between the fibroblast activation signature and reprogramming efficiency in single-cell RNA-seq data.

a, Single-cell RNA-seq of young (n = 30 cells), good and bad old (n = 31 cells) fibroblast cultures. t-SNE was performed on the VST-transformed read counts of all detected genes (analysed using DESeq2). Each dot represents a single fibroblast transcriptome. b, PAGODA of single-cell RNA-seq data from young and old fibroblasts at passage 3 performed using raw expression counts and all KEGG pathways, the in vitro fibroblast ageing, the fibroblast activation and de novo gene sets. Hierarchical clustering is based on 97 significantly overdispersed gene sets and the 405 genes driving the significantly overdispersed gene sets. Top, heat map of single cells from young and old fibroblast cultures. Middle, heat map of the separation of cells based on their principal component scores for the significantly overdispersed gene sets. Top heat map, PAGODA clustering of cells. Maroon and blue colours indicate increased and decreased expression of the associated gene sets, respectively. c, PAGODA as described in Extended Data Fig. 4c. Top, heat map of single cells from young and old fibroblast cultures. Middle, heat map of separation of cells based on their principal component scores for the significantly overdispersed gene sets. Top heat map, PAGODA clustering of cells. Maroon and blue colours indicate increased and decreased expression of the associated gene sets, respectively. d, PAGODA as described in c. Middle, heat map of separation of cells based on their principal component scores for the in vitro fibroblast ageing signature. Bottom, heat map of the expression of the genes that are part of the in vitro fibroblast ageing signature, and decrease with age; expression is shown as VST-transformed read counts, scaled row-wise. The scale for expression fold changes is indicated on the right. The bottom heat map indicates the cells that originate from good and bad old cultures. e, PAGODA as described in c. Middle, heat map of separation of cells based on their principal component scores for the in vitro fibroblast ageing signature. Bottom, heat map of expression of the genes that are part of the in vitro fibroblast ageing signature, and increase with age; expression is shown as VST-transformed read counts, scaled row-wise. The scale for expression fold changes is indicated on the right. The bottom heat map indicates the cells that originate from good and bad old cultures. f, PAGODA as described in c. Middle, heat map of the separation of cells based on their principal component scores for the fibroblast activation signature. Bottom, heat map of the expression of the genes that are part of the fibroblast activation gene set; expression is shown as VST-transformed read counts, scaled row-wise. The scale for expression fold changes is indicated on the right. The bottom heat map indicates the cells that originate from good and bad old cultures. g, PAGODA as described in c. Middle, heat map of the separation of cells based on their principal component scores for the KEGG cytokine–cytokine receptor interaction gene set. The heat map shows expression of the top 30 overdispersed genes in the KEGG cytokine–cytokine receptor interaction pathway; expression is shown as VST-transformed read counts, scaled row-wise. The scale for expression fold changes is indicated on the right. The bottom heat map indicates the cells that originate from good and bad old cultures.

Extended Data Fig. 5 Old fibroblast cultures are enriched for activated (THY1+PDGFRα+) fibroblasts, which are intrinsically poor at reprogramming but facilitate reprogramming extrinsically via secretion of cytokines.

a, Left, representative immunofluorescence images of young and old fibroblasts at passage 3 stained for αSMA (which is encoded by the Acta2 gene) . Right, quantification of the percentage of αSMA+ cells in young (3 months, n = 8) and old (29 months, n = 8) fibroblasts at passage 3 (1 experiment). Based on reprogramming efficiency, old cultures are shown as good old, bad old or old. Each dot represents cells from one mouse. Lines depict median. P value, two-tailed Wilcoxon rank-sum test. b, Left, representative immunofluorescence images of young and old fibroblasts at passage 3 incubated with EdU for 4 h, then stained for αSMA (red), EdU (green) and DAPI (blue). White arrows indicate an EdU-positive activated and a non-activated cell. Right, FACS quantification of the percentage of EdU-positive THY1PDGFRα+ (THY1) and THY1+PDGFRα+ (THY1+) cells in young (3 months, n = 5; 3 independent experiments) and old (29 months, n = 5; 1 experiment) cultures at passage 3. Dots and lines as in a. P values, two-tailed Wilcoxon rank-sum test. c, PAGODA of single-cell RNA-seq from young and old fibroblasts described in Extended Data Fig. 4c. Top heat map, PAGODA clustering of cells. Maroon and blue colours indicate increased and decreased expression of the associated gene sets, respectively. Middle heat map, expression of genes in the GO cellular senescence gene set, for which expression is shown as VST-transformed read counts, scaled row-wise. Bottom heat map, cells that originate from good and bad old cultures. The scale for expression fold changes is indicated on the right. d, RT–qPCR of p16Ink4a expression in cultures of THY1PDGFRα+ (THY1) and THY1+PDGFRα+ (THY1+) young and old cells at passages 4–6. Results are shown as fold change in expression over THY1PDGFRα+ cells. Data are from young (3 months, n = 5) and old (29 months, n = 5) cultures (6 independent experiments). One young and three old cultures were used in 2–3 independent experiments. In this case, an average of the measurements was determined. Dots and lines as in a. P values, two-tailed Wilcoxon signed-rank test. e, Percentage of SA-β-galactosidase-positive cells in young (3 months, n = 11) and old (28–29 months, n = 22) fibroblast cultures at passage 3 (3 independent experiments). log2-transformed fold change in SA-β-galactosidase-positive cells over median of young fibroblasts. Line indicates median. P values, two-tailed Wilcoxon ranked sum test. f, g, RT–qPCR of old THY1PDGFRα+ (THY1) and THY1+PDGFRα+ (THY1+) cells at passage 4–6 (3 independent experiments) untreated (f) (n = 6, 3 independent experiments) or treated with the indicated shRNA constructs for 72h (g) (n = 5, 4 independent experiments). Box-and-whisker plot of fold change in expression over THY1PDGFRα+ populations originating from the same culture (f) or over shLuciferase (shLuc) treated cells (g). Box plots depict the median and interquartile range, with whiskers indicating minimum and maximum values. *P < 0.06, one-tailed Wilcoxon signed-rank test with Benjamini–Hochberg correction. h, RT–qPCR on young (3 months, n = 5) fibroblasts after overexpression of Ebf2 for 48 h (2 independent experiments). Box-and-whisker plot of fold change in expression over cells treated with empty vector. Box plots as in f. *P < 0.06, one-tailed Wilcoxon signed-rank test with Benjamini–Hochberg correction. i, Heat map of significantly differentially expressed genes (determined by DESeq2) between freshly FACS-sorted THY1PDGFRα+Lin (THY1) and THY1+PDGFRα+Lin (THY1+) cells described in Fig. 3b and enriched KEGG pathways. Expression is shown as VST-transformed read counts, scaled row-wise. The scale for expression fold changes is indicated on the left. All depicted KEGG pathways were significantly enriched (one-sided Fisher’s exact test with Benjamini–Hochberg correction, FDR-adjusted P < 0.05). For a complete list of KEGG terms, see Supplementary Table 4e. j, Pathway enrichment analysis of KEGG pathways associated with ageing in dataset described in Fig. 3b. For a complete list of KEGG terms, see Supplementary Table 4g. **P = 0.01, ***P = 0.001; two-sided nominal P value with Benjamini–Hochberg correction. k, Top, ranked fold change (old/young) in levels of the indicated cytokines in plasma (see Extended Data Fig. 1a). Bottom, ranked fold change (old/young) in expression for the indicated cytokines in freshly FACS-sorted  THY1PDGFRα+Lin (THY1) and THY1+PDGFRα+Lin (THY1+) cells from young and old ears. See ‘Cytokine profiling analysis on plasma and conditioned medium using Luminex multi-analyte’ for calculation of ranked fold changes. Gene expression related to wounded fibroblasts is from datasets described in Extended Data Fig. 7e–g. ln, Correlation between the proportion of THY1+PDGFRα+ (THY1+) fibroblasts in old cultures (29 months, n = 23) (l), young (3 months, n = 21) and bad old (29 months, n = 6) cultures (m), and young (3 months, n = 21) and good old (29 months, n = 6) cultures (n), and the reprogramming efficiency of the culture (3 independent experiments). Dots as in a. P values, two-sided algorithm AS 89 in R. The y axis denotes the fold change in the proportion of THY1+PDGFRα+ fibroblasts relative to the median of young mice, and x axis denotes the fold change in reprogramming efficiency of the culture relative to the median of young mice. o, p, Correlation between the proliferation rate (o) or the percentage of SA-β-galactosidase-positive cells (p) of a given fibroblast culture and reprogramming efficiency of the culture. Proliferation rate was determined by calculating the growth slope of young (3 months, n = 15), middle-aged (12 months, n = 10) and old (28–29 months, n = 27) ear fibroblast cultures at passage 3 (4 independent experiments). Senescence was assessed by SA-β-galactosidase staining of young (3 months, n = 11), middle-aged (12 months, n = 11) and old (28–29 months, n = 22) ear fibroblast cultures at passage 3 (3 independent experiments). Dots as in a. P values, two-sided algorithm AS 89 in R. The y axis denotes the fold change in the proliferation rate or percentage of SA-β-galactosidase-positive cells relative to the median of young mice, and x axis denotes the fold change in reprogramming efficiency of the culture relative to the median of young mice. q, Reprogramming efficiency of FACS-sorted young (3 months, n = 8) and old (29 months, n = 7) THY1PDGFRα+ (THY1) and THY1+PDGFRα+ (THY1+) fibroblasts at passages 4–6, assessed using AP staining (3 independent experiments). log2-transformed fold change in reprogramming efficiency of the cells relative to the median of young THY1PDGFRα+ fibroblasts. One old culture was used in two independent experiments. In this case, an average of the measurements was determined. Dots and lines as in a. P values, two-tailed Wilcoxon signed-rank test. r, Reprogramming efficiency of FACS-sorted young (3 months, n = 5) THY1PDGFRα+ (THY1) and THY1+PDGFRα+ (THY1+) fibroblasts at passages 4–6, assessed as in q (3 independent experiments). Reprogramming was induced using a non-lentiviral piggyBac transposon system. Results are shown as number of AP+ colonies. Dots and lines as in a. P values, one-tailed Wilcoxon rank-sum test. s, t, Cytokine profiles of conditioned medium collected from cultures of old (s) (29 months, n = 6, 3 independent experiments) and young (t) (3 months, n = 6, 2 independent experiments) THY1PDGFRα+ (THY1) and THY1+PDGFRα+ (THY1+) fibroblasts at passages 4–6 . Comparisons were made between THY1PDGFRα+ and THY1+PDGFRα+ from the same original culture. Based on cytokines that are significantly different in conditioned medium from fibroblasts (Fig. 1b). Box-and-whisker plot of log2-transformed fold change in mean fluorescence intensity (MFI) over THY1PDGFRα+ fibroblasts. Box plots as in f. *P < 0.05, one-tailed Wilcoxon rank-sum test with Benjamini–Hochberg correction. Exact P values are in Supplementary Table 4h. u, Reprogramming efficiency, assessed as in q, of FACS-sorted young (3 months, n = 6) THY1PDGFRα+ (THY1) and THY1+PDGFRα+ (THY1+) fibroblasts at passages 4–6 treated with fresh conditioned medium daily starting from day 1 after infection (3 independent experiments). Conditioned medium was collected daily from the THY1PDGFRα+ or THY1+PDGFRα+ fibroblasts from the same original culture. log2-transformed fold change in reprogramming efficiency relative to the reprogramming efficiency of THY1PDGFRα+ fibroblasts treated with conditioned medium from THY1PDGFRα+ fibroblasts. One young culture was used in two independent experiments. In this case, an average of the measurements was determined. Dots and lines as in a. P values, two-tailed Wilcoxon signed-rank test. For individual experiments in b, dh, l-u and exact P values, see Supplementary Table 7.

Extended Data Fig. 6 Old fibroblasts secrete cytokines, including IL-6 and TNF, that induce inflammatory signalling pathways and modulate reprogramming efficiency.

a, Reprogramming efficiency, assessed using AP staining, of FACS-sorted young (3 months, n = 8) fibroblast cultures at passage 3 that were treated with conditioned medium from day 1 after infection (4 independent experiments). Conditioned medium was collected from young, good old or bad old fibroblast cultures. Results are shown as fold change in reprogramming efficiency relative to young fibroblasts treated with young conditioned medium. Each dot represents cells from one mouse. Lines depict median. P values, two-tailed Wilcoxon signed-rank test with Benjamini–Hochberg correction. b, Reprogramming efficiency, assessed as in a, of bad old (left, n = 7) or good old (right, n = 6) fibroblast cultures at passage 3 treated with conditioned medium from day 1 after infection (5 independent experiments). Conditioned medium was collected from good or bad old fibroblast cultures. Results are shown as fold change in reprogramming efficiency over old fibroblasts treated with bad old conditioned medium. Dots and lines as in a. P values, two-tailed Wilcoxon signed-rank test. c, Different representation of the data from Fig. 3e; each diamond represents the fold difference in reprogramming efficiency between a unique pair of good and bad old cultures (n = 8 pairs of good and bad old cultures, 5 independent experiments). Dots and lines as in a. P values, two-tailed Wilcoxon signed-rank test with Benjamini–Hochberg correction. d, Western blot analysis of young fibroblasts at passage 3 treated with the indicated cytokines at the concentration of 10 ng ml−1 for 30 min. Representative of 3 independent experiments. e, Reprogramming efficiency, assessed as in a, of young fibroblast cultures (3 months, n = 10) at passage 3 (3 independent experiments). Cells were treated with the indicated cytokines from day 1 after infection at the concentration of 10 ng ml−1. Results are shown as log2-transformed fold change in reprogramming efficiency over untreated cells. Dots and lines as in a. P values, two-tailed Wilcoxon signed-rank test with Benjamini–Hochberg correction. f, Reprogramming efficiency, assessed by SSEA1 staining, of young fibroblast cultures (3 months, n = 4) at passage 3 treated with the indicated cytokines from day 1 after infection at the concentration of 10 ng ml−1 (2 independent experiments). Results are shown as log2-transformed fold change in reprogramming efficiency over untreated cells. Dots and lines as in a. P values, one-tailed Wilcoxon signed-rank test with Benjamini–Hochberg correction. g, Reprogramming efficiency, assessed as in a, of old fibroblast cultures (29 months, n = 7) at passage 3 treated with the indicated cytokines from day 1 after infection at the concentration of 10 ng ml−1 (3 independent experiments). Results are shown as log2-transformed fold change in the reprogramming efficiency over untreated cells. Note that 1 old culture was used in 2 independent experiments. In this case, an average of the resultant measurements was determined. Dots and lines as in a. P values, two-tailed Wilcoxon signed-rank test with Benjamini–Hochberg correction. h, Reprogramming efficiency, assessed as in f, of old fibroblast cultures (29 months, n = 3) at passage 3 treated with the indicated cytokines from day 1 after infection at the concentration of 10 ng ml−1 (2 independent experiments). Results are shown as log2-transformed fold change in the reprogramming efficiency over untreated cells. Dots and lines as in a. P values, one-tailed Wilcoxon signed-rank test with Benjamini–Hochberg correction. i, Western blot analysis using the indicated antibodies of young fibroblasts at passage 3 treated with the indicated cytokines (10 ng ml−1) and blocking antibodies (8 μg ml−1) for 30 min. Cytokines were pretreated with either IgG or their corresponding blocking antibodies for 1 h before treatment. Representative of 2 independent experiments. j, Western blot analysis of old fibroblasts at passage 3 treated with the indicated cytokines at a concentration of 10 ng ml−1 for 30 min. k, Reprogramming efficiency, assessed as in a, of young fibroblast cultures (3 months, n = 4) at passage 3 treated with the indicated conditions from day 1 after infection (2 independent experiments). Cytokines (10 ng ml−1) were pretreated with either IgG or their corresponding blocking antibody (8 μg ml−1) for 1 h before treatment. Results are shown as log2-transformed fold change in reprogramming efficiency over untreated cells. Dots and lines as in a. P values, one-tailed Wilcoxon signed-rank test with Benjamini–Hochberg correction. l, m, Reprogramming efficiency, assessed as in a, of young fibroblast cultures (3 months, n = 6) at passage 3 treated with the indicated conditions from day 1 after infection (3 independent experiments). Conditioned medium was pretreated for 1 h with the indicated blocking antibody before administration. Results are shown as log2-transformed fold change relative to conditioned medium treated with IgG. Dots and lines as in a. P values, one-tailed Wilcoxon signed-rank test with Benjamini–Hochberg correction. n, Different representation of the data from Fig. 3f; each diamond represents the fold difference in reprogramming efficiency between a unique pair of good and bad old cultures (n = 6 pairs of good and bad old cultures, 4 independent experiments). Line marks median. P values, two-tailed Wilcoxon signed-rank test with Benjamini–Hochberg correction. o, Heat map showing the Spearman rank correlation coefficients between the levels of individual cytokines (top row), the ratio of the levels of TNF to other cytokines (middle row), the ratio of the levels of IL-6 to other cytokines (bottom row), and reprogramming efficiency in young (3 months, n = 19) and old (29 months, n = 18) cells (2 independent experiments). *P < 0.05, **P < 0.01; two-sided algorithm AS 89 in R; TNF:IL-6, P = 0.040; IL-6:IFNγ, P = 0.010; IL-6:IL-18, P = 0.040; IL-6:TNF, P = 0.040; IL-6:IL-10, P = 0.040; IL-6:CSF2, P = 0.005. The remaining P values can be found in Supplementary Table 7. p, q, The levels of IL-6 (p) or TNF (q) are not correlated with reprogramming efficiency. The y axis denotes the fold change in the levels of the indicated cytokine relative to the median of young mice and the x axis denotes the fold change in reprogramming efficiency of the culture relative to the median of young mice. Data are from young (3 months, n = 19) and old (29 months, n = 18) mice (2 independent experiments). P values, two-sided algorithm AS 89 in R. For individual experiments in ac, eh, kq, see Supplementary Table 7.

Extended Data Fig. 7 Ageing is associated with an increased variability in wound healing between old mice, and old fibroblasts in wounds are distinct from primary fibroblasts derived from healthy ear skin.

a, Example images of ear wounds of young mice, fast-healing old mice (fast old) and slow-healing old mice (slow old) at the indicated time points (2 independent experiments). Ink circles depict initial size of wounds. b, Ear wound healing curve from young (3–4 months, n = 26) and old (24–26 months, n = 28) mice (2 independent experiments). Full thickness wounds were induced on the dorsal side of both ears (see ‘Wounding and wound healing experiments’ for details) and the size of the wounds was assessed by imaging ear wounds every second day for 20 days. For each mouse, the average of both ear wounds was calculated. Graph depicts the average percentage of wound area remaining at the indicated time points. Data are mean ± s.e.m. c, Ear wound healing curves of the five fastest and the five slowest healing young and old mice. Graph depicts the average average of wound area remaining at the indicated time points. Data are mean ± s.e.m. d, FACS analysis as described in Fig. 3b to assess the percentage of THY1+PDGFRα+Lin (THY1+) cells in ears of young and old mice during basal conditions and at 7 days after induction of wounds. Results are shown as a percentage of THY1+PDGFRα+Lin cells over PDGFRα+Lin cells. Data shown are from young basal (3–4 months, n = 9 replicates, each with 2–3 mice), young wounded (3–4 months, n = 8 replicates, each with 2–3 mice), old basal (24–26 months, n = 10 replicates, each with 2–3 mice) and old wounded (24–26 months, n = 8 replicates, each with 2–3 mice) (3 independent experiments). Each dot represents a replicate with cells pooled from 2–3 mice. Line depicts median percentage. P values, two-tailed Wilcoxon rank-sum test. Note that the percentage of THY1+PDGFRα+Lin in young and old basal conditions is also presented in Fig. 3b. e, Pathway enrichment analysis based on population RNA-seq of young wounded (3–4 months, n = 6 replicates, each with 2–3 mice) and old wounded (24–26 months, n = 6 replicates, each with 2–3 mice) THY1PDGFRα+Lin and THY1+PDGFRα+Lin cells in vivo (1 experiment). The graph shows a subset of KEGG pathways that were found to be significantly enriched (FDR-adjusted P < 0.05). For a complete list of differentially expressed genes and pathways, with corresponding specific P values, see Supplementary Table 5c, d. **P < 0.01, ***P < 0.001; two-sided nominal P value with Benjamini–Hochberg correction. f, Comparison between the transcriptomic changes that occur in fibroblasts with age in vitro (as described in Fig. 2b) and in vivo (as described in Fig. 3b), as well as changes that occur upon wounding in young and old ears (as described in d). The heat map depicts the enrichment of the KEGG pathways that are present in at least two of the conditions described. For the complete list of differentially expressed genes and significant KEGG terms with specific P values, see Supplementary Tables 2b, c, 4f, g, 5a, b. The scale for enrichment is indicated on the left. g, Heat map of expression of a subset of cytokine genes from population RNA-seq of fibroblasts from young and old ears during basal and wounded conditions. Expression is shown as VST-transformed read counts, scaled row-wise. The scale for expression fold changes is indicated on the left. Basal and wound signatures refer to the average expression of the genes that are significantly downregulated or upregulated with wounding, respectively, in this dataset.

Extended Data Fig. 8 Single-cell RNA-seq analysis of fibroblasts in wounds from young and old mice and single-cell RNA-seq analysis of entire wounds from old slow- and fast-healing mice at day 7.

a, Quality control for 10x Genomics single-cell RNA-seq data of freshly isolated PDGFRα+Lin (CD45CD31EpCAMTER119TIE2) cells from wounds of young mice (3–4 months, cells pooled from n = 10 mice) or old mice (24–26 months, cells pooled from n = 10 mice), 7 days after induction of wounds. Number of unique genes (left), percentage of mitochondrial genes (middle) and number of unique molecular identifiers (UMIs, right) for each cell are shown, separated by age group. Each dot represents a single cell. b, Seurat analysis of all live high-quality PDGFRα+Lin cells described in a (3,036 cells in total) identified two main clusters of cells. Heat map depicts the expression of the top 10 marker genes for each significant cell cluster identified by Seurat, which are defined as the genes that are most specific to each population. The cell subpopulation identity assigned to each cluster is indicated below each column. c, PAGODA of the single-cell RNA-seq dataset described in a. PAGODA was performed using all KEGG pathways, and the in vitro fibroblast ageing and the fibroblast activation signatures (see Supplementary Table 2b, f). Top, heat map of single cells from wounds from young and old mice and cell clusters identified by Seurat and PAGODA analyses. Bottom, heat map of the separation of cells based on their principal component scores for the significantly overdispersed gene sets. Top heat map, PAGODA clustering of cells. Maroon and blue colours indicate increased and decreased expression of the associated gene sets, respectively. Bottom, log2-transformed fold change in the subpopulations between young and old wounds at day 7. d, PAGODA as described in c. Middle, heat map of the separation of cells based on their principal component scores for the fibroblast activation signature. Bottom, heat map of the expression of the genes that are part of the fibroblast activation signature (see Supplementary Table 2f); expression is shown as log-transformed and normalized gene expression values as calculated by Seurat and scaled row-wise. The scale for expression fold changes is indicated on the right. e, PAGODA as described in c. Middle, heat map of the separation of cells based on their principal component scores for the KEGG cytokine–cytokine receptor interaction gene set. Bottom, heat map of the expression of the genes that are part of the KEGG cytokine–cytokine receptor interaction; expression is shown as log-transformed and normalized gene expression values as calculated by Seurat and scaled row-wise. The scale for expression fold changes is indicated on the right. f, PAGODA as described in c. Middle, heat map of the separation of cells based on their principal component scores for the KEGG TNF signalling pathway gene set. Bottom, heat map of the expression of the genes that are part of the KEGG TNF signalling pathway; expression is shown as log-transformed and normalized gene expression values as calculated by Seurat and scaled row-wise. The scale for expression fold changes is indicated on the right. g, Representative images of the ears of the two slow-healing and two fast-healing old mice used for single-cell RNA-seq at day 7 after wounding (1 experiment). h, Ear wound healing curves of the ears of the two slow-healing and two fast-healing old mice used for single-cell RNA-seq. The percentage of the wound area that was not healed at the day 6 after induction of the wounds is indicated in parentheses. i, Quality control for 10x Genomics single-cell RNA-seq data of freshly isolated live cells from the ear wounds of slow-healing (n = 2) and fast-healing old mice (n = 2), 7 days after induction of wounds. Number of genes (left), percentage of mitochondrial genes (middle) and number of unique molecular identifier (right) for each cell are shown, separated by mouse. Each dot represents a single cell. j, Seurat analysis of all live high-quality cells described in i (10,797 cells in total) identified seven main clusters of cells. Heat map depicts the expression of the top 10 marker genes for each significant cell cluster identified by Seurat, which are defined as the genes that are most specific to each population. The cell subpopulation identity assigned to each cluster is indicated below each column.

Extended Data Fig. 9 Seurat and PAGODA single-cell RNA-seq analyses of fibroblasts identify distinct fibroblast subpopulations associated with fast- or slow-healing trajectories.

ag, Analysis of cells identified as fibroblasts from the single-cell RNA-seq dataset described in Fig. 4c. a, t-SNE clustering of cells identified as fibroblasts (2,678 cells in total) coloured by significant clusters identified using a k-nearest neighbour (KNN) graph-based algorithm as implemented by Seurat, or by mouse. b, log2-transformed fold change in the number of cells in each of the three subpopulations identified by Seurat between fast-healing old wounds and slow-healing old wounds. c, Seurat analysis of fibroblasts (2,678 cells in total) identified three main clusters. Heat map depicts the expression of the top 10 marker genes for each significant cell cluster identified by Seurat, which are defined as the genes that are most specific to each population. The identity of each cell subpopulation assigned to each cluster is indicated below each column. d, PAGODA of fibroblasts. PAGODA was performed using raw expression counts and all KEGG pathways, and the in vitro fibroblast ageing and the fibroblast activation signatures (see Supplementary Table 2b, f). Top, heat map of single cells from wounds of old mice and cell clusters identified by Seurat and PAGODA analyses. Bottom, heat map of the separation of cells based on their principal component scores for the significantly overdispersed gene sets. Top heat map, PAGODA clustering of cells. Maroon and blue colours indicate increased and decreased expression of the associated gene sets, respectively. e, PAGODA as described in d. Bottom, heat map of the expression of the genes that are part of the fibroblast activation signature (see Supplementary Table 2f); expression is shown as log-transformed and normalized gene expression values as calculated by Seurat and scaled row-wise. The scale for expression fold changes is indicated on the right. f, Expression of the genes that are part of the KEGG cytokine–cytokine receptor interaction gene set as in e. g, Expression of the genes that are part of the KEGG TNF signalling pathway as in e. h–l, Analysis of the combined single-cell RNA-seq datasets described in Fig. 4b, c. h, Seurat analysis of combined datasets clusters fibroblasts from both datasets together. t-SNE clustering of all live, high-quality cells from both datasets (13,833 cells in total) coloured by significant clusters identified using a KNN graph-based algorithm as implemented by Seurat, or by mouse. i, t-SNE clustering of combined fibroblasts from the datasets described in Fig. 4b (PDGFRα+Lin) and Fig. 4c. Combined fibroblasts (5,716 cells in total) are coloured by significant clusters identified using a KNN graph-based algorithm as implemented by Seurat, or by mouse. j, Seurat analysis of combined fibroblasts (5,716 cells in total) identified three main subpopulations. Heat map depicts the expression of the top 10 marker genes for each significant subpopulation identified by Seurat, which are defined as the genes that are most specific to each population. The cell subpopulation identity assigned to each cluster is indicated below each column. k, PAGODA of combined fibroblasts. PAGODA was performed using Seurat normalized counts and all KEGG pathways, the in vitro fibroblast ageing, the fibroblast activation signatures (see Supplementary Table 2b, f). Top, heat map of single fibroblasts from wounds of young and old mice or wounds from old fast- or slow-healing mice, and cell clusters identified by Seurat and PAGODA analyses. Bottom, heat map of separation of cells based on their principal component scores for a subset of the top significantly overdispersed gene sets. Top heat map, PAGODA clustering of cells. Maroon and blue colours indicate increased and decreased expression of the associated gene sets, respectively. Note that fibroblast subpopulation B did not contain cells from old/young in the combined analysis. This is probably owing to the fact that this subpopulation of fibroblast has some markers of the haematopoietic lineage, and is probably depleted in the PFGDRα+Lin FACS-sorting scheme used to isolate fibroblasts from the wounds of young and old mice. l, log2-transformed fold change in each of the three combined fibroblast subpopulations identified by PAGODA between wounds of old fast- and slow-healing mice, or between wounds from young and old mice, at day 7.

Extended Data Fig. 10 FACS schematic for in vitro and in vivo fibroblast analysis and sorting and full western blot membranes.

a, FACS schematic for analysis and sorting of THY1PDGFRα+ and THY1+PDGFRα+ cells in young and old cultures at passage 3. Gates shown on each plot are indicated above the plot. Marker and fluorophore are shown on each axis. FMO, fluorescence minus one. b, FACS schematic for analysis and sorting of live THY1PDGFRα+Lin and THY1+PDGFRα+Lin cells from young and old fresh tissues (ears), used for population RNA-seq and single-cell RNA-seq analyses. Gates shown on each plot are indicated above the plot. Marker and fluorophore are shown on each axis. c, Full western blot membranes from Extended Data Fig. 6d, i, j. Boxes indicate the cropped area.

Supplementary information

Reporting Summary

Supplementary Table 1

Cytokine profiles of young and old plasma, and of conditioned media from young and old fibroblasts at passage 3.

Supplementary Table 2

Analysis of transcriptomic, epigenomic, and metabolomic profiles of young and old fibroblasts at passage 3.

Supplementary Table 3

Differential gene expression and pathway enrichment analyses of good and bad reprogramming cultures.

Supplementary Table 4

Differential gene expression and pathway enrichment analyses of non-activated and activated fibroblasts in vitro and in vivo.

Supplementary Table 5

Differential gene expression and pathway enrichment analyses of in vivo fibroblasts during basal and wound condition.

Supplementary Table 6

List of primers used in this study, and information on RNA-seq data used in this study.

Supplementary Table 7

Experimental list.

Rights and permissions

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Mahmoudi, S., Mancini, E., Xu, L. et al. Heterogeneity in old fibroblasts is linked to variability in reprogramming and wound healing. Nature 574, 553–558 (2019). https://doi.org/10.1038/s41586-019-1658-5

Download citation

  • Received:

  • Accepted:

  • Published:

  • Issue Date:

  • DOI: https://doi.org/10.1038/s41586-019-1658-5

This article is cited by

Comments

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.

Search

Quick links

Nature Briefing

Sign up for the Nature Briefing newsletter — what matters in science, free to your inbox daily.

Get the most important science stories of the day, free in your inbox. Sign up for Nature Briefing