Functional coordination of non-myocytes plays a key role in adult zebrafish heart regeneration

Abstract

Cardiac regeneration occurs primarily through proliferation of existing cardiomyocytes, but also involves complex interactions between distinct cardiac cell types including non-cardiomyocytes (non-CMs). However, the subpopulations, distinguishing molecular features, cellular functions, and intercellular interactions of non-CMs in heart regeneration remain largely unexplored. Using the LIGER algorithm, we assemble an atlas of cell states from 61,977 individual non-CM scRNA-seq profiles isolated at multiple time points during regeneration. This analysis reveals extensive non-CM cell diversity, including multiple macrophage (MC), fibroblast (FB), and endothelial cell (EC) subpopulations with unique spatiotemporal distributions, and suggests an important role for MC in inducing the activated FB and EC subpopulations. Indeed, pharmacological perturbation of MC function compromises the induction of the unique FB and EC subpopulations. Furthermore, we developed computational algorithm Topologizer to map the topological relationships and dynamic transitions between functional states. We uncover dynamic transitions between MC functional states and identify factors involved in mRNA processing and transcriptional regulation associated with the transition. Together, our single-cell transcriptomic analysis of non-CMs during cardiac regeneration provides a blueprint for interrogating the molecular and cellular basis of this process.

Keywords: Topologizer; heart regeneration; non-myocytes; scRNA-Seq; zebrafish.

Figure 1. Single‐cell transcriptomics reveals heterogeneity of zebrafish cardiac non‐myocytes

1. A

   Experimental workflow of non‐CM isolation from zebrafish hearts and scRNA‐seq (10× Genomics).

2. B

   scRNA‐seq data of adult zebrafish cardiac non‐CM visualized on tSNE and colored by cell types. EC, endothelial cells; FB, fibroblasts; MC, macrophages; Mes, resident mesenchymal cells; T/NK/B, T/NK/B cells; Neutro, neutrophils; Eryth, erythrocytes; Throm, thrombocytes.

3. C

   Violin plots showing expression of canonical markers for each cell type. GO analysis (DAVID) of upregulated genes in each population was performed and representative GO terms were listed on the right.

4. D–F

   Non‐CMs expressing canonical EC markers in panel b were zoom‐in analyzed with LIGER. (D) Cells visualized on tSNE and colored by cell types. eEC, endocardial EC; lEC, lymphatic EC; cEC, coronary EC; Mural, mural cells. (E) Pie chart showing contribution of each cell type. (F) Expression of newly identified markers of each EC subpopulation shown on tSNE.

5. G

   Pie chart showing non‐CM composition (erythrocytes and thrombocytes excluded).

6. H

   Dotplot showing expression of top eight positive markers identified for each non‐CM population.

Figure 2. Transcriptome dynamics of macrophages during zebrafish heart regeneration

1. A, B

   Joint analysis of non‐CM scRNA‐seq data from uninjured hearts and hearts at 2, 7, and 14 dpi with LIGER. (A) All non‐CM visualized on tSNE. (B top) Non‐CM from each time point visualized on the tSNE embedding in (A). Data were down‐sampled to the same cell number at each time point. (B bottom) Pie charts showing contribution of different cell types at each time point.

2. C–F

   Zoom‐in analysis of macrophages identified in (A). (C) tSNE plot colored by MC subpopulations. (D) Expression levels of representative markers of each MC subpopulation color‐coded and mapped to tSNE embeddings in (C) (top), and corresponding representative GO terms for each subpopulation (bottom). (E) Dynamics of MC subpopulations during heart regeneration as shown by the proportion of each subpopulation in total non‐CM or total MC. (F) Expression of tnfa and csf3b in MC at each time point.

3. G

   Fluorescent in situ hybridization for cd74a, tnfa, and ctsc at 7 dpi, respectively. mpeg1 was used to label all MC cells. The blue boxed region is highlighted in the zoom‐in image to the right. Scale bar = 50 μm. A stands for atrium, V stands for ventricle, and OFT stands for outflow tract. White dashed lines outline the heart and the yellow dashed lines indicate approximate resection plane.

Figure 3. Transcriptome dynamics of fibroblasts and endothelial cells during zebrafish heart regeneration

1. A–C

   Zoom‐in analysis of fibroblasts identified in Fig 2. (A) tSNE plot colored by FB subpopulations. (B) Dot plot showing expression of ECM genes in FB subpopulations. (C) Bar plot showing relative proportion of FB subpopulations in FB at each time point.

2. D, E

   Expression of fn1a, col12a1a (D) and tagln (E) in FB at each time point.

3. F

   In situ hybridization showing temporal spatial expression patterns of postnb. The white boxed region is shown in zoom‐in images at the bottom. Scale bar = 25 μm. White dashed lines outline the heart, and the yellow dashed lines indicate approximate resection plane. A stands for atrium, V stands for ventricle, and OFT stands for outflow tract.

4. G

   Expression of tnfrsf1a and tgfbr2a in FB at each time point.

5. H

   Fluorescent in situ hybridization for postnb, tnfα, and tnfrsf1a in the injury area at 7 dpf, respectively. The white boxed regions are shown in their respective zoom‐in images to the right. White dashed lines outline cardiac apex, and the yellow dashed lines indicate approximate resection plane. Scale bar = 25 μm.

6. I

   tSNE plot colored by eEC subpopulations. Bar plot showing relative proportion of eEC subpopulations in eEC at each time point.

7. J

   Violin plots showing expression of cxcl18b, sele, atf3, and fosl1a in each eEC subpopulation. In panels (E) and (I), “un” stands for uninjured.

8. K

   Fluorescent in situ hybridization for fosl1a at 7 dpi. The white boxed region is shown in a higher magnification image to the right for the injury area. cdh5 was used to label all EC cells and the red boxed region is shown in a higher magnification image to highlight a pan‐EC cdh5 expression in green. White dashed lines outline the heart and the yellow dashed lines indicate approximate resection plane. Scale bar = 25 μm. A stands for atrium, V stands for ventricle and OFT stands for outflow tract.

Source data are available online for this figure.

Figure 4. Loss of kit function impairs heart regeneration

1. A

   Wild‐type or kit mutant stained with AFOG at 14 days and 30 dpi. Sham‐operated zebrafish hearts serve as control. Scale bar = 100 μm.

2. B

   Quantification of scar area at 30 dpi. N = 10 hearts. Data are presented as Mean ± SEM. P‐value calculated with two‐tailed Student's t test. **P < 0.01.

3. C

   PCNA/Mef2 double staining showing the proliferating CMs in wild‐type and mutant hearts at 7 dpi. The white boxed regions are shown in zoom‐in images to their right. Blue dashed lines indicate approximate resection plane. White arrows point to the PCNA/Mef2 double positive nuclei. Scale bar = 50 μm.

4. D

   Quantification of % PCNA⁺ CMs (Mef2⁺) in (C). N = 10 hearts. Data are presented as Mean ± SEM. P‐value calculated with two‐tailed Student's t test. **P < 0.01.

5. E, F

   Analysis of bulk RNA‐seq data of freshly isolated CMs at different time points during heart regeneration. (E) Two representative gene sets from GSEA analysis, which are enriched in DEGs between mutant and wild‐type CMs at 2 dpi. (F), Heatmap showing expression levels of cell cycle regulators in CMs from wild‐type and kit mutant hearts.

6. G, H

   Joint analysis of kit mutant non‐CM scRNA‐seq data from 0, 2, 7, and 14 dpi with LIGER. (G) All non‐CM visualized on tSNE. (H) Bar plots showing comparison of non‐CM cell type contribution between wild‐type and kit mutant at each time point.

7. I

   Immunostaining for MC marker IB4 in WT (left) and mutant hearts (middle) at 7 dpi. Quantification of IB4⁺ cell number per unit area (right). The white boxed regions are shown in their respective zoom‐in images to the right. White dashed lines indicate approximate resection plane. Scale bar = 50 μm. N = 6 hearts.

Data are presented as Mean ± SEM. P‐value is calculated with two‐tailed Student's t test. *P < 0.05. In panels (F) and (H), “un” stands for uninjured.

Figure 5. Loss of kit function altered transcriptome dynamics of cardiac non‐CM during zebrafish heart regeneration

1. A–G

   Joint analysis of WT and kit mutant MC. (A) tSNE plot colored by WT and kit mutant MC subpopulations. (B) Representative GO terms of genes upregulated in MC1 (WT, red bars) or KMC1 (kit mutant, blue bars) from the uninjured hearts, respectively. Dotted lines indicate P = 0.05. (C) Bar plots showing proportion of each WT and kit mutant MC subpopulation in total non‐CM at each time point. (D, F) Violin plots showing the expression of il1b in MC1 (D), and ctsd expression in MC (F). (E, G) Expression of il1b (E) and ctsd (G) in WT and kit mutant hearts determined by qRT‐PCR.

2. H–N

   Joint analysis of WT and kit mutant FB. (H) tSNE plot colored by WT and kit mutant FB subpopulations. (I) Bar plots showing proportion of each WT and kit mutant FB subpopulation in total non‐CM at each time point. (J, L) Violin plots showing expression of col12a1a in FB3 (J) and mif expression in FB (L). (K, M) Expression of col12a1a (K) and mif (M) in WT and kit mutant hearts determined by qRT‐PCR. (N) Bar plots showing relative proportion of each WT and kit mutant eEC subpopulation in eEC at each time point.

3. O

   Violin plots showing expression of junba in eEC.

4. P

   Expression of junba in WT and kit mutant hearts determined by qRT‐PCR.

5. Q

   Representative image of hearts from 7 dpi tcf21:nucGFP transgenic fish treated by PBS liposome, Clodronate liposome, and Dex, respectively. The red boxes mark the peripheral area, and the white boxes mark the injury area and are shown in their respective zoom‐in images to the right. White dashed lines indicate approximate resection plane. Scale bar = 50 μm.

6. R

   Quantification of the number of tcf21:nucGFP‐positive cells in the boxed areas in Q. N = 5 hearts.

Data are presented as Mean ± SEM. Two‐tailed Student's t test. ***P < 0.001. Data information: (E, G, K, M, P). Mean ± SEM shown. N = 3. P‐value calculated with two‐tailed Student's t test. *P < 0.05, **P < 0.01, ***P < 0.001. In panels (C–G, I–P), “un” stands for uninjured. Clod lipo stands for clodronate liposomes. Dex for dexamethasone. Source data are available online for this figure.

Figure 6. Combined Topologizer and RNA velocity analyses identified molecular regulatory mechanism associated with MC functional state transition during heart regeneration

1. Schematic of Topologizer.

2. Topological structure of WT MC subpopulations.

3. Vector field of RNA velocity projected onto the Topologizer structure for MCs collected from uninjured WT hearts or WT hearts at 2, 7, and 14 dpi. Arrows indicate the direction and “speed” of the velocity at each node.

4. Quantification of arrow length in (C). P‐value is calculated with one‐tailed Wilcoxon rank sum test. *P < 0.05.

5. Transcriptional similarity between MC1 and MC2 subtypes from the uninjured (un) hearts or hearts at 2, 7, and 14 dpi, respectively.

6. Hierarchical tree showing the relationship between the MC subtypes and the myeloid cells from kidney marrow.

7. Three gene clusters identified associated with MC2–MC1 transition.

8. Feature significant gene ontology (GO) terms (adjusted P‐value < 0.05) with representative genes. The number of genes is shown in parentheses. Dotted lines indicate P = 0.05.

9. Topological structure of WT FB subpopulations.

10. Vector field of RNA velocity projected onto the Topologizer structure for FBs collected from uninjured WT hearts or WT hearts at 2, 7, and 14 dpi. Arrows indicate the direction and “speed” of the velocity at each node.

11. Topological structure of WT eEC subpopulations.

12. Vector field of RNA velocity projected onto the Topologizer structure for eECs collected from uninjured WT hearts or WT hearts at 2, 7, and 14 dpi. Arrows indicate the direction and “speed” of the velocity at each node.