Comparative single-cell profiling reveals distinct cardiac resident macrophages essential for zebrafish heart regeneration

Zebrafish exhibit a robust ability to regenerate their hearts following injury, and the immune system plays a key role in this process. We previously showed that delaying macrophage recruitment by clodronate liposome (CL) treatment compromises neutrophil resolution and heart regeneration, even when the infiltrating macrophage number was restored within the first-week post injury (Lai et al., 2017). Here, we examined the molecular mechanisms underlying the cardiac repair of regenerative PBS-control hearts vs. non-regenerative CL-treated hearts. Bulk transcriptomic analyses revealed that CL-treated hearts exhibited disrupted inflammatory resolution and energy metabolism during cardiac repair. Temporal single-cell profiling of inflammatory cells in regenerative vs. non-regenerative conditions further identified heterogenous macrophages and neutrophils with distinct infiltration dynamics, gene expression, and cellular crosstalk. Among them, two residential macrophage subpopulations were enriched in regenerative hearts and barely recovered in non-regenerative hearts. Early CL treatment at 8 days or even 1 month before cryoinjury led to the depletion of resident macrophages without affecting the circulating macrophage recruitment to the injured area. Strikingly, these resident macrophage-deficient zebrafish still exhibited compromised neovascularization and scar resolution. Our results characterized the inflammatory cells of the zebrafish injured hearts and identified key resident macrophage subpopulations prerequisite for successful heart regeneration.

1d_CL treatment led to disrupted inflammatory resolution, ROS homeostasis, and energy 23 metabolism during cardiac repair. Comparative single-cell RNAseq profiling of inflammatory cells 24 from regenerative vs. non-regenerative hearts further identified heterogeneous macrophages and 25 neutrophils, showing alternative activation and cellular crosstalk leading to neutrophil retention and 26 chronic inflammation. Among macrophages, two residential subpopulations (hbaa + Mac 2 and 27 timp4.3 + Mac 3) were enriched only in regenerative hearts and barely recovered after -1d_CL 28 treatment. To deplete the resident macrophage without delaying the circulating macrophage 29 recruitment, we established the resident macrophage-deficient model by administrating CL earlier at 30 8 days (-8d_CL) before cryoinjury. Strikingly, resident macrophage-deficient zebrafish still exhibited 31 defects in revascularization, cardiomyocyte survival, debris clearance, and ECM remodeling/scar 32 resolution without functional compensation from the circulating/monocyte-derived macrophages. 33 Our results characterized the diverse function and interaction between inflammatory cells and 34 identified unique resident macrophages prerequisite for zebrafish heart regeneration. 35

36
Introduction 37 Heart failure is a major cause of morbidity and mortality, in part due to the inability of the 38 human heart to replenish lost cardiomyocytes following myocardial infarction (MI). Unlike adult 39 mice and humans, many vertebrates, including certain fish and amphibians, are capable of 40 endogenous heart regeneration throughout life. As an example, zebrafish (Danio rerio) display a 41 distinct ability to regenerate their heart following injury. However, this ability is not shared by 42 another teleost, the medaka (Oryzias latipes) (1,2). Even mammals can regenerate their heart during 43 embryonic and neonatal stages, despite this capacity is quickly lost postnatally (3, 4). Comparative 44 studies between neonatal and adult mice (5, 6), and between phylogenetically related species such as 45 zebrafish and medaka (1,2) have suggested that the capacity for regeneration does not solely rely on 46 genetic makeup, environmental conditions (e.g., hypoxia), or evolutionary complexity; instead, the 47 type and extent of the immune responses to cardiac injury seem to be a major difference between 48 these regenerative and non-regenerative models (1,7), and may largely influence the recovery post-49 experimental MI, as well as clinical prognosis (8). 50 In our earlier study of reciprocal analyses in zebrafish and medaka, we observed delayed and 51 reduced macrophage recruitment in medaka compared to zebrafish following cardiac injury. 52 Furthermore, delaying macrophage recruitment in zebrafish by intraperitoneal (IP) injection of 53 clodronate liposome (CL) 1 day prior to cryoinjury compromised neovascularization, neutrophil 54 clearance, cardiomyocyte proliferation, and scar resolution, even though the number of infiltrating 55 macrophages recovered to the control levels in the first-week post injury (-1d_CL, macrophage-56 delayed model hereafter). These previous results indicate that late macrophages in -1d_CL-treated 57 zebrafish were different in their identity and regenerative potential, presenting a great opportunity to 58 further compare and learn the molecular properties of regeneration-associated macrophages. Recent 59 studies in zebrafish also identified novel macrophage subpopulations by gene-specific reporters and 60 their functions in cardiac repair and regeneration, hinting at a heterogeneous spectrum of cardiac 61 macrophages (9-13). However, beyond these observations, the overall heterogeneity and function of 62 these inflammatory cells during cardiac regeneration remain unclear. 63 In the present study, we perform comparative bulk and single-cell transcriptomic profiling to 64 investigate the global influence of macrophage pre-depletion and comprehensively analyze the 65 heterogeneity, dynamic, and function of both macrophages and neutrophils in regenerative zebrafish 66 hearts (-1d_PBS) vs. non-regenerative macrophage-delayed hearts (-1d_CL). Bulk RNAseq analysis 67 indicated prolonged and unresolved inflammatory response and mis-regulated energy metabolism in 68 -1d_CL-treated zebrafish until 3 weeks post cardiac injury, while cardiomyocyte replenishment and 69 scar resolution took place extensively in regenerative PBS-treated control hearts. Single-cell analyses 70 further revealed diverse macrophage subpopulations with potential functions in phagocytosis, 71 neutrophil recruitment, reactive oxygen species (ROS) homeostasis, angiogenesis, extracellular 72 matrix (ECM) remodeling, and inflammatory regulation during the first week post cardiac injury. 73 Comparative analyses between regenerative and non-regenerative hearts led to the identification of 74 unique cardiac resident macrophage subpopulations expressing timp4.3 and hemoglobin genes that 75 potentially function in ECM remodeling, inflammatory resolution, and ROS homeostasis. Pre-76 depleting these resident macrophages a week or even a month prior to cardiac injury significantly 77 impaired heart regeneration without affecting macrophage recruitment from circulation, suggesting 78 that these resident macrophages determine the regenerative capacity and cannot be replaced by 79 circulating/monocyte-derived macrophages. Altogether, these results unravel the heterogeneity and 80 function of inflammatory cells during cardiac repair and highlight the indispensable role of cardiac 81 resident macrophages in zebrafish heart regeneration. 82 83 Results 84 Delayed macrophage recruitment results in prolonged expression of genes related to 85 inflammation and disrupted energy metabolism during cardiac repair 86 In zebrafish, macrophage pre-depletion at 1-day prior to cardiac injury (-1d_CL) delayed 87 macrophage recruitment and compromised heart regeneration, even though the overall macrophage 88 numbers were restored within the first week (1). To investigate the global transcriptomic changes 89 under these regenerative (PBS-treated, normal macrophage recruitment) and non-regenerative (-90 1d_CL-treated, delayed macrophage recruitment) conditions, we isolated zebrafish hearts at 7 and 21 91 day-post cryoinjury (dpci), corresponding to the time when cardiomyocytes proliferate and replace 92 the scar tissue during heart regeneration, and subjected the hearts to bulk RNAseq analyses ( Figure  93 1A). We first plotted our data from 7 and 21 dpci against published data points from 6 hours post 94 cryoinjury (hpci) to 5 dpci in a Principal Component Analysis (PCA, Figure 1B). Transcriptomes of 95 the PBS-treated samples at 7 and 21 dpci nicely fit into the trajectory between 5 dpci and untouched 96 hearts, suggesting that transcriptomic changes in control hearts coincide with the recovery toward the 97 naïve state (purple dots and trail, Figure 1B). However, the transcriptome of non-regenerative CL-98 treated hearts at 21 dpci (CL21d) was in proximity to that of PBS-control hearts at 7 dpci (PBS7d), 99 suggesting a delayed transcriptomic response in the non-regenerative CL-treated hearts (green dots, 100 Figure 1B). 101 We next performed hierarchical clustering of all differentially expressed genes (DEGs) that 102 behave differently under regenerative PBS vs. non-regenerative -1d_CL conditions, ( Figure 1C and 103 Figure 1-source data 1). These DEGs were grouped into 22 clusters with similar expression 104 dynamics corresponding to time points and treatments. Among regeneration-associated gene clusters, 105 we found genes associated with ROS homeostasis and injury repair were up-regulated only in PBS-106 hearts at 7 dpci (Cluster 2, C2), including serpine1, havcr2, tnfaip6, and hmox1a. These genes may 107 contribute to heart regeneration through ROS regulation and inflammatory modulation (14-19). At 21 108 dpci, DEGs encoding various ribosomal subunits were up-regulated only in PBS-control hearts (C9-109 C12), suggesting active production of building blocks for replenishing lost tissues ( Figure 1C). 110 Interestingly, some DEGs function in various metabolic processes (C19 and C20) were active in 111 uninjured hearts, downregulated at 7 dpci, and reactivated only in PBS-treated but not CL-treated 112 hearts at 21 dpci ( Figure 1C). This observation corresponds nicely to the metabolic switches of CMs 113 during cardiac regeneration. CMs in adult animals adopt oxidative metabolism after differentiation 114 and maturation to meet the high energy demands from constant beating (20). Upon myocardial 115 infarction, mature CMs switch back to use glucose (glycolysis) instead of fatty acid (oxidative 116 phosphorylation) as the main substrate for energy (21,22). Interestingly, this metabolic switch was 117 also observed during zebrafish regeneration, when pyruvate metabolism and glycolysis are beneficial 118 for CM dedifferentiation and proliferation (23).

119
On the other hand, the majority of the DEGs associated with non-regenerative -1d_CL condition 120 are involved in immune-related processes, including damage-associated patterns, inflammatory 121 cytokines, phagocytosis, and apoptotic cell death (C3-C6). These DEGs were generally expressed at 122 7 dpci in both PBS and CL-treated hearts, but their activation was prolonged and even intensified in 123 CL-treated hearts at 21 dpci ( Figure 1C). Prolonged inflammation may prevent tissue repair and 124 regeneration (24). Inflammation also induces continuous neutrophil migration and infiltration and 125 prevents their clearance by apoptosis (25). Especially in C5 and C6, DEGs involved in "Immune 126 response", "Cytokine-cytokine receptor interaction", and "Apoptosis" were overexpressed at 21 dpci 127 only in CL-treated hearts suggesting that cardiac cells (including CMs) continuously undergo 128 programmed cell death potentially suffering from ROS stress (ferroptosis) and inflammatory 129 microenvironment (pyroptosis) (26, 27). Continuous cell death subsequently triggers the prolonged 130 activation of C-type lectin receptor singling pathway in active phagocytes (28, 29). As a result, 131 proteolysis of engulfed debris and cell adhesion molecules were activated under the same conditions 132 ( Figure 1C). Altogether, the results from hierarchical clustering of DEGs and associated GO analysis 133 suggest that delayed macrophage recruitment disrupts inflammatory resolution, ROS homeostasis, 134 and energy metabolism during cardiac repair. The full list of DEGs in each cluster was summarized 135 in Figure 1-source data 1.

136
To identify the canonical pathways and potential upstream regulators associated with aberrant 137 regeneration, these DEGs were further analyzed by Ingenuity Pathway Analysis (IPA) ( Extravasation Signaling" pathway, we found several integrin genes usually expressed on the 144 leukocyte plasma membrane, including itga4, itgal and itgb2, which are involved in leukocyte 145 enrolling on endothelial cells and transmigration (30, 31). In addition, we found several matrix 146 metalloproteinases (MMPs) including mmp9, mmp13, and mmp25, which might be involved in ECM 147 remodeling and leukocyte recruitment during inflammation ( Figure 1-figure supplement 1C) (32-34).

148
These results support our previous observation of continuous neutrophil infiltration and retention in 149 the CL-treated hearts even until 30 dpci (1). Lastly, among genes of the "NO and ROS production in 150 Macrophages" pathway, we found continuous activation of DAMP/PAMP receptor tlr2, neutrophil 151 cytosolic factors (ncf1, ncf2, and ncf4), myeloid cell-lineage committed gene spi1 and its 152 downstream target ptpn6, in addition to the macrophage differentiation marker irf8 in CL-treated 153 hearts at 21 dpci ( Figure 1-figure supplement 1D). Collectively, these results indicate that 154 macrophages may play roles in regulating ROS homeostasis, immune cell dynamics, and 155 inflammation resolution, as these processes were mis-regulated under -1d_CL treatments and 156 associated with impaired heart regeneration.

158
Single-cell analyses reveal the heterogeneous landscape and dynamic changes of inflammatory 159 cells during cardiac repair 160 Since macrophage properties may be altered upon -1d_CL treatment and thus fail to support 161 heart regeneration by resolving neutrophil infiltration and inflammation, we analyzed and compared 162 the potential identity and function of these inflammatory cells by single-cell transcriptomic profiling 163 ( Figure 2A). Adopting the same macrophage-delayed model (-1d_CL treatment), double transgenic 164 zebrafish Tg(mpx:EGFP;mpeg1:mCherry) were IP injected with CL or PBS at 1 day before 165 cryoinjury, and the EGFP-expressing neutrophils and mCherry-expressing macrophages were 166 isolated by fluorescence-activated cell sorting (FACS) from uninjured hearts, as well as regenerative 167 PBS-control and non-regenerative -1d_CL treated hearts at 1, 3 and 7 dpci ( Figure 2A and rapidly, and divergent cell composition and numbers were observed in PBS vs. CL-treated hearts 175 over time. While macrophage numbers increased after cardiac injury until 7 dpci, neutrophils peaked 176 at 3 dpci and gradually decreased between 3 and 7 dpci in PBS-control hearts, corresponding to the 177 inflammatory resolution phase previously described (Figure 2-figure supplement 1A) (10). In CL-178 treated hearts, macrophages progressively increased, similar to control hearts, but neutrophil 179 numbers became much higher than controls at both 3 and 7 dpci (  Inflammatory cells were subjected to droplet-based high-throughput single-cell RNA 189 sequencing (scRNAseq). To visualize the dataset, sequencing reads were mapped to the zebrafish 190 genome, assigned to each cell, and then processed by the Uniform Manifold Approximation and 191 Projection (UMAP) for dimension reduction and unbiased clustering using Seurat R package ( Figure  192 2B and Figure  cells has been previously shown to express mpeg1 and observed in mpeg1:mCherry fish (40). Minor 202 mpeg1 -/mpxclusters might come from contamination during FACS, even though stringent gating 203 strategies were applied, and will not be further analyzed in this study (41, 42).

204
Besides common lineage markers, heterogeneous macrophage subsets exhibited cluster-205 enriched/specific gene expression ( Figure 2C). Unlike macrophages, neutrophils were classified into 206 only two populations Neu 1 and Neu 2, which both seem to be mature and express granular genes 207 and integrins ( Figure 2C) and Temporal cell proportion analyses identified specific resident macrophage subsets associated 217 with heart regeneration 218 To dissect the dynamic changes of these inflammatory cell clusters under regenerative PBS vs. 219 non-regenerative -1d_CL conditions, we generated split UMAPs for each time point and condition 220 ( Figure 3A). In uninjured heart/naïve state, Mac 2 and 3 represented the major resident macrophage 221 clusters followed by Mac 1, 4 and 8 ( Figure 3B). The proliferating macrophage cluster Mac 5 was 222 also observed in the uninjured hearts, corresponding to those residing nearby the epicardial tissue, 223 suggesting that some resident macrophages may self-renew through local proliferation ( quickly at 1 dpci and gradually reduced back to a steady state at 7 dpci, while Mac 2 and 3 expanded 226 substantially over the first-week post injury ( Figure 3B). On the contrary, we noticed a dramatic 227 reduction of Mac 2 and retention of Mac 1, 4, 5 and 6 in the -1d_CL-treated hearts over time after 228 injury ( Figure 3B). Lastly, the minor resident populations Mac 8 and 9 were diminished after cardiac 229 injury and barely recovered in both conditions ( Figure 3B). We then calculated the cluster 230 contribution toward regenerative vs. non-regenerative conditions in a cell proportion analysis ( Figure  231 3C and D). Coincidently, macrophage clusters enriched in regenerative conditions are the major 232 resident clusters Mac 2 and 3, which were either dramatically decreased or barely recovered in CL-233 treated hearts ( Figure 3C). To examine the key enriched genes of major resident clusters Mac 2 and On the other hand, Neu 1 and Neu 2 were two heterogeneous neutrophil clusters actively 244 recruited to hearts after cardiac injury. We observed a decrease of both clusters from 3 to 7 dpci in 245 the PBS-control hearts, while they retained in the CL-treated hearts by 7 dpci ( Figure 3A). The 246 proportion of Neu 1 was slightly higher than Neu 2 under regenerative conditions throughout the first 247 week post cardiac injury ( Figure 3D). These results delineate the dynamic changes of each 248 inflammatory cell cluster and identified regeneration-associated resident macrophages preferentially 249 enriched in PBS-control hearts which might play pivotal roles in cardiac regeneration.

251
Alternative activation of inflammatory cells during cardiac repair under regenerative vs. 252 macrophage-delayed conditions 253 To investigate the potential function of these heterogeneous inflammatory cells, we extracted 254 and analyzed the overall DEGs of each cluster (Table 1, cluster-enriched genes, full-listed in Figure  255 2-source data 2) and their differential changes toward regenerative conditions (PBS-enriched) vs. 256 non-regenerative conditions (CL-enriched) (Figure 4, condition-enriched genes, full-listed in Figure  257 4-source data 1). We aimed to identify cluster-enriched/specific markers of each macrophage 258 subpopulation, and their function associated with regenerative vs. fibrotic repair by condition- cluster-enriched and condition-enriched genes suggest that circulation-recruited/monocyte-derived 299 macrophages may exhibit distinct functional polarization toward continuous neutrophil recruitment 300 and inflammation upon -1d_CL pre-depletion, reflecting the complex nature of macrophage 301 polarization and functions in cardiac repair.

302
Alternative genes activation under regenerative vs. non-regenerative conditions was also 303 observed in neutrophil clusters, especially for Neu 1 (enriched under the regenerative condition in 304 Figure 3D). Neu 1 upregulated retinoic acid receptor raraa, necroptosis genes fth1a, caspa, and 305 hmgb1a, and inflammatory modulation genes cul1b and spop under regenerative condition (Figure 4, 306 D14), which corresponds nicely with the fact that programmed cell death and phagocytosis are the 307 main mechanisms of neutrophil clearance and inflammatory resolution(72-74). On the other hand, 308 both Neu 1 and Neu 2 expressed genes related to typical neutrophil functions under non-regenerative 309 condition, including vesicle transport, phagocytosis, energy metabolisms, regulation of proteolysis, debris clearance, inflammatory modulation, and turn on programmed cell death for its own clearance. 320 On the contrary, Neu 2 enriched in the macrophage-delayed model and functions in inflammatory 321 propagation and recruitment of more inflammatory cells, nicely explained the continuous neutrophil 322 recruitment and retention that we observed in the CL-treated heart. 323

324
Cellular crosstalk analysis indicates that resident macrophages mediate ECM remodeling and 325 phagocytic clearance of neutrophils 326 Since our previous study and current findings indicate that neutrophils retain in the injury 327 associated with unresolved inflammation when macrophage properties change in -1d_CL-treated To reveal the potential progression from resident macrophages from naïve state to 378 heterogeneous subpopulations post cardiac injury, we performed cell trajectory analysis using 379 Monocle3 ( Figure 6). We set the root from Mac 3, the main resident population identified in the 380 uninjured hearts, which highly expresses progenitor/lineage markers spi1b, coro1a, and irf8. 381 Macrophage clusters and neutrophil clusters were both subjected to the analysis due to the MN 382 clusters constitute of both cell characteristics. As a result, macrophages occupied the pseudotime 383 point 0 to pseudotime point 11 and arranged on the radial trajectories from Mac 3 while Cluster MN 384 was the intermediary cluster along the trajectory to neutrophil population ( Figure 6A). On the other 385 hand, Neu 2 and Neu 1 occupied the pseudotime point 13 to pseudotime point 27 ( Figure 6A). Four 386 radial trajectories/routes were depicted, suggesting that macrophages were highly plastic and 387 simultaneously polarized into heterogeneous subtypes with specialized functions during cardiac 388 repair. The results also suggested resident macrophages may play a sentinel role by detecting tissue 389 damage and recruiting neutrophils to the injured site as previously reported (98). In combination with 390 the profiling results, the resident Mac 3 mainly functions in immune system process and cell 391 differentiation, suggesting they might involve in monocyte recruitment and differentiation into 392 macrophages (Table 1). Route I represents transition from Mac 3 to Mac 1, which is mainly involved 393 in iron homeostasis, angiogenesis and anti-inflammatory activities with vigorous endocytosis 394 activities under the regenerative condition (Table 1 and Figure 6A). Route II represents transition 395 from Mac 3 to pro-inflammatory clusters MN and show the regenerative-associated polarization 396 ( Figure 6A) that MN clusters were constituted of more macrophage-like cells in regenerative hearts 397 (purple route II) and more neutrophil-like cells in non-regenerative hearts (green route II). Route III 398 represents transition from Mac 3 to proliferating Mac 5 via Mac 4, 6 and 7 for responding to 399 oxidative stress, oxidative phosphorylation and neutrophil infiltration ( Figure 6A). Interestingly, the 400 last route IV reveals transition from Mac 3 to Mac 2, and some of them shifted back with their gene 401 signature similar to Mac 3, suggesting potential transition of these clusters during heart regeneration 402 ( Figure 6A). Coincidently, Mac 2 is the most expanded cluster in regenerative condition ( Figure 3B), 403 while this plastic transition was absent in CL7d ( Figure 6B). Next, we identified the DEG changes 404 dynamically over the pseudotime of cardiac repair and analyzed the GO BP and KEGG pathways in 405 which they are involved ( Figure 6C). At early pseudotime points 0-4, genes related to phagocytosis 406 (marco, mrc1b, and havcr1), complement system (c1qc and c1qb), and antigen sensing/presentation 407 pathogen receptor (cd74a, mhc2dab mfap4) were highly enriched. These expression pattern fit nicely 408 with the role of resident macrophages in cardiac homeostasis in removing damaged cells/components 409 and serving as sentinel cells for antigen presentation to other immune cells ( Figure 6C) (99, 100). 410 Later at pseudotime points 5-9, macrophages enhanced their phagocytosis and lysosomal capacity for 411 debris clearance in response to DAMPs in order to prevent further tissue damages ( Figure 6C). 412 Interestingly, a small set of genes peak in the middle pseudotime points 10-14, which correspond to 413 diverse function in regeneration (grn1 and grna), oxygen transport (habb1 and hbba1.1), cell 414 proliferation (top2a, mki67, and cdk1), and inflammation/ribosome function (tnfb, tlr5b, and rpls) 415 ( Figure 6C). These results represent a diverse polarization of different macrophage subpopulations 416 observed in respective routes ( Figure 6A). Lastly, neutrophil-specific genes (lect2l, tcnbb and scpp8) 417 were highly expressed after pseudotime point 14. Glycolysis related genes associated with 418 inflammatory activities in neutrophils were consistently identified in Neu 2, while npsn and lyz were 419 enriched at the latest pseudotime during neutrophil progression in Neu 1 ( Figure 6C). Taken together, 420 these results depict the main trajectories/routes of how macrophage might polarize and progress to 421 gain different function during homeostasis and post cardiac injury, and highlight the regenerative 422 route constitute of mainly resident clusters Mac 2 and 3.

424
Depletion of resident macrophage compromises heart regeneration 425 Based on the single-cell profiling results, regeneration-associated macrophages were mainly 426 resident macrophages Mac 2 and Mac 3, which were substantially enriched in regenerative PBS-427 control hearts compared with non-regenerative -1d_CL-treated hearts. To test the functional 428 significance of these resident macrophages without disrupting the circulation/monocyte-derived 429 macrophage recruitment, we perform CL depletion earlier at 8 days prior to cardiac injury (-8d_CL, we observed from single-cell profiling ( Figure 3A). Next, we performed cryoinjury on these fish 443 after -8d_CL treatment and examined the macrophage content at 1-day post injury (1 dpci; Figure 7-444 figure supplement 1B). Surprisingly, the proportion of mCherry + cells was even higher in -8d_CL 445 group (1.14%) compared to PBS-controls (0.83%). Since scRNAseq results showed specific loss of 446 Mac 2 and Mac 3, we examined their cluster-specific gene hbaa1 and timp4.3 expression by qPCR to 447 test whether the Mac 2 and Mac 3 could recover from -8d_CL depletion. Consistent with our 448 scRNAseq results, the expression of both genes in mpeg1:mCherry + cells were significantly reduced 449 in -8d_CL hearts compared to PBS-controls (Figure 7-figure supplement 1C). In situ hybridization 450 also showed drastic reduction of timp4.3 + cells in -8d_CL hearts than in PBS-controls (Figure 7-451  figure supplement 1D). These results support that early CL administration (-8d_CL) depletes resident 452 macrophages without affecting overall macrophage infiltration after cardiac injury. Replenished 453 macrophages after depletion may potentially derive from circulating monocyte, but resident 454 macrophage Mac 2 and Mac 3 remained diminished without replenishment from those 455 mpeg1:mCherry-low monocyte-derived cells (Figure 7-figure supplement 1A). Of note, we observed 456 an overshoot of macrophage infiltration to the -8d_CL injured hearts, consistent with previously 457 published results showing that -1d_CL treatment actually led to more macrophages infiltrated in the 458 injured heart at 7 dpci (1). This observation suggests an intrinsic role of resident macrophages in 459 modulating inflammation and immune cell recruitment after cardiac injury. 460 To determine the functional requirement of resident macrophages in cardiac regeneration, we 461 then examined the key processes of successful heart regeneration, including revascularization, CM 462 proliferation and scar resolution in -8d_CL and PBS-control hearts ( Figure 7A). Consistent with our 463 previous study, we used Tg(fli1:EGFP;myl7:DsRed-NLS) fish to visualize the vascular endothelial 464 cells in green and the nuclei of CMs in red, and used Acid Fuchsin-Orange G (AFOG) staining to 465 reveal the scar size and composition in the injured hearts (1). Fast revascularization of the injured 466 area in the first-week post injury is essential to support zebrafish heart regeneration and was 467 compromised in -1d_CL hearts (1, 57). We examined revascularization by ex-vivo imaging of GFP + 468 endothelial cells in whole mount hearts and found that revascularization of the injured area in -469 8d_CL hearts were significantly decreased than PBS-controls at 7 dpci ( Figure 7A and Figure 7-470 source data 1). Neutrophils numbers were also higher in the injured area of -8d_CL hearts, 471 suggesting a similar neutrophil retention phenotype ( Figure 7B and Figure 7-source data 1). We next 472 examined and quantified proliferating CMs by EdU incorporation, since dedifferentiation and 473 proliferation of the existing CMs are the major source to replenish the lost myocardial tissue from 474 injury (101). Unlike previously reported in -1d_CL hearts, we did not observe a significant decrease 475 of CM proliferation in -8d_CL hearts ( Figure 7C and Figure 7-source data 1). Instead, we noticed 476 small and round CM nuclei with weaker DsRed signals in the border zone of -8d_CL hearts 477 compared to controls ( Figure 7C; I-IV). Correspondingly, the density of CM nuclei was significantly 478 lower in the border zone of the -8d_CL hearts ( Figure 7C). Furthermore, injured areas were also 479 larger in -8d_CL hearts ( Figure 7C). Together, these results suggest that less CMs survived from the 480 initial injury and were able to re-enter cell cycle and proliferate. To test this possibility, we 481 performed TUNEL assay to label damaged nuclei of dying cells. Indeed, we found significantly more 482 TUNEL positive cells within the injured area of -8d_CL hearts ( Figure 7D and Figure 7-source data 483 1). As a result, -8d_CL hearts exhibited larger/unresolved scar tissues composed of both collagen and 484 fibrin than -8d_PBS hearts at 30 dpci, reflecting compromised heart regeneration ( Figure 7E and 485 Figure 7-source data 1). 486 Since we observed a replenishment of macrophages after -8d_CL treatment (Figure 7-figure  487 supplement 1A), we further tested whether the regenerative capacity may be restored after longer 488 recovery and perform cryoinjury 1 month after CL depletion (-1m_CL, Figure 7-figure supplement 2 489 and Figure 7-source data 2). Strikingly, -1m_CL hearts still failed in regeneration, exhibiting 490 significant defects in revascularization, neutrophil retention, and scar resolution (Figure 7-figure  491 supplement 2). These results suggest that resident macrophages are prerequisite for heart 492 regeneration in modulating the revascularization, CM survival, and the resolution of inflammation 493 and fibrotic scars, which cannot be replenished or functionally replaced by circulation/monocyte-494 derived macrophages upon depletion. 495

496
The resident macrophage population Mac 2 expresses hmox1a, involved in heme clearance 497 To gain mechanistic insight into how resident macrophages participate in heart regeneration, 498 particularly in CM survival, we further characterized hmox1a, a differentially expressed gene in 499 resident macrophage Mac 2 ( Figure 8A). Heme released by damaged erythrocytes (hemoglobin) and 500 CMs (myoglobin) is a major source of ROS stress after tissue injury, often lead to secondary damage 501 and extended cell death (19). hmox1a encodes Heme oxygenase-1, which degrades heme into 502 cardioprotective and anti-inflammatory effectors such as biliverdin/bilirubin, carbon monoxide (CO), 503 and free iron/ferritin (51). Biliverdin/bilirubin and CO further exert antioxidant and anti-504 inflammatory activities, upregulate genes such as il10, and promote M2 macrophage polarization 505 (50). Heme oxygenase-1 may also regulate iron homeostasis and ferroptosis, the main type of cell 506 death after myocardial infarction (19, 50). We first validated the expression of hmox1a in both -507 8d_PBS and -8d_CL hearts at 7dpci by ISH, and observed the hmox1a + cells with macrophage-like 508 morphology. These hmox1a + cells were located in the injured area and around the epicardial layer in 509 the -8d_PBS-control hearts, and largely reduced with weaker expression in resident macrophage-510 deficient hearts ( Figure 8B). We further confirmed their expression in macrophage by 511 immunostaining against Hmox1 in the mpeg1:mCherry background ( Figure 8C). Corresponding to 512 the scRNAseq profiling results and the ISH results, a portion of mpeg1:mCherry + macrophages 513 indeed express Hmox1 in the -8d_PBS-control hearts, while double positive cells were significantly 514 less in the resident macrophage-deficient hearts ( Figure 8C). Since Hmox1 functions in metabolizing 515 heme, we examined heme accumulation in the injured hearts by O-dianisidine staining ( Figure 8D). 516 Indeed, heme released in the injured area was largely cleared in the -8d_PBS-control hearts at 7 dpci 517 as a part of debris clearance and microenvironment remodeling during cardiac repair. However, 518 excessive heme accumulated in the injured area of the resident macrophage-deficient hearts, 519 potentially affect the CMs survival as we observed previously ( Figure 8D and Figure 7C). These 520 results suggest that resident macrophages accelerate cardiac repair partly by clearing heme and 521 preserve the regenerative CMs from ROS stress and extended cell death. cardiac injury by CL-mediated predepletion compromised zebrafish heart regeneration in terms of 528 revascularization, neutrophil retention, CM proliferation and scar resolution, even though the 529 macrophage numbers gradually recovered before 7 dpci (1). These data indicated a certain degree of 530 functional divergence in the late infiltrating macrophages compared with the macrophages in the 531 PBS-control hearts. Dynamics of inflammatory cells were also investigated during cardiac repair in 532 mammals. Upon myocardial infarction, reperfusion is the standard practice in clinics, which salvages 533 some CMs from ischemic death and thus reduce the infarcted area. In the mouse ischemic-534 reperfusion (IR) model, the vessel occlusion is released 30 mins after ligation to allow blood flow 535 and immune cell trafficking, as compared to the permanent ligation (MI) model. Interestingly, 536 inflammatory cell dynamic in the ischemic-reperfusion model resembles the regenerative zebrafish, 537 when macrophage infiltration peak at an earlier time point, and neutrophils also resolve faster than in 538 the permanent MI model, even though the cell numbers might not be comparable due to the 539 differences in cell death and infarct size (100). Furthermore, macrophages from the IR model seem to 540 repolarize and become inflammation-resolving M2 type faster than those from the permanent MI 541 model, suggesting a progressive inflammation resolution (100). In the present study, we further 542 characterized the PBS vs. -1d_CL models in zebrafish and found the transcriptomic changes upon 543 macrophage delay. -1d_CL-treated hearts showed prolonged inflammation and dysregulated 544 metabolism even until 21 dpci, suggesting that macrophages may play important roles in modulating 545 neutrophil and inflammation resolution, which was further supported by the scRNAseq profiling of 546 these inflammatory cells (discussed below). However, whether and how macrophages may regulate 547 cardiac metabolism remain unclear. Metabolic shift in utilizing different metabolites for energy 548 production is a critical process during cardiac repair and regeneration (20-22). Compared with PBS-549 control hearts, we observed downregulated genes predicted to be involved in both glycolysis and 550 oxidative phosphorylation in CL-treated hearts from 7 to 21 dpci. The aberrant metabolism might 551 relate to macrophage phagocytosis function, since cardiac macrophages have been shown to preserve 552 the metabolic stability of CMs by actively clearing the dysfunctional mitochondria and other waste 553 via phagocytosis during homeostasis (99), while mitochondrial dysfunction leads to impaired 554 revascularization and fibrosis formation in mammals (102-104). Whether similar events occur during 555 cardiac repair and how macrophages may influence CM metabolism during dedifferentiation, 556 proliferation, and redifferentiation await future investigation. 557 558 ScRNAseq profiling revealed heterogeneous cardiac resident macrophages during steady-state 559 and repair/regeneration 560 Since macrophages and neutrophils seem to be the main players in the regenerative PBS vs. 561 non-regenerative -1d_CL conditions, we further profiled these inflammatory cells by scRNAseq. 562 During cell isolation, we noticed a portion of mpeg1:mCherry + /mpx:EGFP + cells in steady-state and 563 enriched in PBS-treated hearts post injury. These double positive cells were confirmed as 564 macrophages by cell sorting and Giemsa staining and were also described previously in the fin-fold 565 amputated larvae (36, 105). From the scRNAseq profiling, we only observed a corresponding 566 population of mpeg1 + /mpx + macrophages in Mac 2 at PBS7d. We speculate that EGFP proteins are 567 stable in sorted mpeg1:mCherry + /mpx:EGFP + macrophages and the mpx gene only needs to be turned 568 on after injury and the replenishment of Mac 2 in PBS-treated heart. Myeloperoxidase-expressing 569 macrophages were observed in human atherosclerotic plaque and increased during inflammation 570 (106). Among macrophage clusters, Mac 1, 2, 3, 4, 8 are the major subsets found in the 571 naïve/untouched hearts ( Figure 3B). The cluster enriched genes show similar functions to murine 572 cardiac resident macrophages including phagocytosis/homeostasis, angiogenesis, antigen 573 presentation, and sentinel function in monocytes and neutrophil chemotaxis (42, 107). However, the 574 cluster/lineage markers were not comparable between fish and murine cardiac resident macrophages. 575 The origin/lineage of fish cardiac macrophage is also difficult to determine without comparable 576 markers and lineage tracing tools as in rodents. Therefore, new genetic tools are needed for further 577 investigating these cardiac macrophage clusters. Nevertheless, judging by the dynamic change post 578 injury and their replenishment in CL-treated hearts, Mac 2, 3, and 8 might be originated from the 579 embryonic lineage, as they show minor or no replenishment after depletion and injury. On the other 580 hand, Mac 1 and 4 showed active recruitment to the injured heart and stronger ccr2 expression (Mac 581 1), suggesting they might be derived from monocytes ( Figure 3B). After cardiac injury, the dynamic 582 changes in numbers and gene expression among heterogeneous macrophage clusters reflect the swift 583 response to injury and high plasticity of polarization, making it difficult to define regeneration-584 associated macrophage properties. Thus, we combined multiple analyses, including cell proportional 585 analysis, GO analysis on cluster-enriched and condition-enriched DEGs, and crosstalk analysis in the 586 present study. Interestingly, we observed some macrophage clusters other than Mac 2 and 3 turn on 587 and off specific sets of genes under regenerative vs. non-regenerative conditions without affecting 588 their overall identity (in clustering, Figure 4). This phenomenon of gene alternative activation 589 suggests that functional polarization of circulation/monocyte-derived macrophages might be 590 regulated by resident macrophages. Interestingly, in addition to retention, neutrophils also exhibited 591 alternative property upon macrophage pre-depletion, with tissue repairing/regeneration associated 592 and self-clearance genes expression under regenerative condition and inflammatory exacerbation 593 genes under regenerative condition. These results highlight the importance of interactions between 594 neutrophils and both resident and recruited macrophages and how they coordinately modulate 595 inflammation and tissue repair. 596 597 Zebrafish cardiac resident macrophages are indispensable for heart regeneration 598 Our functional results reveal the essential roles of cardiac resident macrophages in zebrafish 599 heart regeneration, as they are required to support revascularization, CM survival and replenishment, 600 and eventually scar resolution. However, we did not observe cardiac resident macrophages directly 601 contributing to CM dedifferentiation and proliferation. CM repopulation is affected possibly due to 602 the reduced pool of regenerative cells ( Figure 7C), compromised revascularization (108), and the 603 abnormal ECM remodeling which leads scar deposition without affecting CM proliferation (61, 109). 604 From our results, we identified that resident Mac 2 expresses genes in regulating ROS homeostasis, 605 including hmox1a, which functions in heme degradation and alleviating ROS stress for CM survival. 606 In addition, we found that resident macrophages express ECM components and remodeling genes, 607 including fibronectin, collagen, ADAM metallopeptidases, and tissue inhibitor of metalloproteinase. 608 Dynamic remodeling of ECM is associated with immune cell behavior, fibrosis, tissue stiffness for 609 CM proliferation and migration, inflammation, together provide a pro-regenerative 610 microenvironment for cardiac repair. On the other hand, alternative activation of monocyte-derived 611 Mac 1 may contribute to CM proliferation directly as we found genes associated with angiogenesis 612 and cardiovascular development under regenerative condition. Particularly, vegf has been associated 613 with CM proliferation in both larval and adult heart regeneration ( Figure 4) (110, 111). 614 Methodologically, CL was administrated by IP, so we also cannot rule out the potential systemic 615 effects of losing CL sensitive cells elsewhere in the body on cardiac repair, despite we examined the 616 macrophage recruitment to the injured hearts after a long recovery period. In addition, resident 617 macrophages are usually confirmed through lineage tracing experiments to understand their 618 embryonic origin. Future studies using more sophisticated methods are required to confirm the 619 embryonic origin of resident macrophages and functionally validate the specific mechanisms of how 620 these resident macrophages may directly or indirectly contribute to CM proliferation and zebrafish 621 heart regeneration.

623
In summary, our study characterized the plasticity and heterogeneity of both macrophages and 624 neutrophils at steady-state and during the first-week post injury. In addition, we identified cardiac 625 resident macrophages and functionally characterized their roles in heart regeneration. New genetic 626 tools and models are required to further investigate the origin and function of distinct macrophage 627 subsets and the molecular mechanism underlying how cardiac macrophages contribute to heart 628 regeneration, which may gain knowledge and provide new insights for developing therapeutic 629 strategies for cardiac repair.  Cryoinjury was performed as previously described in zebrafish (115-117). In brief, fish were 647 anesthetized in 0.04% tricaine (Sigma, St Louis, MO, USA) and immobilized in a damp sponge with 648 ventral side up. A small incision was created through the thoracic wall by microdissection scissors. A 649 stainless steel cryoprobe precooled in liquid nitrogen was placed on the ventricular surface until 650 thawing. Fish were then moved back to a tank of freshwater for recovery, and their reanimation could 651 be enhanced by pipetting water onto the gills after surgery.

653
Cryosections and histological analyses 654 Zebrafish hearts were extracted and fixed in 4% (wt/vol) paraformaldehyde (Alfa Aesar, MA, 655 USA) at room temperature for 1 hour (hr). Collected hearts were subsequently cryopreserved with 656 30% (wt/vol) sucrose, followed by immersed in OCT (Tissue-Tek, Sakura Finetek, Torrance, CA) 657 and stored at −80°C immediately. 10 μm cryosections were collected for histological analysis. measuring the scar area in each heart from 5 discontinuous sections including the one with the largest 666 scar as well as the two sections at the front and the two sections at the back as previously described 667 (1). Statistic was performed on Prism 9 using student's t-test.

669
Quantitative polymerase chain reaction (qPCR) 670 For collecting single cells from hearts, 45-60 ventricles were isolated in each experiment at First-strand cDNA was synthesized using SuperScript™ III First-Strand Synthesis System (Life 679 Technologies Invitrogen, CA, USA) with oligo (dT) 18 primer. The qPCR analysis was carried out 680 using DyNAmo ColorFlash SYBR Green qPCR Kit (Thermo Scientific™, USA) on a LightCycler-® 681 480 Instrument II (Roche). The relative gene expression was normalized using eef1a1l1 as an 682 internal control. Oligonucleotide sequences for qPCR analysis were listed in Figure 7-source data 2. 683 684 Immunostaining and imaging 685 For immunofluorescence, slides were washed twice with PBS and twice with ddH 2 O, and then 686 incubated in blocking solution [1×PBS, 2% (vol/vol) sheep serum, 0.2% Triton X-100, 1% DMSO]. 687 Slides were then incubated in primary antibodies overnight at 4°C, followed by three PBST (0.1% 688 Triton X-100 in 1×PBS) washes and incubation with secondary antibodies for 1.5 hr at 28°C. Slides 689 were washed again with PBST, and stained with DAPI (Santa Cruz Biotechnology, Texas, USA) 690 before mounting. Antibodies and reagents used in this study included anti-Mpx (Genetex, San 691 Antonio, TX) at 1:500, anti-mCherry (Abcam, Cambridge, UK) at 1:250 and anti-HMOX1 (Aviva, 692 California, USA) at 1:100. EdU staining was performed using the Click-iT™ EdU Cell Proliferation 693 Kit for Imaging (Thermo Fisher, Massachusetts, USA) following the manufacturer's instructions. 694 EdU (25 μg/fish) was IP injected 24 hr before extracting the heart. TUNEL assay was performed 695 using In Situ Cell Death Detection Kit, TMR red (Sigma, St Louis, MO, USA) following the 696 manufacturer's instructions. 697 Imaging of whole-mount hearts and heart sections was performed using Nikon SMZ25 and 698 Zeiss LSM 800, respectively. Proliferating CM and the density of CM nuclei in the 200 μm border 699 zone directly adjacent to the injured area were quantified as previously described (57). 700 Revascularization was examined by live imaging of the endogenous fluorescence from vessel 701 reporter fish. Re-vascularized vessel density in the whole injured area was measured and quantified 702 by ImageJ. The student's t-test was applied to assess all comparisons by Prism 9. 703 704

In situ hybridization 705
In situ hybridization was performed on cryosections according to standard procedures. Briefly, 706 the templates of antisense digoxigenin (DIG)-labeled riboprobes for hbaa1, timp4.3 and hmox1a 707 were generated by PCR with the additional T7 promoter sequence on the primers. Oligonucleotide 708 sequences for probe synthesis are listed in Key Resources Quantification of the red color developed in the injured area was performed by ImageJ and the 728 student's t-test was applied for the statistics. 729 730 Next-generation RNA sequencing analysis 731 For each time point and condition, zebrafish ventricles from 3 fish hearts were used as 732 biological duplicates for the RNA-seq experiments. RNA extraction was performed as previously 733 described with minor modifications (1). Briefly, RNA isolation was done using the miRNeasy micro 734 Kit (Qiagen). RNA quality analysis was done using Qubit and Bioanalyzer at the NGS High 735 Throughput genomics core, Academia Sinica. Sequencing was performed on the HiSeq Rapid using Ingenuity Pathway Analysis (IPA) software (Qiagen) (122) following manufacturer protocols. 752 Log2 ratio was input from the normalized read counts in zebrafish and defined as DEGs at log2FC 753 above ±1.

755
Single-cell RNA-sequencing (scRNAseq) and bioinformatic analysis 756 Heart dissociation was followed the same procedures in qPCR. For each scRNA-seq sample, 757 35-45 cryoinjured ventricles were collected from each experiment at respective time-points and 758 conditions, pooled together and digested using LiberaseDH (Sigma, St Louis, MO, USA) at 28℃ in 759 a 24-well plate. Cell suspension was filtered through 100 μm, 70 μm and 40 μm cell strainers (SPL, 760 Gyeonggi-do, Korea) and centrifuged at 200g for 5 min. Cell pellet was resuspended in 0.04% 761 BSA/PBS, stained with DAPI and filtered through 35 μm Flowmi cell strainer (BD, New Jersey, 762 USA). Cells were sorted by FACS (BD Facs Aria) and subjected to single-cell sequencing 763 (scRNAseq) following cell counting with countess II automated cell counter (Invitrogen). 764 scRNAseq library was generated using single cell 3' Reagent kits in chromium platform (10x 765 Genomics). Cell ranger software suite (10x Genomics) was utilized for processing and de-766 multiplexing raw sequencing data (123). The raw base files were converted to the fastq format and 767 the subsequent sequences were mapped to zebrafish Ensembl genome assembly for processing. 768 Downstream analysis of the gene count matrix generated by CellRanger (10x Genomics) was 769 performed in R version of Seurat 3.1 (39, 124). The raw gene count matrix was loaded into Seurat 770 and a Seurat object was generated by filtering cells that expressed more than 400 nUMI counts and heatmaps (top enriched genes) was performed using functions within the Seurat package. Gene 792 enrichment analysis were determined using the above described RNA-seq pipeline. In brief, the 793 DEGs across cell types were used in WebGestalt to generate the cluster-specific biological processes 794 of GO and KEGG pathways (127). 795 In addition to basic scRNA-seq analysis steps, cell-cycle scoring and regression were also 796 performed. Zebrafish's cell-cycle related orthologs to human were identified based on Ensembl 797 annotation. All cells were assigned a cell-cycle score by "CellCycleScoring" function. Effects of cell-798 cycle heterogeneity were regressed out by "ScaleData" function with cell-cycle scores. 799 Dimensionality reduction of PCA and UMAP, and clustering were performed with the same 800 functions and parameters in the basic analysis steps (Figure 2-source data 1).

802
Pseudobulk analysis 803 We generated reference matrixes of the unsupervised and un-normalized mean counts computed 804 for each gene across all individual cells from each cell type using the Seurat and the 805 SingleCellExperiment package (128). The raw matrices were split by cell type and the associated 806 conditions of zebrafish (PBS-treated vs CL-treated). Normalization was performed on the sparse 807 gene matrix using total read counts and the total cell numbers for each cell type. The normalized 808 count matrix was applied for differential expression analysis using NOISeq (120). Furthermore, 809 average expression of genes across all the cell types detected by the scRNA-seq data from each time 810 point was calculated using Seurat function "AverageExpression", which was used as a reference 811 dataset for query. Ligand-receptor interaction network analysis 815 The zebrafish ligand-receptor pairs were derived from Human ligand-receptor pair database as 816 described previously (129). We mapped the ligand-receptor pair orthologous to the human database 817 and generated the fish database for use in interaction network. We kept the ligand-receptor pairs that 818 showed expression in upper-quantile normalization as the cutoff to define the expressed data. Next, 819 the macrophages and neutrophils that expressed the genes of ligands and receptors were determined 820 based on the counts from pseudobulk analysis with an upper-quantile expression at the given time 821 point. The total number of the potential ligand-receptor interaction events between each cluster of 822 macrophages and neutrophils were visualized using Cytoscape (130). 823

824
Pseudotemporal trajectory analysis 825 We adopted the macrophage and neutrophil populations from the datasets respectively after 826 importing the Seurat object using the as.cell_data_set and chooseCells function in Monocle3 (131, 827 132). UMAP dimensional reduction was applied to project the cells using the plotCells and 828 clusterCells function. To learn the pseudotemporal trajectory, we then used the learnGraph and 829 orderCells functions using modified parameters. We determined the DEGs over pseudotime using 830 the topMarkers function. The top 3 DEGs were selected by greater than 0.1 fraction of cells 831 expressing the given genes and plotted by plot_genes_by_group function. The following gene 832 expression of the cell types were also plotted as function of pseudotime to generate cell trajectory 833 analysis across pseudotime using the plot_genes_in_pseudotime function in Monocle3 ( Figure 6-834 source data 1  Clodronate liposomes (CL) one-day before cardiac cryoinjury. Injured hearts were collected at 7-and 1230 21-day post cryoinjury (dpci), respectively. Uninjured hearts were collected as the control of the 1231 infarcted hearts. Yellow highlights the cluster-enriched genes with gene names listed on the left. 1270 DEGs, differentially expressed genes. 1271 1272   was observed in zebrafish hearts (PBS or CL injection at 8 days before cardiac injury, -8d_PBS or -1423