Heterogeneous pdgfrb⁺ cells regulate coronary vessel development and revascularization during heart regeneration

Coronary heart disease is the leading cause of human mortality worldwide. Understanding the mechanisms of coronary vessel development and revascularization has thus drawn much attention in the hope of advancing the treatment of heart disease. In contrast to mammals, zebrafish is a well-established, genetically tractable model organism that shows remarkable regenerative capacity in the adult heart. Fast vascularization is essential for regeneration of damaged zebrafish hearts (Marin-Juez et al., 2016). The presence of well-structured coronary vasculature and ease of imaging the developing and adult hearts make zebrafish an ideal system for exploring cellular and molecular mechanisms of coronary vessel formation. We previously reported that the Cxcr4a-Cxcl12b chemokine axis guides newly emerged endothelial sprouts from the atrioventricular canal (AVC) to undergo angiogenesis and gradually cover the juvenile heart ventricle during development (Harrison et al., 2015). However, the cellular and molecular mechanisms that govern the maturation and maintenance of the developed coronary network remain unclear.

Mural cells (including both pericytes and smooth muscle cells) are a collection of diverse supporting cells that are recruited onto endothelial cells and cover the circulatory vessels as single or multiple cell layers. They regulate vessel development, stability, and physical functions such as vessel contractility (Armulik et al., 2011). Endothelial cells express Platelet-derived Growth Factor b (Pdgf-b) and its receptor β (Pdgfrβ) is expressed by mural cells; this signaling regulates mural cell recruitment onto endothelial cells (Armulik et al., 2005, 2011; Hellstrom et al., 1999; Lindahl et al., 1997; Lindblom et al., 2003; Winkler et al., 2010). Genetic disruption of Pdgfb or Pdgfrb significantly decreases mural cell coverage on vessels throughout the mouse embryo (Hellstrom et al., 1999). In the central nervous system (CNS), loss of mural cell coverage in the vasculature makes the blood vessels hyperplastic (significantly more endothelial cells per vessel), dilated, and susceptible to hemorrhage (Lindahl et al., 1997).

In mouse hearts, mural cells originate from epicardial (Cai et al., 2008; Mellgren et al., 2008) and endocardial (Chen et al., 2016) cells. During development, pericytes are recruited to microvessels and mature further to become smooth muscle cells on coronary arteries. Smooth muscle cell differentiation is induced with the commencement of blood flow through the arterial vessel, which in turn stimulates Notch activation in the pericytes (Volz et al., 2015). PDGFRβ signaling is essential for coronary smooth muscle development in both mouse and avian hearts (Mellgren et al., 2008; Smith et al., 2011; Van Den Akker et al., 2005, 2008) and pericytes are decreased in Pdgfrb null mice (Volz et al., 2015). It is not yet clear how diverse the cardiac mural cell populations are and how PDGFRβ signaling regulates different mural cell populations.

Here, we describe pdgfrb expression patterns and mutant phenotypes, and heterogeneity of pdgfrb⁺ cells during zebrafish coronary vessel development and heart regeneration. We find a subpopulation of pdgfrb⁺ mural cells associated with coronary arteries that also expresses cxcl12b. These pdgfrb; cxcl12b double-positive coronary arterial mural cells express markers of smooth muscle cells. Interestingly, these arterial mural cells remain associated with the coronary arteries whereas other (non-arterial) coronary vessels lose pdgfrb⁺ mural cells in the pdgfrb mutant. Single-cell RNA sequencing (scRNAseq) analyses reveal decreased smooth muscle gene expression, suggesting differentiation defects in these coronary arterial mural cells of pdgfrb mutants. Using a novel fluidic device-based culture and live-imaging system, we have further demonstrated that pdgfrb expression changes dynamically in the epicardium-derived cells (EDPCs) and that pre-existing mural cells migrate with angiogenic endothelial cells during heart regeneration. scRNAseq further reveal distinct populations of pdgfrb⁺ EDPCs and pre-existing mural cells, and differential gene expression analysis suggests that these cells might play important roles during heart regeneration. Furthermore, pdgfrb mutants show defects in mural cell association with coronary endothelial cells during heart regeneration; this severely compromises the regenerative response. We have identified mdka as one of the top differentially expressed genes in pdgfrb⁺ cells during heart regeneration. However, mdka mutants do not show severe defects in heart regeneration after amputation. Our results suggest that heterogeneous pdgfrb⁺ cells are essential for coronary development and heart regeneration.

To characterize spatiotemporal expression patterns of pdgfrb in the developing zebrafish heart, we utilized a pdgfrb:mCitrine transgenic reporter (Vanhollebeke et al., 2015). pdgfrb expression was consistently observed in the bulbus arteriosus (BA) and around the AVC at the late larval stage [24 days post-fertilization (dpf), ∼7 mm in standard body length (SL)], before any coronary endothelial cells emerged (Fig. 1A). In early juvenile fish, pdgfrb expression spread along the AVC at 31 (∼10.5 mm SL) and 43 dpf (∼17 mm SL). The coronary endothelial cells [marked by Tg(fli1a: DsRed) in Fig. 1A″] emerged on the ventricle around 36-49 dpf (∼14-18 mm SL). At the late juvenile stage by 55-74 dpf, pdgfrb expression became localized in the mural cells accompanying the nascent coronary endothelial cells sprouting out from the AVC (Fig. 1A‴, Fig. S1A).

As growing coronary vessels covered the ventricle, mural cells remained associated with large and small coronary vessels (Fig. 1A‴). Throughout development, mural cell density remained greater on the coronary vessels closer to the AVC (defined as the area proximal to the AVC) than the growing end of the vessels (regions distal to the AVC) (Fig. 1B, Fig. S1A). The colocalization of pdgfrb:mCitrine and tcf21:DsRed (which marks epicardium) around the AVC at 33 dpf (Fig. S1B) suggests that the origins of pdgfrb⁺ cells might be the epicardium. To confirm this, we performed lineage-tracing experiments using the tcf21:CreERT2 line combined with pdgfrb:Citrine and the ubi:Switch reporter. We observed that 76.15% of pdgfrb⁺ cells were co-­labeled with mCherry, indicating that most pdgfrb⁺ cells were derived from the tcf21⁺ epicardial lineage (Fig. 1Ca-Cb′, white arrowheads in Cc-Cc″, Fig. S1C). 23.85% of pdgfrb⁺ cells (Fig. S1C) were not co-labeled with mCherry, suggesting that either CreER-loxP recombination was not efficient in these cells or these pdgfrb⁺ cells were derived from other sources (Fig. 1Cd). Furthermore, we observed tcf21:CreERT2-labeled cells that did not express pdgfrb (Fig. 1Cc-Cc″, blue arrowheads), consistent with the notion that epicardium can also contribute to other cell types, such as fibroblasts in zebrafish (Sanchez-Iranzo et al., 2018) and mice (Ivey et al., 2018).

Next, we investigated whether the appearance of pdgfrb-expressing mural cells depends on blood circulation through the developing coronary vessels. We performed coronary angiography using Tg(pdgfrb:Citrine; fli1a:DsRed) juvenile zebrafish. pdgfrb:Citrine⁺ mural cells were observed at 43 dpf on nascent coronary vessel sprouts in the absence of blood flow (Fig. 1D). By 50 dpf, initiation of blood flow through major vessels was observed in some hearts. However, in these juvenile fish pdgfrb:Citrine⁺ mural cells were already attached to the immature vessel plexus without any detectable blood circulation (Fig. 1D′). All major vessels had blood flow by 64 dpf (Fig. 1D″). These results indicated that blood flow is not necessary for pdgfrb:Citrine⁺ mural cell association with coronary vessels.

To determine the role of Pdgfrβ during zebrafish coronary vessel development, we examined pdgfrb^um148 homozygous mutants (pdgfrb^−/−; Kok et al., 2015). Consistent with findings in Pdgfrb knockout mice (Lindahl et al., 1997), the brain vasculature of adult pdgfrb^−/− fish became significantly dilated with reduced branching density and decreased mural cell association (Fig. S2A). We did not observe hemorrhage in the heart ventricle. In pdgfrb^−/− and pdgfrb heterozygous (pdgfrb^+/−) mutants, coronary vessels did develop, although not as efficiently as in wild-type (WT; pdgfrb^+/+) fish. Both pdgfrb^−/− and pdgfrb^+/− ventricles had reduced coverage of coronary vessels at 98 dpf (Fig. 2A, Fig. S2B′). Unlike in pdgfrb^+/− fish, in which coronary vessel development was delayed, these defects in coronary vessel coverage remained by 167 dpf in pdgfrb^−/− homozygous mutants (hereafter referred to as pdgfrb mutants) (Fig. 2A,A′). This phenotype is consistent with what was reported recently using a different pdgfrb^sa16389 allele (Ando et al., 2021). Isolated endothelial cells, which fail to form continuous vessel networks, were observed in pdgfrb mutants (Fig. 2A, Fig. S2B, asterisks). Mural cell association with small vessels decreased significantly, but association with some of the large vessels remained unaffected in the pdgfrb mutants (Fig. 2B,B′, Fig. S2C). Further analysis revealed that among all large vessels, narrow large vessels (i.e. coronary arteries; Harrison et al., 2015) maintained pdgfrb⁺ mural cell association (Fig. 2B, Fig. S2C, white arrowhead), whereas the wide large vessels (vein-like) lost a significant number of pdgfrb⁺ mural cells (Fig. 2B, Fig. S2C, yellow arrowhead). Interestingly, in contrast to brain vessels that became dilated (Fig. S2A), the diameter of both wide large vessels and coronary arteries significantly decreased in pdgfrb mutants (Fig. S2D), suggesting a potential spatiotemporally specific mechanism of blood vessel formation/maturation regulation by pdgfrb⁺ mural cells. Taken together, these data suggest that different mechanisms might be utilized to regulate mural cell recruitment or maintenance along different subtypes of coronary vessels.

We previously reported that epicardium-derived cxcl12b:Citrine-expressing mural cells surround cxcr4a⁺ arterial endothelial cells (Harrison et al., 2015). Here, we also observed pdgfrb⁺ cells covering coronary arteries (Fig. 2B, Fig. S2C,C′, narrow large vessels, white arrowheads). Therefore, we examined whether cxcl12b:Citrine⁺ mural cells also express pdgfrb:EGFP. We observed that most (82.39%) of the mural cells on the coronary arteries express both cxcl12b and pdgfrb reporters (Fig. 3A′-A‴, blue arrowheads, 3A″″). pdgfrb:EGFP⁺ only (Fig. 3A′-A‴, green arrowheads) and cxcl12b:mCitrine⁺ only (Fig. 3A′-A‴, red arrowheads) mural cells were less abundant on coronary arteries (12.88% and 4.73%, respectively) (Fig. 3A″″).

We examined whether these cxcl12b⁺ cells along the coronary artery are affected in pdgfrb mutants. Consistent with the finding using the pdgfrb:EGFP reporter (Fig. 2B, Fig. S2C), we found that cxcl12b:Citrine⁺ (also pdgfrb:EGFP⁺) mural cells remained associated with coronary arteries whereas other coronary vessels lost mural cells. Furthermore, the diameter of the coronary artery was significantly decreased in pdgfrb mutants (Fig. 3B). These results suggest that pdgfrb⁺ mural cells are heterogeneous on the heart ventricle. The non-arterial mural cells express only pdgfrb and most of the mural cells on the coronary arteries are double-positive with cxcl12b. They also respond differently to the loss of pdgfrb.

To investigate further the gene expression signature in these pdgfrb:EGFP⁺ and cxcl12b:Citrine⁺ double-positive (pdgfrb⁺/cxcl12b⁺) mural cells and other pdgfrb:EGFP⁺ only cells in WT controls versus pdgfrb mutants, we performed scRNAseq. We sorted pdgfrb:EGFP⁺ cells from pdgfrb mutants and sibling WT control fish by fluorescence-activated cell sorting (FACS) to enrich the mural cells, although different cell types (especially cardiomyocytes) were also collected, likely owing to autofluorescence (Fig. S3A,B,D, Table S1). Uniform manifold approximation and projection (UMAP) analysis revealed that the pdgfrb⁺ and egfp⁺ cells are mainly in clusters 3 and 6, whereas cxcl12b⁺ cells are specifically in cluster 6. Furthermore, more cells expressed egfp transcripts at a higher level in cluster 6 (Fig. 3C, Fig. S3C). Differential gene expression analyses showed that the epicardial markers tcf21 and tbx18 are specifically found in cluster 3, suggesting that this cluster consists of epicardium and EPDCs. Cells in cluster 6 expressed more mural cell and smooth muscle cell markers (acta2, myh11a, tagln) (Fig. 3C′, Fig. S3C, Table S2). Furthermore, more cells in cluster 6 express higher notch3, which regulates vascular pericyte differentiation into smooth muscle cells (Volz et al., 2015) (Fig. 3C′, Fig. S3C′). Among other established mural cell markers (He et al., 2016; Whitesell et al., 2019), rgs5a and cd248a as well as a new marker, ndufa4I2a, were expressed in cluster 6 (Fig. 3C′) whereas kcne4 was expressed in both clusters (Fig. 3C, Fig. S3C). Other brain pericyte markers, abcc9, desma, cspg4, anpepa (He et al., 2016), showed minimal expression in very few cells, indicating overall heterogeneity of the coronary mural cell population and their differences from brain mural cells. Furthermore, these data suggest that cluster 6 is likely composed of coronary mural cells. We validated expression of some of these mural cell markers by performing RT-PCR. The smooth muscle cell marker acta2 showed increased expression in the pdgfrb:EGFP high cells, which likely represent the cxcl12b⁺/pdgfrb⁺ cells (Fig. 3C‴, Fig. S4). Because these cluster 6 mural cells express notch 3, we determined Notch activity using the Tp1bglb:Venus-PEST reporter (Ninov et al., 2012) and found that the coronary arterial mural cells show high Notch activity (Fig. 3C″, white arrowheads). To characterize the arterial mural cells further, we selected cxcl12b⁺ cells from cluster 6 and compared them with cxcl12b⁻ (pdgfrb⁺ only) cells. Differential gene expression analysis (Table S3) showed significantly higher smooth muscle marker gene (acta2, tagln) expression in the cxcl12b⁺ cells, indicating that the coronary arterial mural cells are smooth muscle-like cells (Fig. 3D).

In the pdgfrb mutant (pdgfrb^−/−), cluster 6 was most affected compared with other cell clusters (Fig. S5A). Cluster 6 contains ∼11% of all WT cells and only ∼2% of all pdgfrb^−/− cells (Fig. S5A′). The WT cells and pdgfrb mutant cells took distinct positions in the UMAP of the cluster 3 and 6, reflecting their overall gene expression differences (Fig. S5B). The comparative scatterplot of average gene expression across WT and pdgfrb^−/− conditions revealed that cluster 6 cells have more differential gene expression than cluster 3 cells (Fig. S5C). The mural cell markers were differentially expressed, with kcne4 the most significantly different, in the WT cluster 6 cells compared with the pdgfrb^−/− mutant cells (Fig. S5C, Table S4). These results indicate impaired smooth muscle cell differentiation in pdgfrb mutants even though low level expression of the smooth muscle and mural cell markers suggests that these cells still maintain mural cell identity. In contrast, the gene expression changes in the epicardial cluster (cluster 3) were less significant (Figs S3C and S5C). Among the epicardial markers, tbx18 and tcf21 showed prominent expression in the cluster 3 (epicardial cells) in both pdgfrb mutants and WT controls and much less expression in the cluster 6 (mural cells) (Fig. 3C′, Fig. S3C). These results suggested that in pdgfrb mutant heart, remaining epicardial cells and EPDCs (cluster 3) are less affected compared with differentiated mural cells (cluster 6), which decrease in number and change gene expression (Fig. S5A,A′,C).

Subclustering of the mural cells (cluster 6) revealed that subcluster 1 and 2, which have relatively higher pdgfrb expression and mural cell marker expression, are absent in pdgfrb mutants. Subcluster 0 cells have low pdgfrb but comparable cxcl12b expression with subcluster 1 and 2. Smooth muscle cell marker expression (tagln, acta2, myh11a) in subcluster 0 was also decreased in pdgfrb mutants (Fig. S6A,B, Table S5). A few cells in the epicardial cluster 3 also expressed some mural cell marker genes at a low level (kcne4, tagln, notch3) (Fig. 3C′, Fig. S3C). Further subclustering showed that subclusters 0, 1 and 3 express these mural cell markers along with epicardial markers and the fibroblast marker pdgfra. Subcluster 2 only expressed epicardial markers and likely are epicardial cells. Subcluster 4 expressed the mural cell marker cd248a and kcne4, along with the fibroblast marker postnb and pdgfra (Rajan et al., 2020) (Fig. S6C,D, Table S6). Thus, the epicardial pdgfrb-expressing cells are heterogeneous. Subclusters 0, 1, 3 and 4 are likely EPDCs showing genetic signatures of both mural cells and fibroblasts.

For functional analysis of pdgfrb-expressing cells, we performed gene ontology (GO) term analysis for differentially expressed genes in cluster 3 (epicardial cells) and cluster 6 (mural cells) cells compared with each other. The mural cells showed enrichment for GO terms related to development and morphogenesis, including ‘neuronal development’ and ‘blood vessel development/angiogenesis’. The epicardial cells showed GO term enrichment for ‘protein/peptide biosynthesis’ and ‘collagen and extracellular matrix synthesis’ (Fig. S7, Tables S7 and S8).

Our previous observation of pdgfrb upregulation and the presence of mural cells in the regenerating area of the zebrafish heart (Kim et al., 2010) prompted us to characterize pdgfrb expression patterns during regeneration. The apical regions of Tg(pdgfrb:Citrine;fli1a:DsRed) fish were injured and imaged at different days post-amputation (dpa). At 1 dpa, there were no obvious pdgfrb-expressing regions near the wound area except in mural cells around the pre-existing coronary vessels. Starting from 3 dpa, patchy pdgfrb expression was observed around the border of the amputated area (Fig. 4A). This expression gradually expanded significantly at 7 dpa and 10 dpa, plateaued at 14 dpa, and decreased at 30 dpa (Fig. 4A,A′, Fig. S8A). The mural cells associated with pre-existing coronary endothelial cells continued to express pdgfrb (Fig. S8A). This patchy pdgfrb expression was likely in epicardium or EPDCs because they were also positive for the epicardial marker tcf21 (Fig. 4B). These pdgfrb⁺ epicardial cells or EDPCs migrated and enveloped the regenerating area whereas the pdgfrb-expressing mural cells remained associated with the coronary endothelial cells migrating into the regenerating area.

We developed a novel fluidic device-based long-term explant culture for live imaging (Yip et al., 2020) to confirm the dynamic changes in pdgfrb expression patterns and potential interactions between pdgfrb⁺ cells with other cell types. Hearts from transgenic zebrafish Tg(pdgfrb:EGFP; fli1a:DsRed) were injured and allowed to recover in vivo for 5-10 days, dissected out from the fish, then placed in the device for live imaging for 72-120 h. pdgfrb:EGFP expression at 5-6 dpa showed diffuse expression in the epicardium within most of the wound site and more localized expression in mural cells towards the edge and outside the wound site (Movie 1). As regeneration proceeded, pdgfrb⁺ epicardial cells closed in to cover the entire wound area by 6 dpa. Following this, punctate, mural cell-like expression was observed within the regenerating region. The epicardial cells, or the diffuse expression observed within them, were highly dynamic within the wound area compared with the more stable behavior of existing mural cells (Movie 2). These mural cells migrated into the wound site alongside the individual endothelial cell to which they were attached (Movies 2 and 3). However, it was the dynamic epicardial expression or migration that preceded the endothelial cell migration into the wound site. The pre-existing mural cell expression moved more slowly alongside the cell body of the endothelial cells, behind the filopodial extension (Movies 2 and 3). The imaging data described above suggest that two different pdgfrb⁺ populations are present at the wound site during zebrafish heart regeneration: the pre-existing mural cells and a second epicardial cell-derived population. Furthermore, a subcomponent of the epicardium may be dynamically changing gene expression and transforming cell identity.

We hypothesized that pdgfrb⁺ populations with essential roles in zebrafish heart regeneration might have more active gene expression. To characterize further the diverse pdgfrb⁺ cell populations and determine their gene expression signatures, we performed FACS using WT adult Tg(pdgfrb:EGFP; cxcl12b:Citrine) zebrafish and FACS-sorted EGFP-positive cells from ventricles of uninjured and injured hearts at 7 dpa to perform scRNAseq. The sorted EGFP⁺ cells form 15 clusters based on their gene expression differences (Fig. S9A) and differentially expressed genes were identified (Fig. S9B, Table S9). The epicardial/EPDCs/mural cell cluster was identified based on marker gene expression (Table S9, Fig. S9C,D). Three-hundred and twelve cells from uninjured hearts and 199 cells from 7 dpa hearts (clusters 5, 6, 10 and 12 combined) were analyzed further (Fig. 5A,A′). We found that pdgfrb expression is mainly detected in clusters of epicardial cells/EPDCs/mural cells (clusters 5, 6, 10 and 12) (Fig. S9C,D). The heterogeneity of tcf21:nucEGFP⁺ epicardial cells from uninjured hearts was previously reported (Cao et al., 2016). Therefore, we focused our gene signature analysis on the regenerating hearts.

By epicardial marker (tcf21, tbx18) and mural cell marker (rgs5a, cd248a, kcne4) (Cho et al., 2003; Venero Galanternik et al., 2017) gene expression, clusters 5, 6 and 10, were found to be epicardial/EPDC clusters, which have very few cells with mural cell marker expression, and cluster 12 was found to be the main mural cell cluster (Fig. S9D). Interestingly, when pdgfrb expression was checked across these clusters, the mural cell cluster (cluster 12) showed high pdgfrb expression but there was little to no difference in expression level between mural cells from uninjured and injured hearts. In all epicardial clusters, more cells from the injured hearts showed pdgfrb expression than uninjured heart cells (Fig. S9C). This is consistent with our observation that dynamic expression changes of pdgfrb in the injured hearts occur in the epicardial cells (Fig. 4). To characterize further overall gene expression changes in response to heart injury, differential gene expression analysis was performed combining all pdgfrb-expressing cluster cells (combining clusters 5, 6, 10 and 12 cells) for 7 dpa hearts compared with uninjured hearts (Fig. 5A′-C, Table S10).

Functional characterization of these pdgfrb-expressing cells was carried out by GO term analysis based on the differentially expressed genes in the injured hearts. After categorization of all enriched GO terms, it was found that functionalities related to extracellular matrix (ECM), collagen (16% of all GO terms), regeneration and development (12%), endopeptidase inhibitor activities (8%), supramolecular polymer changes in cytoskeleton (8%), extracellular space (6%), and response to oxygen level/ hypoxia (6%) are enriched in the differentially expressed genes (Fig. 5D, Table S11). These functions are well-aligned with the requirements of a regenerating heart where ECM deposition occurs immediately after amputation. As the healing process progresses, cell migration and wound response/regeneration occur.

Fig. 5B shows the top 15 differentially expressed genes among the genes upregulated in the 7 dpa injured hearts based on the most significant adjusted P-values; some of them (anxa2a, mdka, pdgfrl, si:ch211-198c19.3) were validated by RT-PCR using FACS-sorted pdgfrb:EGFP⁺ cells from 7 dpa hearts compared with uninjured hearts (Fig. S10B). The top 5 differentially expressed genes were: zgc:152791, annexin A2a (anxa2a), fibronectin 1b (fn1b), midkine a (mdka) and periostin b (postnb) (Fig. 5B). anxa2a has been recently shown to be required in zebrafish caudal fin regeneration (Quoseena et al., 2020). The ECM proteins Fibronectin 1b and Periostin b were previously shown to be upregulated/involved zebrafish heart regeneration (Rodius et al., 2014; Sanchez-Iranzo et al., 2018; Wang et al., 2013). mdka, encoding a growth factor, was previously shown to be induced after zebrafish heart injury, but its functions in heart regeneration after amputation have not been characterized (Lien et al., 2006).

Other upregulated genes at 7 dpa encode proteins including Hyaluronan and Proteoglycan link protein 1a (hapln1a), which is predicted to have hyaluronic acid-binding activity and has been shown to be involved in fin regeneration (Ouyang et al., 2017), epicardial EMT and heart regeneration (Missinato et al., 2015). Collagen type 1, alpha 1a (col1a1a) has also been shown to be involved in fin development and regeneration (Duran et al., 2015, 2011; Padhi et al., 2004). Thymosin beta 1 (tmsb1) and Chemokine ligand 8a (cxcl8a) are differentially expressed in clusters 5, 10 and 12, respectively (Table S9). Translation initiation factor (eif4ebp3l) and lipoprotein lipase (lpl) are two highly expressed genes in the uninjured hearts (Fig. S11A). Consistent with our previous findings (Kim et al., 2010), expression of epithelial-to-mesenchymal transition (EMT) genes (snai1a, snai1b, snai2, twist1b) were induced during heart regeneration (Fig. S11B). Overall, these data suggested that in response to zebrafish heart amputation the heterogeneous pdgfrb-expressing EPDCs and mural cells (clusters 5, 6, 10, 12) actively express genes with potential roles in heart regeneration (Fig. 5C). Further functional characterization by GO term enrichment analysis for individual pdgfrb⁺ clusters indicated that the mural cell cluster (cluster 12) and epicardial/EPDC/fibroblast cluster (cluster 6) have the most potential roles in supporting tissue development and morphogenesis (17% of all GO term enrichment for each cluster). The mural cell cluster had the highest GO term enrichment for the circulatory system development (10% of all GO terms for cluster 12) followed by cluster 6 cells (7% of all GO terms for cluster 6). For clusters 5 and 10, 5% of each cluster's GO terms were related to circulatory system development. For cluster 10, which is the cluster expressing postnb in uninjured heart, 8% of the GO terms were related to immune system regulation (Fig. S12, Tables S12-S15).

We previously reported that fish treated with a Pdgfrβ inhibitor showed defects in revascularization during heart regeneration (Kim et al., 2010). Here, we utilized adult pdgfrb mutants to examine further the requirement of pdgfrb during zebrafish heart regeneration. We observed very few blood vessels in the regenerating area of pdgfrb mutant hearts at 21 and 34 dpa. These vessels were large, had significantly fewer branches, and were mostly devoid of pdgfrb⁺ mural cells. In contrast, sibling control fish had dense networks of pdgfrb⁺ mural cell-covered coronary vessels in the regenerating region (Fig. 6A). These data suggest that Pdgfrβ signaling plays an essential role during revascularization. The pdgfrb mutants failed to regenerate after ventricular resection and maintained a fibrotic scar (Fig. 6B), confirming that Pdgfrβ signaling is essential for heart regeneration. mdka was found as one of the top differentially expressed genes (Fig. 5B) and showed increased expression in 7 dpa hearts. We confirmed the upregulation of mdka in EDPCs by qRT-PCR using FACS-sorted pdgfrb:EGFP⁺ cells and by in situ hybridization. qRT-PCR using 7 dpa whole heart ventricles further validated the upregulation of mdka whereas mdkb expression did not increase (Fig. S13A-C). Midkine (Mdk) has been shown to play important roles in tissue regeneration in many different organs in multiple species (Ang et al., 2020; Ikutomo et al., 2014; Nagashima et al., 2020; Tsai et al., 2020). However, in contrast to pdgfrb mutant hearts, mdka mutants did not show significant defects in heart regeneration (Fig. S13D) and collagen deposition (Fig. S13E,E′), despite its strong expression.

The cell compositions, functions and origins of the mural cells in the hearts remain incompletely understood. Focusing on pdgfrb⁺ cells, we determined the heterogeneity of these cardiac mural cells and EPDCs in zebrafish heart during development and regeneration and investigated how Pdgfrβ signaling affects these different cell populations. We found that, during heart development, pdgfrb and cxcl12b double-positive (pdgfrb⁺;cxcl12b⁺) cells line coronary arteries whereas mural cells expressing only pdgfrb surround non-arterial vessels. Furthermore, pdgfrb⁺ EPDCs are less affected than mural cells in pdgfrb⁺ mutant hearts. In adult regenerating hearts, pdgfrb⁺ cells were identified in pre-existing mural cells and EPDCs. The use of scRNAseq and confocal and live imaging allowed us to delineate unique gene expression signatures in different populations, revealing their different functions.

Although very few pdgfrb⁺;cxcl12b⁺ double-positive cells were captured in our scRNAseq analysis, smooth muscle gene expression in these double-positive cells suggests that they are more differentiated than the pdgfrb⁺ only cells. Interestingly, the association of these pdgfrb⁺;cxcl12b⁺ mural cells with coronary arteries does not depend on Pdgfrβ signaling and all non-arterial vessels lose pdgfrb⁺ mural cells in pdgfrb mutants. Recent scRNAseq data of mouse brain revealed two subclasses of mural cells; pericytes are in a continuum with venous smooth muscle cells, which are distinct from arteriole or arterial smooth muscle cells (Vanlandewijck et al., 2018). It is possible that pdgfrb⁺;cxcl12b⁺ mural cells are similar to arteriole or arterial smooth muscle cells, which are different from the pdgfrb⁺ only cells surrounding the wide large vessels that are likely the coronary veins. In mice, smooth muscle cells along the coronary arteries are derived from Pdgfrb⁺ pericyte progenitors after the blood flow starts but the smooth muscle cells around the coronary arteries fail to form in Pdgfrβ knockout mice (Volz et al., 2015). However, these pdgfrb⁺;cxcl12b⁺ double-positive cells remained associated with coronary arteries in pdgfrb mutant fish, suggesting that a previously unknown mechanism exists to maintain the mural cell association with coronary arteries. One likely mechanism is that Cxcl12-Cxcr4 signaling (Stratman et al., 2020) might activate an alternative signaling pathway in the arterial mural cells that acts in parallel with PDGF signaling for mural cell recruitment. The scRNAseq identified several signaling pathways that could mediate this Pdgfrβ-independent mural cell association with coronary endothelial cells and this will be pursued in future investigation. Furthermore, our scRNAseq analyses also revealed new mural cell markers, kcne4 and ndufa4I2a, that are affected in pdgfrb mutants, consistent with recently reported results in zebrafish embryos (Shih et al., 2021).

Pdgfrβ signaling and pdgfrb⁺ cells have been implicated in the development and regeneration of different tissues and organs. Consistent with our finding, intra-myocardial delivery of PDGF-BB provides myocardial protection and improves ventricular functions (Hsieh et al., 2006a,b). Nonetheless, our current findings cannot distinguish between a phenotype in EPDCs versus a developmental phenotype in coronary vessels or revascularization that resulted in the impaired heart regeneration in pdgfrb mutants. A temporospatial-specific knockdown of pdgfrb during heart regeneration will be performed in the future to elucidate this mechanism further. Recently, it was demonstrated that ECM derived from pdgfrb⁺ myoseptal and perivascular cells prevents scarring and promotes axon regeneration of zebrafish spinal cord in a Pdgfrβ signaling-dependent manner (Tsata et al., 2021). Consistent with this finding, our scRNAseq analyses identified significant changes (16% of all GO term enrichment) in genes encoding ECM components. How Pdgfrβ signaling regulates the ECM in the regenerating heart and shapes the regenerative environment will be of interest for future studies.

Taking a candidate approach to examine further the differentially expressed genes identified from our scRNAseq, we decided to first assess the roles of mdka, which was also discovered by microarray gene expression profiling previously (Lien et al., 2006). mdka has been shown to play an important role in neural regeneration in the retina by regulating cell cycle progression of Müller glia (Nagashima et al., 2020), epimorphic regeneration (Ang et al., 2020) and regeneration after skeletal muscle injury (Ikutomo et al., 2014). Furthermore, it was recently reported that axolotl Midkine (Mdk) could regulate wound epidermis development and inflammation during the initiation of limb regeneration (Tsai et al., 2020). mdka is highly expressed in epicardium and EPDCs after amputation and this led us to examine its role during heart regeneration. We did not observe any significant defects in heart regeneration and increased fibrotic scarring remained at the injury site at 34 dpa compared with controls. However, it was reported that loss of mdka after cryoinjury decreased proliferation of endothelial cells and retention of a collagen scar (Grivas et al., 2021). We cannot exclude the possibility that different mdka mutant alleles and injury models might reveal regeneration defects. Examining the phenotypes of different mdka mutant alleles in cryoinjury versus amputation injury models is of interest for future investigation; however, is beyond the scope the current study. Nonetheless, our scRNAseq data provide new candidate genes that may have a role in coronary vessel development and heart regeneration.

The following zebrafish lines were raised and maintained at Children's Hospital Los Angeles (CHLA) under standard conditions of care and with CHLA IACUC oversight. IACUC approved all experimental procedures used in this study. Both males and females were utilized in this study. The age and length of the fish are specified in the results. The fish lines used in this study were: Tg(fli1a.ep:DsRedEX)^um13 (also known as fli1a:DsRed; Covassin et al., 2007), Tg(fli1a:EGFP) (Lawson and Weinstein, 2002), TgBAC(pdgfrb:EGFP) (Ando et al., 2016), TgBAC(pdgfrb:Citrine) (Vanhollebeke et al., 2015), pdgfrb^um148 (Kok et al., 2015), Mdka^mi5001(Nagashima et al., 2020), Tg(cxcl12b:Citrine) (Bussmann et al., 2011), TgBAC(tcf21:Cre-ERT2)^pd42 (Kikuchi et al., 2011), Tg(-3.5ubi:loxP-EGFP-loxP-mCherry) (Mosimann et al., 2011) and Tp1bglb:Venus-PEST (Ninov et al., 2012).

Lineage tracing tcf21:CreERT2 was performed as described by Harrison et al. (2015) with pdgfrb:Citrine as the mural cell marker. Tg(pdgfrb:Citrine; tcf21:CreERT2; ubi:loxp-EGFP-loxp-mCherry) fish were generated and the embryonic epicardial cells were labeled by mCherry by activating tcf21:CreERT2 by administering 4-hydroxytamoxiphen (4-OHT) during the first 5 days of embryonic development. 4-OHT (13 µg/ml) was added to E3 medium and the medium was changed daily. After treatment, the fish were raised into adulthood. The adult hearts carried all the progenies derived from the embryonic-labeled epicardium. pdgfrb:Citrine was used to identify the mural cells. The percentage of mCherry-labeled Citrine-positive cells was quantified.

Dextran Cascade Blue dye was injected through the retro-orbital sinus as described (Harrison et al., 2015). Briefly, injection was performed 2.5-3.0 h before heart extraction. Zebrafish were anesthetized in a tricaine solution and placed horizontally on a sponge so that the right eye was facing up. A Hamilton syringe loaded with 4 pl Dextran Cascade Blue was inserted at 30° to the plane of the sponge, about 2 mm deep at the 7 o'clock position of the right eye socket (Pugach et al., 2009). The amputation experiments were performed as described (Poss et al., 2002).

Hearts were immobilized in 1% low melting point agarose in PBS. The heart was oriented with the ventricle at the center of the image, the bulbus arteriosus at the top, and the atrium on the left. z-stack images were collected using a Zeiss 710 confocal microscope. z-stack images were then converted to a maximum intensity projection before further processing/quantification. ImageJ/Fiji software was used for quantification.

Hearts from transgenic zebrafish [Tg(pdgfrb:EGFP; fli1a:DsRed)] were amputated and allowed to recover in vivo for 5-10 days, after which hearts were removed into imaging media [L15 (300 ml), 10% fetal calf serum/fetal bovine serum, 100 μg/ml Primocin, 1.25 mM CaCl₂, 800 mg/l glucose, Pen/Step]. Hearts were cleaned of any external blood or attached tissue debris. The imaging microfluidic device was prepared as described (Yip et al., 2020) and four hearts were then placed into the imaging wells under Ringer's solution. The device was sealed and mounted onto a Zeiss cell observer system equipped with a Hamamatsu ORCA-flash4.0LT and Colibri 7 LED light source for live imaging. The imaging device was maintained at a temperature of 28.5°C with a PeCon atmospheric control stage and cover. Imaging media was perfused through the imaging wells at a rate of 0.5 ml/h and imaging was carried out for 72-120 h. Acquisitions were carried out to capture a z-stack of 35-50 images 6 μm apart every 30 min. Collected images were cropped (in z-plane and time) and deconvolved (using AutoQuant). ImageJ/Fiji was used to produce the final moves with the aid of a Gaussian-based focusing macro and a z sub-selection macro (available on request).

Mural cell association with the coronary vessels was quantified as the mural cell number per unit vessel length. Using ImageJ software, freehand lines were drawn along the length of randomly selected vessels of certain size (e.g. Fig. 2B′) or from certain regions of the heart ventricles (e.g. Fig. 1B). The length (by pixel numbers) of the freehand line was measured by ImageJ and pdgfrb⁺ cells along the line were counted manually. Then, the pdgfrb⁺ cell numbers were divided by the freehand line's length to get the measurement of the mural cell number per unit vessel length. Paired t-test was used to quantify the significance of the value differences between different samples.

Using ImageJ freehand selection tools, areas covered by coronary vessels following the tips of the vessels were selected in heart ventricles. Respective areas were measured and the areas following the coronary vessels were divided by whole ventricle area and corresponding percentages were quantified (Fig. 2A′). One-way ANOVA was used to quantify the significance of the value differences between pdgfrb^+/+, pdgfrb^+/− and pdgfrb^−/− hearts.

pdgfrb expression was quantified by measuring the area of the pdgfrb-expressing region as a percentage of the area of the whole apical view of the ventricle (Fig. 4A′). The area of the pdgfrb-expressing region (marked by the white small-dotted line in Fig. 4A) was quantified using the freehand selection tool of ImageJ software. One-way ANOVA was used to quantify the significance of the measurement difference between 3 dpa and 7 dpa.

The heart ventricle was immersed in the modified Tyrode's solution (Tessadori et al., 2012) on ice. After washing blood from the tissue, the ventricle was torn open and into small pieces by forceps and transferred into Ca²⁺-free modified Tyrode's solution (Tessadori et al., 2012) at 30°C. The tissue was then washed twice with ice-cold Ca²⁺-free modified Tyrode's solution. The tissue was then treated with the digestion mix (500 µl per five hearts) containing Liberase (500 CDU); Elastase (3.1 U) and DNaseI (32 U) at 33-35°C for ∼15 min with continuous stirring. The digestion was stopped by adding ice-cold Ca²⁺-free modified Tyrode's solution containing 10% fetal bovine serum and DNaseI (32 U). The solution was then filtered through a 40-μm sieve and the cells were pelleted by centrifuging at 500 g for 5 min at 4°C. The cells were washed with ice-cold Ca²⁺-free modified Tyrode's solution containing 30% fetal bovine serum. The cells were then precipitated again and redissolved in an appropriate volume of ice-cold Ca²⁺-free modified Tyrode's solution containing 30% fetal bovine serum to maintain a suitable cell density for single-cell cDNA library preparation or qRT-PCR.

pdgfrb:EGFP cells were FACS-sorted from Tg(pdgfrb:EGFP; cxcl12b:Citrine) 6-month-old WT and pdgfrb^−/− fish. pdgfrb:EGFP cells were FACS-sorted from approximately ten Tg(pdgfrb:EGFP; cxcl12b:Citrine) 18-month-old uninjured or injured (7 dpa) fish. Heart ventricles were dissociated into single-cell suspension following the abovementioned method for each sample. pdgfrb:EGFP cells were FACS-sorted using a BD FACSAria cell sorter. The live cells were sorted by filtering out dead cells with DAPI staining. EGFP only, DsRed only transgenic samples were used as a positive control, and WT fish without any transgenic reporter were used as the negative control. During sorting, the FACS gates were set to maximize inclusion of all GFP-positive cells (low to high EGFP signals) to isolate all possible EGFP-expressing cell populations. This wide gate set-up led to the inclusion of other cell types (e.g. cardiomyocytes) that do not express pdgfrb:EGFP, likely due to background fluorescence. After sorting, the cells were collected into Ca²⁺-free modified Tyrode's solution containing 30% fetal bovine serum and kept on ice. This solution was then centrifuged at 500 g for 5 min at 4°C to concentrate the sample for downstream procedures (e.g. single-cell RNA sequencing, qRT-PCR). For qRT-PCR, the FACS gates for sorting GFP⁺ cells were divided to isolate cells expressing GFP at a high level separately from those expressing GFP at a moderate level based on our confocal imaging and scRNAseq observations that coronary arterial mural cells are GFP^high cells whereas GFP^low cells include mostly epicardial cells.

For the generation of single-cell gel beads in emulsion, cells were loaded on a Chromium single cell instrument (10X Genomics) following the manufacturer's protocol. In brief, a single-cell suspension of cells in 0.4% bovine serum albumin in PBS were added to each channel on the 10X chip. Cells were partitioned with gel beads into emulsion in the Chromium instrument in which cell lysis and barcoded reverse transcription of RNA occurred following amplification. scRNAseq libraries were prepared using the Chromium single cell 3′ library and gel bead kit v3 (10X Genomics). Sequencing was performed on HiSeq platform (Illumina), and the digital expression matrix was generated using the Cell Ranger pipeline (10X Genomics). In total, 1833 (99,466 mean reads per cells) and 1504 (67,463 mean reads per cell) pdgfrb:EGFP cells were sequenced for control and pdgfrb mutant hearts, respectively. After removing low-quality cells, we analyzed 955 and 1126 cells from control and pdgfrb mutant hearts, respectively. From a combined dataset of control and pdgfrb mutant hearts (955+1126 cells), 274 pdgfrb only and 126 pdgfrb; cxcl12b cells were identified and analyzed from 5826 cells from uninjured hearts and 3127 cells from regenerating WT adult hearts with an average of 29,327 and 66,193 reads per cell, respectively. After removing low-quality cells, finally we analyzed 1631 cells from uninjured hearts and 1097 cells from 7 dpa hearts. From combined dataset of uninjured and injured hearts (1631+1097 cells), 511 pdgfrb cells or epicardial/EPDC/mural cells were identified and analyzed using the Seurat R package.

To identify different cell types and find signature genes for each cell type, the R package Seurat (version 3.2.3) was used to analyze the digital expression matrix. Cells with <100 and >2500 unique feature count and >25% mitochondrial expression were removed from further analysis. The Seurat function NormalizeData was used to normalize the raw counts. Variable genes were identified using the FindVariableGenes function. The Seurat ScaleData function was used to scale and center expression values in the dataset for dimensional reduction. Principal component analysis, t-distributed stochastic neighbor embedding and UMAP were used to reduce the dimensions of the data, and the first two dimensions were used in plots. The FindClusters function was later used to cluster the cells. The FindAllMarkers function was used to determine the marker genes for each cluster, which were then used to define cell types. Also, known cell type marker expression was determined across different clusters to assign the cell type to a cluster. We integrated the WT and pdgfrb^−/− dataset together. Also, we integrated uninjured and injured (7 dpa) datasets using pairwise anchors by Seurat. After clustering cells based on the differentially expressed genes, cell types were identified by marker gene expression. For GO term analysis, differentially expressed genes in a cluster were selected with adjusted P-value<0.1. The selected genes were then put in a list form at http://www.zebrafishmine.org for GO term analysis. Enriched GO terms were generally selected after running Holm–Bonferroni or Benjamini–Hochberg corrections with P-value<0.05. All the GO terms were then categorized based on common criteria (e.g. similar expression location, functional similarities, structural similarities, etc.) and their relative percentages were presented as a pie chart. Different plots (e.g. violin plots, dot plots, feature plots) were made following available default Seurat codes. The Subset function was used for isolating and comparing different cell populations either from different clusters (Fig. 3C, Table S2, cluster 6 versus cluster 3 cells), or expressing/not expressing certain genes (Fig. 3D, Table S3, cxcl12b⁺ versus cxcl12b⁻ cells from cluster 6) or from WT versus pdgfrb mutant (from cluster 6, Table S4) or pdgfrb-expressing cells from injured (7 dpa) versus uninjured hearts (combining clusters 5, 6, 10 and 12; Fig. 5B, Table S10). Separate Seurat objects were formed with the isolated cells and merged together to identify differentially expressed genes in the desired group of cells (Fig. 3D, Table S3, cxcl12b⁺ versus cxcl12b⁻ cells). For cxcl12b⁺ cells from cluster 6, cxcl12b expression level >0 was used to subset cxcl12b⁺ cells and cxcl12b≤0 was used to subset cxcl12b⁻ cells (Fig. 3D).

To validate candidate genes found in the mural cell cluster by scRNAseq, qRT-PCR was performed with cDNAs of GFP^high cells compared with GFP^low cells (Fig. S4). For validating upregulated genes in the injured (7 dpa) heart pdgfrb-expressing cells (Fig. 5B), cDNAs from GFP^high cells and GFP^low cells were mixed in a 1:1 ratio for each sample for qRT-PCR (Fig. S10). FACS-sorted pdgfrb:EGFP cells or fli1a:DsRed cells or whole-ventricle samples were collected in Trizol for RNA extraction by the Trizol-chloroform method. The RNAs were precipitated overnight at −20°C in isopropanol and washed with 70% alcohol and dissolved in DEPC-treated water. cDNAs were made from the RNA using the SuperScript^® III First-Strand Synthesis System. qRT-PCR was performed using Applied Biosystems^® SYBR^® Green PCR Master Mix. Gene expression fold changes were normalized to the housekeeping gene rpl13. The mean cycle threshold (Ct) values of each sample triplicate were calculated and 2^−ΔΔCt values were calculated for each sample to determine expression fold changes. See Table S16 for qRT-PCR primers.

We thank Drs Koji Ando and Christer Betsholtz for discussions and communications prior to the submission; Dr Didier Stainier for the pdgfrb:Citrine line; Drs Koji Ando and Hiro Nakajima for the pdgfrb:EGFP line; and Dr Ken Poss for the tcf21:dsRed line. We also thank the TSRI Cellular Imaging, FACS and Single-Cells, Sequencing and CyTOF (SC2) cores.

The peer review history is available online at https://journals.biologists.com/dev/article-lookup/doi/10.1242/dev.199752.