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:

Spatial multi-omic map of human myocardial infarction

Abstract

Myocardial infarction is a leading cause of death worldwide1. Although advances have been made in acute treatment, an incomplete understanding of remodelling processes has limited the effectiveness of therapies to reduce late-stage mortality2. Here we generate an integrative high-resolution map of human cardiac remodelling after myocardial infarction using single-cell gene expression, chromatin accessibility and spatial transcriptomic profiling of multiple physiological zones at distinct time points in myocardium from patients with myocardial infarction and controls. Multi-modal data integration enabled us to evaluate cardiac cell-type compositions at increased resolution, yielding insights into changes of the cardiac transcriptome and epigenome through the identification of distinct tissue structures of injury, repair and remodelling. We identified and validated disease-specific cardiac cell states of major cell types and analysed them in their spatial context, evaluating their dependency on other cell types. Our data elucidate the molecular principles of human myocardial tissue organization, recapitulating a gradual cardiomyocyte and myeloid continuum following ischaemic injury. In sum, our study provides an integrative molecular map of human myocardial infarction, represents an essential reference for the field and paves the way for advanced mechanistic and therapeutic studies of cardiac disease.

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: Spatial multi-omic profiling of human myocardial infarction.
Fig. 2: Characterization of tissue organization using spatial transcriptomics data.
Fig. 3: Characteristic tissue structures inferred from spatial transcriptomics data.
Fig. 4: Sub-clustering of cardiomyocytes.
Fig. 5: Sub-clustering of endothelial cells.
Fig. 6: Characterization of mesenchymal–myeloid interaction.

Similar content being viewed by others

Data availability

Processed snRNA-seq, snATAC-seq, and spatial transcriptomics data are available at cellxgene https://cellxgene.cziscience.com/collections/8191c283-0816-424b-9b61-c3e1d6258a77 and at the Zenodo data archive (https://zenodo.org/record/6578047). Raw data generated by CellRanger and SpaceRanger pipelines are available through the Human Cell Atlas Data Portal at https://data.humancellatlas.org/explore/projects/e9f36305-d857-44a3-93f0-df4e6007dc97 and at the Zenodo data archive (https://zenodo.org/record/6578553, https://zenodo.org/record/6578617 and https://zenodo.org/record/6580069). Source data are provided with this paper.

Code availability

All code used for analysis is available at https://github.com/saezlab/visium_heart and https://github.com/KramannLab/visium_heart.

References

  1. Wong, N. D. Epidemiological studies of CHD and the evolution of preventive cardiology. Nat. Rev. Cardiol. 11, 276–289 (2014).

    Article  PubMed  Google Scholar 

  2. Niccoli, G. et al. Optimized treatment of ST-elevation myocardial infarction. Circ. Res. 125, 245–258 (2019).

    Article  CAS  PubMed  Google Scholar 

  3. Prabhu Sumanth, D. & Frangogiannis Nikolaos, G. The biological basis for cardiac repair after myocardial infarction. Circ. Res. 119, 91–112 (2016).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Litviňuková, M. et al. Cells of the adult human heart. Nature 588, 466–472 (2020).

    Article  ADS  PubMed  PubMed Central  CAS  Google Scholar 

  5. Wang, L. et al. Single-cell reconstruction of the adult human heart during heart failure and recovery reveals the cellular landscape underlying cardiac function. Nat. Cell Biol. 22, 108–119 (2020).

    Article  PubMed  CAS  Google Scholar 

  6. Tucker, N. R. et al. Transcriptional and cellular diversity of the human heart. Circulation 142, 466–482 (2020).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Dodou, E., Verzi, M. P., Anderson, J. P., Xu, S.-M. & Black, B. L. Mef2c is a direct transcriptional target of ISL1 and GATA factors in the anterior heart field during mouse embryonic development. Development 131, 3931–3942 (2004).

    Article  CAS  PubMed  Google Scholar 

  8. Yan, C., Zhu, M., Staiger, J., Johnson, P. F. & Gao, H. C5a-regulated CCAAT/enhancer-binding proteins β and δ are essential in Fcγ receptor-mediated inflammatory cytokine and chemokine production in macrophages. J. Biol. Chem. 287, 3217–3230 (2012).

    Article  CAS  PubMed  Google Scholar 

  9. Kovary, K. & Bravo, R. The jun and fos protein families are both required for cell cycle progression in fibroblasts. Mol. Cell. Biol. 11, 4466–4472 (1991).

    CAS  PubMed  PubMed Central  Google Scholar 

  10. Li, S., Wang, D.-Z., Wang, Z., Richardson, J. A. & Olson, E. N. The serum response factor coactivator myocardin is required for vascular smooth muscle development. Proc. Natl Acad. Sci. USA 100, 9366–9370 (2003).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  11. Ge, Y. et al. Switching macrophage gene expression from inflammation-resolution to hemorrhage-resolution by redirection of activating transcription factor 1 (ATF1) binding by SMARCA4, BACH1 and histone H3K9 acetylation. Atherosclerosis 315, e2 (2020).

    Article  Google Scholar 

  12. Pirruccello, J. P. et al. Analysis of cardiac magnetic resonance imaging in 36,000 individuals yields genetic insights into dilated cardiomyopathy. Nat. Commun. 11, 2254 (2020).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  13. Ruetten, H., Dimmeler, S., Gehring, D., Ihling, C. & Zeiher, A. M. Concentric left ventricular remodeling in endothelial nitric oxide synthase knockout mice by chronic pressure overload. Cardiovasc. Res. 66, 444–453 (2005).

    Article  CAS  PubMed  Google Scholar 

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

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  15. Pakshir, P. et al. Dynamic fibroblast contractions attract remote macrophages in fibrillar collagen matrix. Nat. Commun. 10, 1850 (2019).

    Article  ADS  PubMed  PubMed Central  CAS  Google Scholar 

  16. Aoyagi, T. & Matsui, T. Phosphoinositide-3 kinase signaling in cardiac hypertrophy and heart failure. Curr. Pharm. Des. 17, 1818–1824 (2011).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Mak, T. W., Hauck, L., Grothe, D. & Billia, F. p53 regulates the cardiac transcriptome. Proc. Natl Acad. Sci. USA 114, 2331–2336 (2017).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Bergmann, M. W. WNT signaling in adult cardiac hypertrophy and remodeling: lessons learned from cardiac development. Circ. Res. 107, 1198–1208 (2010).

    Article  CAS  PubMed  Google Scholar 

  19. Ahern, B. M. et al. Myocardial-restricted ablation of the GTPase RAD results in a pro-adaptive heart response in mice. J. Biol. Chem. 294, 10913–10927 (2019).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Lotteau, S. et al. Acute genetic ablation of cardiac sodium/calcium exchange in adult mice: implications for cardiomyocyte calcium regulation, cardioprotection, and arrhythmia. J. Am. Heart Assoc. 10, e019273 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Fernandez-Caggiano, M. et al. Mitochondrial pyruvate carrier abundance mediates pathological cardiac hypertrophy. Nat. Metab. 2, 1223–1231 (2020).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Hama, N. et al. Rapid ventricular induction of brain natriuretic peptide gene expression in experimental acute myocardial infarction. Circulation 92, 1558–1564 (1995).

    Article  CAS  PubMed  Google Scholar 

  23. Waspe, L. E., Ordahl, C. P. & Simpson, P. C. The cardiac β-myosin heavy chain isogene is induced selectively in α1-adrenergic receptor-stimulated hypertrophy of cultured rat heart myocytes. J. Clin. Invest. 85, 1206–1214 (1990).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Beggah, A. T. et al. Reversible cardiac fibrosis and heart failure induced by conditional expression of an antisense mRNA of the mineralocorticoid receptor in cardiomyocytes. Proc. Natl Acad. Sci. USA 99, 7160–7165 (2002).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  25. Bakker, M. L. et al. T-box transcription factor TBX3 reprogrammes mature cardiac myocytes into pacemaker-like cells. Cardiovasc. Res. 94, 439–449 (2012).

    Article  CAS  PubMed  Google Scholar 

  26. Jiang, J. et al. Cardiac myosin binding protein C regulates postnatal myocyte cytokinesis. Proc. Natl Acad. Sci. USA 112, 9046–9051 (2015).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  27. Piroddi, N. et al. Myocardial overexpression of ANKRD1 causes sinus venosus defects and progressive diastolic dysfunction. Cardiovasc. Res. 116, 1458–1472 (2020).

    Article  CAS  PubMed  Google Scholar 

  28. Hill, C. et al. Inhibition of AP-1 signaling by JDP2 overexpression protects cardiomyocytes against hypertrophy and apoptosis induction. Cardiovasc. Res. 99, 121–128 (2013).

    Article  CAS  PubMed  Google Scholar 

  29. Kramann, R. et al. Perivascular Gli1+ progenitors are key contributors to injury-induced organ fibrosis. Cell Stem Cell 16, 51–66 (2015).

    Article  CAS  PubMed  Google Scholar 

  30. Hulsmans, M. et al. Macrophages facilitate electrical conduction in the heart. Cell 169, 510–522.e20 (2017).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Vieira, J. M. et al. The cardiac lymphatic system stimulates resolution of inflammation following myocardial infarction. J. Clin. Invest. 128, 3402–3412 (2018).

    Article  PubMed  PubMed Central  Google Scholar 

  32. Armulik, A., Abramsson, A. & Betsholtz, C. Endothelial/pericyte interactions. Circ. Res. 97, 512–523 (2005).

    Article  CAS  PubMed  Google Scholar 

  33. Deshpande, S. S., Angkeow, P., Huang, J., Ozaki, M. & Irani, K. Rac1 inhibits TNF‐α‐induced endothelial cell apoptosis: dual regulation by reactive oxygen species. FASEB J. 14, 1705–1714 (2000).

    Article  CAS  PubMed  Google Scholar 

  34. Kuppe, C. et al. Decoding myofibroblast origins in human kidney fibrosis. Nature 589, 281–286 (2021).

    Article  ADS  CAS  PubMed  Google Scholar 

  35. Li, Z. et al. Chromatin-accessibility estimation from single-cell ATAC-seq data with scOpen. Nat. Commun. 12, 6386 (2021).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  36. Bugg, D. et al. MBNL1 drives dynamic transitions between fibroblasts and myofibroblasts in cardiac wound healing. Cell Stem Cell 29, 419–433.e10 (2022).

    Article  CAS  PubMed  Google Scholar 

  37. Davis, J. et al. MBNL1-mediated regulation of differentiation RNAs promotes myofibroblast transformation and the fibrotic response. Nat. Commun. 6, 10084 (2015).

    Article  ADS  CAS  PubMed  Google Scholar 

  38. Kramann, R. et al. Pharmacological GLI2 inhibition prevents myofibroblast cell-cycle progression and reduces kidney fibrosis. J. Clin. Invest. 125, 2935–2951 (2015).

    Article  PubMed  PubMed Central  Google Scholar 

  39. Lawler, J. Thrombospondin-1 as an endogenous inhibitor of angiogenesis and tumor growth. J. Cell. Mol. Med. 6, 1–12 (2002).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Alexanian, M. et al. A transcriptional switch governs fibroblast activation in heart disease. Nature 595, 438–443 (2021).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  41. Forte, E. et al. Cross-priming dendritic cells exacerbate immunopathology after ischemic tissue damage in the heart. Circulation 143, 821–836 (2021).

    Article  CAS  PubMed  Google Scholar 

  42. Dick, S. A. et al. Three tissue resident macrophage subsets coexist across organs with conserved origins and life cycles. Sci. Immunol. 7, eabf7777 (2022).

    Article  CAS  PubMed  Google Scholar 

  43. MacDonald, L. et al. COVID-19 and RA share an SPP1 myeloid pathway that drives PD-L1+ neutrophils and CD14+ monocytes. JCI Insight 6, e147413 (2021).

    Article  PubMed Central  Google Scholar 

  44. Morse, C. et al. Proliferating SPP1/MERTK-expressing macrophages in idiopathic pulmonary fibrosis. Eur. Respir. J. 54, 1802441 (2019).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Bevan, L. et al. Specific macrophage populations promote both cardiac scar deposition and subsequent resolution in adult zebrafish. Cardiovasc. Res. 116, 1357–1371 (2020).

    Article  CAS  PubMed  Google Scholar 

  46. DeLeon-Pennell, K. Y. et al. CD36 is a matrix metalloproteinase-9 substrate that stimulates neutrophil apoptosis and removal during cardiac remodeling. Circ. Cardiovasc. Genet. 9, 14–25 (2016).

    Article  CAS  PubMed  Google Scholar 

  47. Ismahil, M. A. et al. Remodeling of the mononuclear phagocyte network underlies chronic inflammation and disease progression in heart failure: critical importance of the cardiosplenic axis. Circ. Res. 114, 266–282 (2014).

    Article  CAS  PubMed  Google Scholar 

  48. Shiraishi, M. et al. Alternatively activated macrophages determine repair of the infarcted adult murine heart. J. Clin. Invest. 126, 2151–2166 (2016).

    Article  PubMed  PubMed Central  Google Scholar 

  49. Bergmann, O. et al. Dynamics of cell generation and turnover in the human heart. Cell 161, 1566–1575 (2015).

    Article  CAS  PubMed  Google Scholar 

  50. Yekelchyk, M., Guenther, S., Preussner, J. & Braun, T. Mono- and multi-nucleated ventricular cardiomyocytes constitute a transcriptionally homogenous cell population. Basic Res. Cardiol. 114, 36 (2019).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  51. Habib, N. et al. Massively parallel single-nucleus RNA-seq with DroNc-seq. Nat. Methods 14, 955–958 (2017).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Koenig, A. L. et al. Single-cell transcriptomics reveals cell-type-specific diversification in human heart failure. Nat. Cardiovasc. Res. 1, 263–280 (2022).

  53. Tippani, M. et al. VistoSeg: processing utilities for high-resolution Visium/Visium-IF images for spatial transcriptomics data. bioRxiv https://www.biorxiv.org/content/10.1101/2021.08.04.452489v2 (2022).

  54. Curaj, A., Simsekyilmaz, S., Staudt, M. & Liehn, E. et al. Minimal invasive surgical procedure of inducing myocardial infarction in mice. J. Vis. Exp. 99, e52197 (2015).

    Google Scholar 

  55. Linkert, M. et al. Metadata matters: access to image data in the real world. J. Cell Biol. 189, 777–782 (2010).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Bahry, E. et al. RS-FISH: precise, interactive, fast, and scalable FISH spot detection. Preprint at bioRxiv https://doi.org/10.1101/2021.03.09.434205 (2021).

  57. Greenwald, N. F. et al. Whole-cell segmentation of tissue images with human-level performance using large-scale data annotation and deep learning. Nat. Biotechnol. 40, 555–565 (2021).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  58. Schneider, C. A., Rasband, W. S. & Eliceiri, K. W. NIH Image to ImageJ: 25 years of image analysis. Nat. Methods 9, 671–675 (2012).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Hao, Y. et al. Integrated analysis of multimodal single-cell data. Cell 184, 3573–3587.e29 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Germain, P.-L., Lun, A., Macnair, W. & Robinson, M. D. Doublet identification in single-cell sequencing data using scDblFinder. F1000Res. 10, 979 (2021).

    Article  PubMed  Google Scholar 

  61. O’Flanagan, C. H. et al. Dissociation of solid tumor tissues with cold active protease for single-cell RNA-seq minimizes conserved collagenase-associated stress responses. Genome Biol. 20, 210 (2019).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  62. Korsunsky, I. et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat. Methods 16, 1289–1296 (2019).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. McCarthy, D. J., Chen, Y. & Smyth, G. K. Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 40, 4288–4297 (2012).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Wolf, F. A., Angerer, P. & Theis, F. J. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 19, 15 (2018).

    Article  PubMed  PubMed Central  Google Scholar 

  65. Granja, J. M. et al. ArchR is a scalable software package for integrative single-cell chromatin accessibility analysis. Nat. Genet. 53, 403–411 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

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

    Article  CAS  PubMed  PubMed Central  Google Scholar 

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

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  68. Li, Z. et al. Identification of transcription factor binding sites using ATAC-seq. Genome Biol. 20, 45 (2019).

    Article  PubMed  PubMed Central  Google Scholar 

  69. Fornes, O. et al. JASPAR 2020: update of the open-access database of transcription factor binding profiles. Nucleic Acids Res. 48, D87–D92 (2020).

    CAS  PubMed  Google Scholar 

  70. Badia-i-Mompel, P. et al. decoupleR: ensemble of computational methods to infer biological activities from omics data. Bioinformatics Adv. 2, vbac016 (2022).

    Article  Google Scholar 

  71. Chang, C. C. et al. Second-generation PLINK: rising to the challenge of larger and richer datasets. Gigascience 4, 7 (2015).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  72. Ulirsch, J. C. et al. Interrogation of human hematopoiesis at single-cell and single-variant resolution. Nat. Genet. 51, 683–693 (2019).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Wolock, S. L., Lopez, R. & Klein, A. M. Scrublet: computational identification of cell doublets in single-cell transcriptomic data. Cell Syst. 8, 281–291.e9 (2019).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Hansen, B. B. & Klopfer, S. O. Optimal full matching and related designs via network flows. J. Comput. Graph. Stat. 15, 609–627 (2006).

    Article  MathSciNet  Google Scholar 

  75. Haghverdi, L., Büttner, M., Wolf, F. A., Buettner, F. & Theis, F. J. Diffusion pseudotime robustly reconstructs lineage branching. Nat. Methods 13, 845–848 (2016).

    Article  CAS  PubMed  Google Scholar 

  76. Schep, A. N., Wu, B., Buenrostro, J. D. & Greenleaf, W. J. chromVAR: inferring transcription-factor-associated accessibility from single-cell epigenomic data. Nat. Methods 14, 975–978 (2017).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Freeman, L. C. Centrality in social networks conceptual clarification. Soc. Networks 1, 215–239 (1978).

    Article  Google Scholar 

  78. Page L. et al. The PageRank Citation Ranking: Bringing Order to the Web. (Stanford IfoLab, 1999).

    Google Scholar 

  79. Hafemeister, C. & Satija, R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biol. 20, 296 (2019).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  80. Kleshchevnikov, V. et al. Cell2location maps fine-grained cell types in spatial transcriptomics. Nat. Biotechnol. 40, 661–671 (2022).

    Article  CAS  PubMed  Google Scholar 

  81. Schubert, M. et al. Perturbation-response genes reveal signaling footprints in cancer gene expression. Nat. Commun. 9, 20 (2018).

    Article  ADS  PubMed  PubMed Central  CAS  Google Scholar 

  82. Holland, C. H., Szalai, B. & Saez-Rodriguez, J. Transfer of regulatory knowledge from human to mouse for functional genomics analysis. Biochim. Biophys. Acta 1863, 194431 (2020).

    Article  CAS  Google Scholar 

  83. Zhu, J., Sun, S. & Zhou, X. SPARK-X: non-parametric modeling enables scalable and robust detection of spatial expression patterns for large spatial transcriptomic studies. Genome Biol. 22, 184 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Liberzon, A. et al. Molecular signatures database (MSigDB) 3.0. Bioinformatics 27, 1739–1740 (2011).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  85. Gillespie, M. et al. The reactome pathway knowledgebase 2022. Nucleic Acids Res. 50, D687–D692 (2021).

    Article  PubMed Central  CAS  Google Scholar 

  86. Grant, A. O. Cardiac ion channels. Circ. Arrhythm. Electrophysiol. 2, 185–194 (2009).

    Article  PubMed  Google Scholar 

  87. Tanevski, J., Flores, R. O. R., Gabor, A., Schapiro, D. & Saez-Rodriguez, J. Explainable multiview framework for dissecting spatial relationships from highly multiplexed data. Genome Biol. 23, 97 (2022).

  88. Pawlowsky-Glahn V. & Buccianti A. Compositional Data Analysis: Theory and Applications. (John Wiley & Sons, 2011).

    Book  MATH  Google Scholar 

  89. Lun, A. T. L., McCarthy, D. J. & Marioni, J. C. A step-by-step workflow for low-level analysis of single-cell RNA-seq data with Bioconductor. F1000Res. 5, 2122 (2016).

    PubMed  PubMed Central  Google Scholar 

  90. Dimitrov, D. et al. Comparison of methods and resources for cell-cell communication inference from single-cell RNA-Seq data. Nat. Commun. 13, 3224 (2022).

  91. Efremova, M., Vento-Tormo, M., Teichmann, S. A. & Vento-Tormo, R. CellPhoneDB: inferring cell–cell communication from combined expression of multi-subunit ligand-receptor complexes. Nat. Protoc. 15, 1484–1506 (2020).

    Article  CAS  PubMed  Google Scholar 

  92. Türei, D. et al. Integrated intra- and intercellular signaling knowledge for multicellular omics analysis. Mol. Syst. Biol. 17, e9923 (2021).

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  93. Nagai, J. S., Leimkühler, N. B., Schaub, M. T., Schneider, R. K. & Costa, I. G. CrossTalkeR: analysis and visualization of ligand–receptor networks. Bioinformatics 37, 4263–4265 (2021).

    Article  CAS  PubMed  Google Scholar 

  94. Alquicira-Hernandez, J. & Powell, J. E. Nebulosa recovers single cell gene expression signals by kernel density estimation. Bioinformatics 37, 2485–2487 (2021).

    Article  CAS  Google Scholar 

  95. Benjamini, Y. & Hochberg, Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J. R. Stat. Soc. 57, 289–300 (1995).

    MathSciNet  MATH  Google Scholar 

Download references

Acknowledgements

This work was supported by grants of the German Research Foundation (DFG: SFBTRR219 to R.K.), CRU344-4288578857858 and CRU5011-445703531 to R.K., by two grants from the European Research Council (ERC-StG 677448, ERC-CoG 101043403), a grant from the Else Kroener Fresenius Foundation (EKFS), the Dutch Kidney Foundation (DKF), TASKFORCE EP1805, the NWO VIDI 09150172010072 and a grant from the Leducq Foundation, all to R.K., a grant from the Germany Society of Internal Medicine (DGIM) and a grant from the European Research Council (ERC-StG 101040726) to C.K., and a grant from the IZKF Faculty of Medicine at the RWTH Aachen University to I.G.C. This work was also supported by the BMBF eMed Consortia Fibromap (to I.G.C., R.K. and R.K.S.). J.S.-R., R.O.R.F. and R.T.L. are supported by Informatics for Life funded by the Klaus Tschira Foundation. D.Schapiro. and F.W. are supported by the German Federal Ministry of Education and Research (BMBF 01ZZ2004). P.B. is supported by the German Research Foundation (DFG, Project IDs 322900939, 454024652). 10X Genomics supported the project via the Visium challenge (to C.K.) and provided free data generation for Visium data of 8 specimens. We thank J. Cool and J. Hilton for supporting data submission to cellxgene Data Portal; E. Sapena Ventura and G. Yordanova for their support; J. Chew and S. Williams (10X Genomics) for their support; A. Valdeolivas, J. Lanzer, A. Gabor, P. Badia i Mompel and D. Dimitrov for fruitful discussions that shaped the computational analysis. We also acknowledge the support from E. Siera, U. Förster, V. Künstler, N. Graff, P. Cappelmann, L. Tenten, I.-K. Härthe and L. Palm. R.O.R.F. and J.S.-R. acknowledge support from the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant no INST 37/935-1 FUGG.

Author information

Authors and Affiliations

Authors

Contributions

C.K. and R.K. planned and designed the study. R.K. obtained funding for the study. R.O.R.F. and Z.L. designed the data and computational analysis, advised by J.S.-R. and I.G.C. C.K., R.O.R.F., Z.L., S.H., R.T.L., M.Halder., M.T.H., D.S., S.M., G.S., K.H., J.T., I.G.C., F.P., N.K., T.S., R.J.A.V., J.A., P.B., J.S.-R. and RK analysed and interpreted the data. C.K., R.O.R.F., Z.L., M.H., I.G.C., J.S.-R. and R.K. wrote the manuscript and organized the figures. R.K.S., K.L., H.M., L.S. and S.S. edited the manuscript and advised on data interpretation. L.W.V.L., A.V., K.K., A.K., J.G., M.M., K.L. and H.M. organized patient tissue collection and biobanking and consent from the patients. N.K. established the PDGFRβ cell line from human hearts. S.Z., X.Z. and C.K. carried out the overexpression studies. C.K. together with X.Z., X.L. and Y.X. carried out all other single-cell and imaging experiments. Image analysis was performed by F.W. and D.Schapiro. R.M.H. and E.M.J.B. validated and sequenced the single-cell libraries. I.G.C., M.C., J.S.N. and Z.L. carried out snATAC-seq data analysis, its integration with snRNA-seq and spatial transcriptomics data, cell sub-state definitions and gene-regulatory networks. R.R., M.H., J.T. and J.S.-R. carried out the analysis and integration of spatial gene expression and snRNA-Seq data. R.O.R.F. designed the multi-omics patient comparison, advised by J.S.-R. C.K., H.M. and R.K. initiated the study. All authors read and approved the manuscript.

Corresponding authors

Correspondence to Julio Saez-Rodriguez or Rafael Kramann.

Ethics declarations

Competing interests

The authors declare no competing interests directly related to this work. The authors however disclose unrelated funding, honorariums and ownership as follows: R.K. has grants from Travere Therapeutics, Galapagos, Chugai and Novo Nordisk and is a consultant for Bayer, Pfizer, Novo Nordisk and Gruenenthal. J.S.-R. reports funding from GSK and Sanofi and fees from Travere Therapeutics and Astex Therapeutics. K.L. has grants from Novartis and Amgen and receives consulting fees from Medtronic and Implicit Biosciences. I.G.C. has a grant from Illumina. L.S. has received grants from Bayer, Boehringer Ingelheim and Nattopharma, is consultant for ImmunoDiagnostics Systems (IDS) and is a shareholder in Coagulation Profile. L.W.V.L. reports consultancy fees to UMCU from Abbott, Medtronic, Vifor and Novartis.

Peer review

Peer review information

Nature thanks Eva van Rooij and the other, anonymous, reviewers for their contribution to the peer review of this work. Peer review reports are available.

Additional information

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

Extended data figures and tables

Extended Data Fig. 1 Computational workflow and snRNA-Seq data analysis.

a, Schematic of the computational workflow for the main analyses of snRNA-seq, snATAC-seq, and spatial transcriptomics data. b, UMAP embedding of snRNA-seq data from all samples for patients, regions, and clusters. c, UMAP embedding of snRNA-seq data showing G2-M-phase cell-cycle score (left) and S-phase score (right). d, Barplots showing cell-type proportions of snRNA-seq data across all samples. Colours indicate different major cell types. e, Comparison of the generated snRNA-seq atlas to a previously published human heart cell atlas (HCA) for molecular profiles (left) and cell-type proportion (right). Left panel shows the adjusted P-value of the overlap of top gene markers of each cell-type between the different data sets (hypergeometric test). Right panel shows the Pearson correlation between the median proportion of each shared cell-type of the reference atlas and our control, border zone, and remote zone samples. f, UMAP embedding and annotation of an external dataset of three ischaemic zone samples following MI. g, Comparison of the generated snRNA-seq atlas to the external ischaemic data for molecular profiles (left) and cell-type proportion (right). Left panel shows the adjusted P-value of the overlap of top gene markers of each cell-type between the different data sets (hypergeometric test). Right panel shows the Pearson correlation between the median proportion of each shared cell-type of the external data set and our ischaemic samples.

Extended Data Fig. 2 snATAC-seq and spatial transcriptomics data analysis.

a, UMAP showing snATAC-seq data for all patients, regions, and clusters. b, Validation of snATAC-seq cell type annotation using snRNA-seq data. Left: UMAP showing the predicted labels for snATAC-seq data. Right:  heatmap showing evaluation results by adjusted rand index (ARI). c, snATAC-seq cell-type proportions of all samples. d, Spearman correlations of cell-type compositions of samples estimated from snRNA-seq and and snATAC-Seq. Box-Whisker plots showing median and IQR (n = 3 for BZ, n = 4 for control, n = 4 for RZ, n = 5 for FZ, and n = 9 for IZ) e, Cell type detection for snRNA-seq and snATAC-seq data across samples. The colour represents in which modality was cell type detected: snRNA-seq and snATAC-seq, snRNA-seq only, or not recovered. f, Transcription factor (TF) binding and TF target expression per major cell type based on the snATAC-seq and snRNA-seq data. g, Overrepresented spatially variable biological processes (myocardium related, immune related and fibrosis related) across regions. Each cell contains the mean adjusted P-value of a hypergeometric test across all spatial transcriptomics samples belonging to each region. h, Spearman correlations of cell-type compositions of samples estimated from spatial transcriptomics and snRNA-seq (left) or snATAC-seq (right). Box-Whisker plots showing median and IQR (n = 3 for BZ, n = 4 for control, n = 4 for RZ, n = 5 for FZ, and n = 9 for IZ). i, GWAS SNP enrichment score across major cell types. j, Visualization of GWAS12 ll SNP enrichment (left ventricular ejection fraction) in spatial transcriptomics data. In d,h, each spot is a patient sample (n = 3 for borderzone (BZ), n = 4 for controls, n = 6 for fibrotic zone (FZ), n = 9 for ischaemic zone (IZ), n = 5 for remote zone (RZ)). Data are represented as boxplots where the middle line is the median, the lower and upper hinges correspond to the first and third quartiles, the upper whisker extends from the hinge to the largest value no further than 1.5 × IQR from the hinge (where IQR is the inter-quartile range) and the lower whisker extends from the hinge to the smallest value at most 1.5 × IQR of the hinge, while data beyond the end of the whiskers are outlying points that are plotted individually. In a-c the number of spots of the bottom panels correspond to the barplots in the upper panel.

Extended Data Fig. 3 Cell-type niches and spatially contextualized views analysis.

a-b, UMAP embedding of spatial transcriptomics (ST) spots based on major cell-type compositions coloured by the patient (a) or region (b). c, Cell-type niche compositions across patient sample. d, UMAP embedding of ST spots based on major cell-type compositions colored by the compositions of cycling cells, pericytes, adipocytes, endothelial cells and myeloid cells and their respective marker genes (RYR2 cardiomyocytes, PDGFRA fibroblasts, MKI67 cycling cells, ABCC9 pericytes, FASN adipocytes, VWF endothelial cells, IL7R myeloid cells, and MYH11 vSMCs). e, Standardized mean PROGENy pathway activities across different niches. f, H&E staining and cell2location cell-type abundance estimations of cardiomyocytes, fibroblasts and myeloid cells. Delineated areas represent myogenic or fibroblast/myeloid cell enriched tissue areas. g, Median standardized importances (> 0) of cell-type abundances in the prediction of other cell types within the immediate neighbourhood (upper part) and the extended neighbourhood (effective radius of 15 spots) (lower part) inferred from spatially contextualized models. Cell-type abundances of the immediate (upper panels) and extended neighbourhood of cardiomyocytes, fibroblasts and myeloid cells. h, Visualization examples of the  dependencies between the abundance of endothelial cells and the abundance of pericytes in the immediate neighborhood, and vascular smooth muscle cells in the extended neighbourhood (effective radius of 15 spots) visualized on three tissues.

Extended Data Fig. 4 Spatial pathway and cell type dependency analysis.

a, Median standardized importances (> 0) of cell-type abundances within the spot and in the extended neighbourhood (effective radius = 15), and PROGENy pathway activities within a spot, and in the immediate or extended neighbourhood on the prediction of pathway activities inferred from spatially contextualized models. b, Standardized importances of cardiomyocyte abundances in the extended neighbourhood (radius of 15 spots) in predicting hypoxia pathway activity. Comparison was performed between two sample groups: control, border and remote zones samples (n = 10), and fibrotic and ischaemic samples (n = 15). Two-sided Wilcoxon rank sum test. c–e, Visualization examples of pathway to pathway and pathway to cell-type dependencies in spatial transcriptomics. (c) p53 and PI3K spatial pathway distribution. (d) JAK-STAT and NFkB activities’ dependencies to myeloid and lymphoid cells’ abundance. (e) Cardiomyocyte to WNT and hypoxia pathway dependency. f,  Hierarchical clustering of the pseudo-bulk transcriptional profile of spatial transcriptomics slides of patient samples. Colour represents the three major sample groups. g, Comparisons of patient mean major cell-type proportions (inferred from snRNA, snATAC-seq data, and spatial transcriptomics) between myogenic, ischaemic and fibrotic groups. Two-sided Wilcoxon rank sum test. Each dot represents a patient sample (n=13 for myogenic group, n=9 for ischaemic group, n=5 for fibrotic group). h, Pairwise comparison between patient groups of the standardized importances of cell-type abundances within the same spot (intra), and in the immediate (juxta) or extended neighborhood (para) to predict other cell types’ abundances (two-sided Wilcoxon rank sum test). Myogenic vs ischaemic (left), myogenic vs fibrotic (middle), fibrotic vs ischaemic (right). * = adj. P-value <= 0.15. i, Comparison of standardized importances of lymphoid cells’ abundances in the immediate neighborhood (juxta) to predict myeloid cells between myogenic, ischaemic, and fibrotic groups (two-sided Wilcoxon rank sum test). Visualization of lymphoid and myeloid cells in an ischaemic (up) and control (down) sample. j, Comparison of standardized importances of pericytes’ abundances to predict cardiomyocytes’ abundances within the same spot (two-sided Wilcoxon rank sum test). Each dot represents a sample (n = 14 for myogenic group, n=9 for ischaemic group, n = 5 for fibrotic group). Spatial distributions of pericyte and cardiomyocytes abundances estimated from deconvolution in a myogenic sample (left) and a fibrotic sample (right). Arrows in the fibrotic enriched sample show a strong colocalization event. k, Comparison of the compositions of cell-type niches between myogenic, fibrotic, and ischaemic groups (Kruskal-Wallis test, line denotes an adj. P-value < 0.1). Comparison of the proportions of cell-type niches 4, 5, 7, and 8 between groups (adj. P-value estimated from two-sided Wilcoxon rank sum test) and visualization example. Each dot represents a patient sample (n for myogenic group = 13, n for ischaemic group = 7, n for fibrotic group = 5). l, Comparison of the proportions of cell-type niche 9 between groups (two-sided Wilcoxon rank sum test). Each dot represents a patient sample (n for myogenic group = 13, n for ischaemic group = 7, n for fibrotic group = 5). m, Pairwise comparison of the compositions of the cell-type niches between control, border and remote zone samples (two-sided Wilcoxon rank sum test). * = reflect changes with adj. P-value <  0.1, n for CTRL = 4, n for RZ = 5, n for BZ = 3. In b, g, i, j, k, l data are represented as boxplots where the middle line is the median, the lower and upper hinges correspond to the first and third quartiles, the upper whisker extends from the hinge to the largest value no further than 1.5 × IQR from the hinge (where IQR is the inter-quartile range) and the lower whisker extends from the hinge to the smallest value at most 1.5 × IQR of the hinge, while data beyond the end of the whiskers are outlying points that are plotted individually.

Extended Data Fig. 5 Characterization of molecular niches.

a, UMAP embedding of ST spots based on gene expression coloured by patients (left) and regions (right). b, Gene expression of RYR2, ABCC9, MYH11, PDGFRa, CD8A and IL7R. Colours refer to gene-weighted kernel density as estimated by using R package Nebulosa. c, Bar plots visualizing molecular niches proportion per patient. d, Standardized mean PROGENy pathway activities across different molecular niches. e, Visualization of molecular niches in control, IZ, and FZ samples. f, Comparison between patient groups of the proportions of molecular niches 5, 6, and 1 (adj. P-value from a two-sided Wilcoxon rank sum test). Each dot represents a patient sample (n for myogenic group = 13, n for ischaemic group = 9, n for fibrotic group = 5). g, Comparison between control (CTRL), border zone (BZ) and remote zone (RZ) of the proportions of molecular niches 1, 2, 3, and 4 (adj. P-value, two-sided Wilcoxon rank sum test, Box-Whisker plots showing median and IQR. Maximum and minimum values as described previously). Each dot represents a patient sample (n for BZ = 3, n for CTRL = 4, n for RZ group = 4). h, Top five differentially upregulated genes between spots belonging to molecular niches enriched with cardiomyocytes. Summary area under the curve (AUC) is used as a size effect of comparing the expression of one molecular niche against the rest. * = reflect FDR < 0.05 and summary AUC > 0.55 (n for mol. niche 0 = 30,058, n for mol. niche 1 = 19,958, n for mol. niche 3 = 7,360). Spatial distribution of the expression of MYLK3 in a control slide (upper) and a border zone (lower)  i, Comparison between patient groups of the proportions of molecular niche 2 (adj. P-value, two-sided Wilcoxon rank sum test, Box-Whisker plots showing median and IQR. Maximum and minimum values as described previously). Each dot represents a patient sample (n for myogenic group = 13, n for ischaemic group = 9, n for fibrotic group = 5). j, Same as (i) for niche 3. In f, g, i, j data are represented as boxplots where the middle line is the median, the lower and upper hinges correspond to the first and third quartiles, the upper whisker extends from the hinge to the largest value no further than 1.5 × IQR from the hinge (where IQR is the inter-quartile range) and the lower whisker extends from the hinge to the smallest value at most 1.5 × IQR of the hinge, while data beyond the end of the whiskers are outlying points that are plotted individually.

Extended Data Fig. 6 Characterization of cardiomyocyte state clusters.

a, UMAP embedding of the integrated snRNA-seq and snATAC-seq data coloured by patients, regions, modalities, and clusters (resolution= 1). b, Dot plot showing the top 10 upregulated genes for each CM-state. c, State-specific pseudo bulk ATAC-seq profiles showing distinct chromatin accessibility for NPPB between different states. d, In situ RNA hybridization of NPPB, ANKRD1 and TNNT2 on human cardiac control tissue (upper panel). Note that only a single NPPB/ANKRD1 positive cardiomyocyte could be detected. (below) Representative image of the same in-situ staining of a human myocardial infarction sample. Note the NPPB and ANKRD1 expression in several cardiomyocytes (arrows). Scale bars: 25 µm (upper) + 50 µm (lower). e, Spearman correlation between the spatial expression of Ion-channel transport and transmural-ion channels related genes to the expression of the marker genes (state scores) of vCM1-3 (two-sided Wilcoxon rank sum test). Each dot represents a visium slide (n = 28). f, Differential gene expression of ion channel related genes in vCM1 contrasted to vCM3 cardiomyocytes (two-sided Wilcoxon rank sum test, area under the curve (AUC) is shown as size effect). Colours refer to the gene set membership of each gene. g, Visualization of vCM1 gene marker expression mapping (state scores) and the expression of RYR2 and ATP2B4 in a border zone sample. h, Comparison of vCM1 and vCM2 cell proportion between patient groups. P-values were calculated using Wilcoxon Rank Sum test (unpaired, two-sided). Each dot represents a sample (n = 13 for myogenic group, n = 7 for ischaemic group, and n = 4 for fibrotic group). i, Distribution of vCM3 marker gene expression (state-score) across molecular niche 1, 2, and 4. Each dot represents a spatial transcriptomics spot belonging to a molecular niche across samples (n = 30,058 for niche 1, n = 19,958 for niche 2, niche 4 = 7,360). Two-sided Wilcoxon rank sum test, adj. p-values (niche 1 vs niche 4 = 0, niche 2 vs niche 4 = 2e-09, niche 1 vs niche 2 = 0). j, Visualisation of the expression of ANKRD1 and NPPB, the spatial distribution of TGFβ, hypoxia, p53 and PI3K signaling activities, the spatial distribution of the expression of vCM-states marker genes (state score) for a border zone sample and the distribution of molecular niches 1,2, and 4. Box-Whisker plots showing median and IQR. Maximum and minimum values as described previously. In e, h, i data are represented as boxplots where the middle line is the median, the lower and upper hinges correspond to the first and third quartiles, the upper whisker extends from the hinge to the largest value no further than 1.5 × IQR from the hinge (where IQR is the inter-quartile range) and the lower whisker extends from the hinge to the smallest value at most 1.5 × IQR of the hinge, while data beyond the end of the whiskers are outlying points that are plotted individually. In a-c the number of spots of the bottom panels correspond to the barplots in the upper panel.

Extended Data Fig. 7 Cardiomyocyte pseudotime and gene-regulatory network analysis.

a, Diffusion map embedding of pseudotime cardiomyocyte clusters vCM1, vCM2 and vCM3. b, Computational workflow for building gene regulatory network (upper) and schematic of linking TF to target genes through peak-to-gene links and predicted TF binding sites (lower). c, Heatmap showing TF binding activity and gene expression along the pseudotime trajectory. d, Heatmap of gene expression across pseudotime vCM1-vCM3. e, Heatmap showing the correlation between TF binding activity and gene expression along the pseudotime trajectory from vCM1 to vCM3. Colors on the top refer to pseudotime point where the TF showed the highest binding activity. Clustering analysis identified three modules. f, ANKRD1 and NPPB CM-stress marker gene expression along pseudotime vCM1-vCM3. g, Peak-to-gene links showing that JDP2 regulates TGFB2. Each loop represents a putative link between TGFB2 and a peak. Loop height represents the significance of the correlation and dash line represents threshold of significance (P = 0.05). ATAC-seq tracks were generated from pseudo-bulk chromatin profiles of vCM1, vCM2, and vCM3. Binding sites of JDP2 are highlighted. h, Line plots showing TF activity and expression after z-score normalization (y-axis) over pseudotime (x-axis) for JDP2 (left), its corresponding target gene expression after z-score normalization (y-axis) over pseudotime (x-axis) for TGFB2 (middle), and visualization of all target genes in the BZ sample (right). i, Comparison of standardized importances of myeloid abundances within the spot and fibroblast abundances in the local neighbourhood (radius of 5 spots) to predict vCM3 between patient groups samples. Each dot represents a sample (n = 9 for myogenic group, n = 7 for ischaemic group, n = 4 for fibrotic group). (lower panel) Spatial distribution of the state score of vCM3 and myeloid cell abundances in a RZ slide with histological evidence of a scar. (Two-sided Wilcoxon rank sum test, Box-Whisker plots showing median and IQR. Maximum and minimum values as described previously.) Data are represented as boxplots where the middle line is the median, the lower and upper hinges correspond to the first and third quartiles, the upper whisker extends from the hinge to the largest value no further than 1.5 × IQR from the hinge (where IQR is the inter-quartile range) and the lower whisker extends from the hinge to the smallest value at most 1.5 × IQR of the hinge, while data beyond the end of the whiskers are outlying points that are plotted individually. In a-c the number of spots of the bottom panels correspond to the barplots in the upper panel.

Extended Data Fig. 8 Endothelial cell heterogeneity.

a. UMAP embedding of integrated snRNA-seq and snATAC-seq data coloured by patients, regions, modalities and clusters (resolution= 0.9). b, Gene expression of NRG3, FLT4, SEMA3G and ACKR1. Colours refer to gene-weighted kernel density as estimated by using R package Nebulosa. c, Sub-cluster-specific pseudo bulk ATAC-seq tracks showing the chromatin accessibility of the genes from (b). d, Marker dot plot showing the DEGs for each endothelial cells state and subtype. e, Gene expression of SEM3G and POSTN. Colors refer to gene-weighted kernel density as estimated by using R package Nebulosa. Immunofluorescence of SEMA3G and ACTA2 (upper image). In situ mRNA (RNAscope) for POSTN and PECAM1 (magenta). Arrows highlight endocardial endothelial cells. Scale bars: 75 µm (upper) and 10 µm (lower). f,  Endothelial cell state compositions across all patient samples. g, Comparison between patient groups of the venous endothelial cell proportion (two-sided Wilcoxon rank sum test). Each dot represents a patient sample (n for myogenic group = 13, n for ischaemic group = 9, n for fibrotic group = 5). h, Mean Pearson correlation between the abundance of major cell-types and endothelial cell-state scores across all spatial transcriptomics slides. i, Visualization of the spatial distribution of the abundances of pericytes and the state score of capillary endothelial cells sample. Arrows point at colocalization events. j, Visualization of the spatial distribution of molecular niches, the state score of capillary endothelial cells and the abundance of pericytes and cardiomyocytes. Arrows point at locations where molecular niche 10 is present together with high abundances of cardiomyocytes and pericytes. k, Endothelial cell state score distributions in all spots belonging to the molecular niche 10 across all slides (n = 1,874, P-value = 3.19e-298, obtained from a two sided Wilcoxon signed-rank test).  l, Mean signalling pathway activities in the molecular niche 10 across different patient groups (adj. P-value of two-sided Wilcoxon Rank Sum test). Each dot represents a slide (n = 14 for myogenic group, n = 9 for ischaemic group, n = 5 for fibrotic group). m, Overrepresented hallmark pathways in the differentially upregulated genes of molecular niche 10 across different patient groups (hypergeometric tests, adj. P-values). In g, k, l data are represented as boxplots where the middle line is the median, the lower and upper hinges correspond to the first and third quartiles, the upper whisker extends from the hinge to the largest value no further than 1.5 × IQR from the hinge (where IQR is the inter-quartile range) and the lower whisker extends from the hinge to the smallest value at most 1.5 × IQR of the hinge, while data beyond the end of the whiskers are outlying points that are plotted individually. In a–c the number of spots of the bottom panels correspond to the barplots in the upper panel.

Extended Data Fig. 9 Fibroblast heterogeneity and PDGFRb lineage tracing in murine MI.

a. UMAP embedding of integrated snRNA-seq and snATAC-seq data coloured by patients, regions, modalities and clusters (resolution = 0.9). b, Marker dot plot showing the DEGs for each fibroblast state. c, Box plots showing module score of NABA-collagen and NABA-core-matrisome score per fibroblast state. P-values were calculated using Wilcoxon rank sum (unpaired and two sided). d, Expression of COL1A1, ACTA2 and FN1 by RNA qPCR after RUNX1 overexpression with and without TGFβ compared to empty vector (EV) (n = 6). One-way ANOVA followed by Bonferroni correction. Error bars = S.D. e, In situ hybridization (RNAscope) on human myocardial tissue of SCARA5 and COL15a1 and quantification and comparison of SCARA5+/POSTN+ cells vs. POSTN+/COL1A1+ cells in human heart failure tissues (n = 7). Mann-Whitney test. Error bars = S.D. f,  Visualization of SCARA5, POSTN, COL1A1 and FN1 on spatial transcriptomics slides from IZ and FZ human MI samples. g, Cell-type state compositions across patient sample. h, Comparison of Fib3 cell proportion between patient groups. P-values were calculated using Wilcoxon Rank Sum test (unpaired, two-sided). Each dot represents a sample (n = 13 for myogenic group, n = 7 for ischaemic group, and n = 4 for fibrotic group). i, Time course of lineage tracing experiment using PDGFRβCreER-tdTomato mice. j, Results of echocardiographic measurements EF (in %), and global longitudinal strain (GLS, in %), for each of the time points. One-way-ANOVA. n=4 day 0, n=3 day 4, n=4 day 7, n=4 day 14. Error bars = S.D.  k, Quality measurements of single-cell RNA-Seq data from mouse MI experiment. l, UMAP representation of time points, cluster and gene expression (POSTN, SCARA5, COL1A1 and FN1) from mouse MI experiment. m, Confusion matrix comparing predicted fibroblasts states and obtained clusters for mouse dataset. n, UMAP representation of species cross-annotation of fibroblast clusters. o, Cell proportion of mouse fibroblast state Fib1 (SCARA5+) and Fib2 (POSTN+) per MI time-point. p, Box plots showing NABA collagens scores and NABA core matrisome score per fibroblast state (mouse Fib1 vs. Fib1) per time point. P-values were calculated using Wilcoxon rank sum (unpaired and two sided) (Sham: n = 578 for Fib1 and n = 685 for Fib2; Day 4: n = 1688 for Fib1 and n = 2498 for Fib2; Day 7: n = 203 for Fib1 and n = 330 for Fib2; Day 14: n = 2232 for Fib1 and n = 7684 for Fib2). q, Gene set enrichment analysis per human fibroblast-cell state. Scale bar: 10 µm. In c, d, h, k, p data are represented as boxplots where the middle line is the median, the lower and upper hinges correspond to the first and third quartiles, the upper whisker extends from the hinge to the largest value no further than 1.5 × IQR from the hinge (where IQR is the inter-quartile range) and the lower whisker extends from the hinge to the smallest value at most 1.5 × IQR of the hinge, while data beyond the end of the whiskers are outlying points that are plotted individually.

Extended Data Fig. 10 Gene regulatory network analysis of human cardiac fibroblasts.

a, Pseudotime heatmap showing gene expression (left) and TF binding activity (right) along the trajectory Fib1 to Fib2 (myofibroblasts). b, Pseudotime heatmap showing highly variable genes along the trajectory. c, Heatmap showing the correlation between TF binding activity and gene expression each TF from (a) and gene from (b). Each column represents a TF and each row represents a TF. Colours on the top refer to pseudotime labels for each TF estimated based on binding activity. d, Upper: line plots showing TF activity and expression after z-score normalization (y-axis) over pseudotime (x-axis) for GLI2 and their corresponding target gene expression after z-score normalization (y-axis) over pseudotime (x-axis) for TGFB1. Lower: visualization of target gene expression of KLF4 and GLI2, and TGFB pathway activity in the BZ sample. e, Line plots showing TF activity and expression after z-score normalization (y-axis) over pseudotime (x-axis) for RUNX2 and its corresponding target gene expression after z-score normalization (y-axis) over pseudotime (x-axis) for COL1A1. Visualization of RUNX2 target genes in the BZ sample (right).

Extended Data Fig. 11 Myeloid cell heterogeneity and spatial mapping in human myocardial infarction.

a, UMAP embedding of integrated snRNA-seq and snATAC-seq data coloured by patient contribution, region, modality and clusters (resolution= 1). b, Marker dot plot showing the DEGs for each fibroblast state. c, Gene expression of MRC1, ITGAX, FLT3 and CD163. Colors refer to gene-weighted kernel density as estimated by using R package Nebulosa. d, Bar plots visualizing myeloid cell sub-population proportion per patient. e, Gene expression of MRC1, ITGAX, SPP1, CCL18, FLT3, ZBTB46 in an external dataset of ischaemic patients. Colors refer to gene-weighted kernel density as estimated by using R package Nebulosa. f, Comparison of myeloid cell proportion between patient groups. P-values were calculated using Wilcoxon Rank Sum test (unpaired, two-sided) (n = 13 for myogenic, n = 8 for ischaemic, and n = 5 for fibrotic group). g, In situ hybridization of CCR2, SPP1 and TREM2 on human myocardial infarction section IZ. Arrows point to phagocytotic vacuoles in the magnification. h, In situ hybridization of SPP1, POSTN and CD163 on human myocardial infarction section IZ. Note that SPP1+ macrophages colocalized in distinct tissue regions and appeared to have an enlarged cell-size. i, In situ hybridization quantification of cell type proportion per tissue section. n=7 patient tissues. Adjusted P-values were calculated using Wilcoxon signed-rank test. j, Median standardized importances (>0) of the cell-state scores of myeloid cells within the spot in the prediction of other myeloid cell-state scores in spatial transcriptomics. k, Mean pearson correlation between myeloid cell’s state scores across all spatial transcriptomics slides. Data in f, i are represented as boxplots where the middle line is the median, the lower and upper hinges correspond to the first and third quartiles, the upper whisker extends from the hinge to the largest value no further than 1.5 × IQR from the hinge (where IQR is the inter-quartile range) and the lower whisker extends from the hinge to the smallest value at most 1.5 × IQR of the hinge, while data beyond the end of the whiskers are outlying points that are plotted individually.

Extended Data Fig. 12 Myeloid cell spatial modelling and cell-cell communication.

a, Median standardized importances (>0) of the cell-state scores of myeloid cells (colocalization) in the prediction of fibroblast cell-state scores in spatial transcriptomics. b, Spatial mapping of cell-type niches in a control and an ischaemic sample. Arrows in the IZ sample show niche 4 and how it surrounds niche 5. c, Distribution of myofibroblasts and SPP1+ macrophages marker gene expression (state-score) across cell-type niche 4 and 5. Each dot represents a spatial transcriptomics spot belonging to a molecular niche across samples (n = 12,427 for niche 4, n = 4,375 for niche 5) (two-sided Wilcoxon rank sum test, adj P-value = 6.02e-133 for myofibroblast, adj P-value = 0 for SPP1+ monocytes). Data are represented as boxplots where the middle line is the median, the lower and upper hinges correspond to the first and third quartiles, the upper whisker extends from the hinge to the largest value no further than 1.5 × IQR from the hinge (where IQR is the inter-quartile range) and the lower whisker extends from the hinge to the smallest value at most 1.5 × IQR of the hinge, while data beyond the end of the whiskers are outlying points that are plotted individually. d, Cell-cell communication network representing number of ligand-receptor interactions (edge richness), expression of ligand-receptor pairs (LR scores; colour gradients) and cell centrality (Pagerank; node size) as estimated in snRNA-seq of myogenic, ischaemic and fibrotic samples. e, Sankey plots summarizing the top 50 ligand-receptor interactions for selected source and target cells and contrast. These ligand-receptor pairs are selected by absolute value of the difference in LRscore, as provided by CellPhoneDB method implemented in Liana. f. Sankey plot summarizing top 50 TGFbeta ligand-receptor interactions from Fib2 to SPP1+ Mac. cells. g, In situ hybridization mRNA (RNAscope) staining of CD163 (myeloid), POSTN (myofibroblast) and SPP1 on human cardiac MI tissue (crop-out in Fig. 6n). Arrows indicate CD163+SPP1+ macrophages near myofibroblasts. Scale: 25 µm.

Supplementary information

Reporting Summary

Supplementary Figure 1

Sampling location, macroscopic and microscopic images. a, Schematic of human cardiac tissue sampling locations. LV-free wall n= 20. Apex n=11. b, Macroscopic images of myocardial tissue specimens. Arrows point to sampled infarcted tissue areas. c, Histological images (trichrome staining) of cardiac tissues included in this study. Arrows highlight the fibrotic tissue areas identified by trichrome staining. Scale bar sizes as indicated.

Supplementary Figure 2

Quality measurements of different data modalities and samples. a, snRNA-Seq QC. Included are the number of nucleus per sample, number of genes (per nucleus), and number of unique molecular identifiers (UMIs) (per nucleus). Samples from different tissue regions are represented in different colors. b, snATAC-Sq QC. Included are the number of nucleus per sample, number of fragments per sample (per nucleus) and transcription start site (TSS) enrichment scores (per nucleus). c, Spatial transcriptomics QC. Included are the number of spots per sample, the number of UMIs (per spot), the number of quantified nucleus (per spot) and the number of genes (per spot). d, Mean “pathway expression” of BioCarta’s “Death Pathway” and Reactome’s “Regulated Necrosis Pathway” per unfiltered slide, and number of recovered nucleus per sample in snRNA-Seq data separated by sample regions. Right panels show for an ischaemic sample the normalized mean gene expression of the previously mentioned pathways per spot and the number of UMIs. n = 4 for control, n = 5 for RZ, n = 3 for BZ, n = 6 for FZ, and n = 9 for IZ. e-g, Spatial transcriptomic histology (hematoxylin and eosin stainings, H&Es) and UMI distribution across the slide for controls (e), remote zones (f) and border zone (g). Box-Whisker plots showing median and inter-quartile range (IQR). Raw data for (a)-(c) are provide in the Supplementary Tables 2-4. Scale bars 8mm.

Supplementary Figure 3

Spatial quality controls. a-b, Spatial transcriptomic histology (hematoxylin and eosin stainings, H&Es) and UMI distribution across the slide for ischaemic zones (a) and fibrotic zone (b) samples. Scale bars 8mm.

Supplementary Figure 4

FACS-gating strategy. a, Gating strategy on DAPI-positive nuclei isolated from frozen human heart tissue.

Supplementary Table 1

Clinical data for all samples used in this study.

Supplementary Table 2

Meta information and quality check for snRNA-seq samples.

Supplementary Table 3

Meta information and quality check for snATAC-seq samples.

Supplementary Table 4

Meta information for visium samples

Supplementary Table 5

Meta information for imaging samples

Supplementary Table 6

Number of nuclei per major cell type for snRNA-seq and snATAC-seq.

Supplementary Table 7

Spatially variable genes

Supplementary Table 8

Cell-type-specific gene expression estimated from snRNA-seq

Supplementary Table 9

Cell-type-specific gene chromatin accessibility estimated from snATAC-seq

Supplementary Table 10

Marker genes for cardiomyocytes states

Supplementary Table 11

Gene regulatory network for cardiomyocytes

Supplementary Table 12

Marker genes for endothelial cell subtypes

Supplementary Table 13

Marker genes for fibroblasts states

Supplementary Table 14

Gene regulatory network for fibroblasts

Supplementary Table 15

Marker genes for myeloid cell subtypes

Supplementary Table 16

Statistical results of the sub-clusters for cardiomyocytes, endothelial cells, fibroblasts, and myeloid cells.

Supplementary Table 17

Ion-channel related gene sets

Supplementary Table 18

qPCR Primer

Peer Review File

Source data

Rights and permissions

Springer Nature or its licensor holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Kuppe, C., Ramirez Flores, R.O., Li, Z. et al. Spatial multi-omic map of human myocardial infarction. Nature 608, 766–777 (2022). https://doi.org/10.1038/s41586-022-05060-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • Issue Date:

  • DOI: https://doi.org/10.1038/s41586-022-05060-x

This article is cited by

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