Single-cell transcriptomics defines heterogeneity of epicardial cells and fibroblasts within the infarcted murine heart

In the adult heart, the epicardium becomes activated after injury, contributing to cardiac healing by secretion of paracrine factors. Here, we analyzed by single-cell RNA sequencing combined with RNA in situ hybridization and lineage tracing of Wilms tumor protein 1-positive (WT1+) cells, the cellular composition, location, and hierarchy of epicardial stromal cells (EpiSC) in comparison to activated myocardial fibroblasts/stromal cells in infarcted mouse hearts. We identified 11 transcriptionally distinct EpiSC populations, which can be classified into three groups, each containing a cluster of proliferating cells. Two groups expressed cardiac specification markers and sarcomeric proteins suggestive of cardiomyogenic potential. Transcripts of hypoxia-inducible factor (HIF)-1α and HIF-responsive genes were enriched in EpiSC consistent with an epicardial hypoxic niche. Expression of paracrine factors was not limited to WT1+ cells but was a general feature of activated cardiac stromal cells. Our findings provide the cellular framework by which myocardial ischemia may trigger in EpiSC the formation of cardioprotective/regenerative responses.


Introduction
Myocardial infarction (MI), still the most frequent cause of death in western societies, is associated with massive activation of cardiac fibroblasts, ultimately resulting in excessive accumulation of extracellular matrix (ECM) components that finally impair cardiac function (Tallquist and Molkentin, 2017). During development, the majority of cardiac fibroblasts are derived from the epicardium, which forms the thin outermost epithelial layer of all vertebrate hearts and exhibits extensive R sa d 2 D kk3 C ox4 i2 R g s5 W if1 U b e 2 c S tm n 1 P rc1 To p 2 a P b k P tn C o l5 a 3 T n c T m sb 4 x C a lcb L g a ls7 K r t1 9 K r t1 8 K r t7 M sln S fr p 2 E ln S e le n o p P d lim 3 M fa p 4 D p e p 1 C fh C le c3 b A p o e Ly6 a C ra b p 2 T n fa ip 6 A q p 1 C d 4 4 G stm 1 S p a rcl1 S m o c2 P i1 6 G sn C xcl5 S p p 1 M yl9 Ta g ln A cta 2 G a tm D m kn Epicardial stromal cells (EpiSC) and activated cardiac stromal cells (aCSC) were simultaneously collected from the surface and the myocardium of the isolated perfused heart by applying mild shear forces to the cardiac surface (Owenier et al., 2020) at 5 days post myocardial infarction (MI) (n = 3). CSC were purified from three non-infarcted control hearts 5 days after sham surgery. Mesh purification, low-speed centrifugation, and cell sorting by flow cytometry were performed to remove cardiomyocytes,  Figure 1E and top 10 heat map in Figure 1-figure supplement 3. Remarkably, transcripts of epithelial cell-associated genes and established epicardial genes (Bochmann et al., 2010) were enriched in both EpiSC-1 (Dmkn, Saa3) and Krt7,Krt18,Krt19,Lgals7), while mesenchymal marker genes were primarily expressed in EpiSC-3 (Gsn, Pi16) and Ly6a). Genes coding for ECM proteins were highly expressed in EpiSC-2 (Spp1), EpiSC-3 (Smoc2, Sparcl1), EpiSC-5 (Clec3b), EpiSC-6 (Mfap4, Eln), and Tnc). Genes encoding contractile proteins were preferentially expressed in EpiSC-2 (Acta2, Myl9, Tagln). Wnt pathway-associated gene transcripts were enriched in EpiSC-6 (Sfrp2) and Dkk3). Genes related to the cellular response to interferon characterized Ifit3,Ifit3b,Cxcl10). Finally, expression of genes associated with high cell cycle activity and mitosis was a feature of Top2a,Prc1,Stmn1,Ube2c).

Hierarchy of EpiSC populations
To better define the cellular relationship between the identified EpiSC populations, we performed RNA velocity analysis (La Manno et al., 2018), which can predict and visualize future cell states based on the ratio of unspliced and spliced mRNA read counts. As shown in Figure 1F, this trajectory inference separated the EpiSC populations into three groups: EpiSC-7 and -1 (group I), EpiSC-2, -4, -8, -11 (group II), and EpiSC-6, -5, -3, -9, -10 (group III). Cell cycle scoring showed that each of the three EpiSC groups identified above was equipped with a cell cluster expressing a number of G 2 /M and S phase genes ( Figure 1G), indicating proliferating cells: within group I, a subpopulation of EpiSC-1, within group II, a subpopulation of EpiSC-2, and within group III, EpiSC-9. This location is very similar to the origin of velocity arrows ( Figure 1F).
The number of expressed genes per cell was lower in group III than in groups I and II ( Figure 1H). Since the expressed genes per cell correlate with developmental potential (Gulati et al., 2020), this may indicate that a large fraction of group III cells are terminally differentiated.
To study cell-cell communication mediated by ligand-receptor interactions between the EpiSC populations, we used CellPhoneDB (Efremova et al., 2020). As shown in Figure 1I, the highest numbers of interactions were predicted between groups I and II. Ligand-receptor pairs potentially involved in the interaction between the groups include IGF1, transforming growth factor beta-2 (TGF-b2), and Wnt4 signaling (Supplementary file 3, Figure 1-figure supplement 4). Source data 1. Source data for EpiSC population cell counts summarized in Figure 1C.

Epicardial marker gene expression
The cellular distribution of well-established epicardial progenitor marker genes such as Wt1, Tbx18, Sema3d, Aldh1a2, Gata5, and Tcf21 (Cao and Poss, 2018) is shown in Figure 2A. Wt1, commonly used for lineage tracing studies (Quijada et al., 2020), was found to be highly enriched in group I EpiSC populations 1 and 7, as well as the adjacent population EpiSC-8. Sema3d, Aldh1a2, and Gata5 were also predominantly expressed in group I, with the highest expression levels in EpiSC-7. Tbx18, on the other hand, was broadly distributed in group I and II populations . In contrast, Tcf21 was strongly expressed in group III EpiSC-3, -5, -6, and -9 and showed no overlap with Wt1-expressing populations (Figure 2A).
Since cells in EpiSC-1 showed substantial inhomogeneity of markers (Figure 2A), we carried out a separate clustering analysis for EpiSC-1, including the two adjacent Wt1-expressing populations EpiSC-7 and -8. We identified five subclusters in EpiSC-1 (EpiSC-1.1, -1.2, -1.3, -1.4, -1.5), with no subclustering in EpiSC-7 and -8 ( We also searched for expression of Tgm2, Sema3f, and Cxcl12, which were recently found to mark three functional different epicardial populations in the developing zebrafish heart (Weinberger et al., 2020). We found Tgm2 and Sema3f preferentially expressed in group I populations EpiSC-1 and -7 (Figure 2-figure supplement 2A-B). However, Mylk, which is an additional marker of the Sema3f-expressing zebrafish epicardial population (Weinberger et al., 2020), was primarily expressed in EpiSC-8. Cxcl12, present only in a small cell population in the zebrafish (Weinberger et al., 2020), was rather broadly expressed with enrichment in EpiSC-2 and subcluster EpiSC-1.1 (Figure 2-figure supplement 2C-D). These findings not only demonstrate some degree of evolutionary preservation in epicardial populations from zebrafish to mice, but also reveal considerable differences in marker gene expression patterns and population sizes.

Spatial mRNA expression of EpiSC population identifiers
To determine the specific location of EpiSC populations in the infarcted heart, RNA in situ hybridization of gene transcripts for selected EpiSC populations was carried out ( Figure 2B, upper panel), using Wt1 and other population-specific identifiers ( Figure 2B, lower panel). Msln expression (EpiSC-7) was detected on the outmost layers of the epicardium, consistent with the epithelial signature of EpiSC-7. Expression of Wt1 (EpiSC-7, -1, -8) and Cd44 (EpiSC-4) similarly labeled cells that were localized in the outer part. In contrast, Ifit3 (EpiSC-11) expression was also found in deeper layers of the activated epicardium and was additionally detected in the myocardium, labeling aCSC. Expression of group III population markers Sfrp2 (EpiSC-6), Dkk3 (EpiSC-10), and Pcsk6 (EpiSC-3) was rather homogenously dispersed throughout the epicardium and the aCSC in the infarct border zone. Top2a, chosen as an identifier of the proliferative (sub)populations of all three groups (EpiSC-1.4, EpiSC-2 subpopulation, EpiSC-9), showed a similar expression pattern. This scattered distribution of proliferating cells in the post-MI heart and especially within the epicardial multi-cell layer has already been confirmed in BrdU incorporation experiments (van Wijk et al., 2012).
To explore whether the dispersed Wt1-negative group II and group III cells might derive from the Wt1-expressing populations located in the outmost layer of the activated epicardium, we combined scRNAseq of EpiSC with lineage tracing using tamoxifen-inducible Wt1-targeted reporter (Wt1 CreER-T2 Rosa tdTomato ) mice ( Figure 2C). Transcriptional profiles of 13,373 cells from two mouse hearts 5 days after MI were assigned to previously defined EpiSC populations ( Figure 2D  Expression of cardiomyogenesis-associated genes, HIF-1-responsive genes, and paracrine factors Lineage tracing experiments have shown that cardiomyogenesis can be initiated in vivo in WT1 + epicardial cells when stimulated with thymosin b4 (10). We found thymosin b4 (Tmsb4x) highly expressed in Wt1-expressing EpiSC-8 ( Figure 3A and B), together with BRG1 (Smarca4), a transcription activator for WT1 (Vieira et al., 2017). Interestingly, many key cardiogenic factors were also expressed within the Wt1-expressing EpiSC populations ( Figure 3A, cardiogenic factors I). These include MESP1, which marks early cardiovascular progenitor specification (Chiapparo et al., 2016), as well as GATA4 ( Figure 3B), which in combination with MEF2C and TBX5 has been shown to facilitate cardiomyogenesis (Ieda et al., 2010). Surprisingly, several Hox family members were expressed in EpiSC-8, for example, Hoxa5 ( Figure 3A and B). HOX transcription factors are downstream effectors of retinoic acid signaling (Nolte et al., 2019), which is required for differentiation of cardiac progenitors during heart development (Roux and Zaffran, 2016). Intriguingly, a second set of cardiogenic factors was found in group III populations EpiSC-3, -5, -6, -9, and -10 ( Figure 3A, cardiogenic factors II). This includes MEF2C ( Figure 3B), Nkx-2.5, as well as BMP2 and BMP4, the latter of which are crucial in the regulation of Nkx-2.5 expression and specification of the cardiac lineage (Brown et al., 2004). In addition, we found low expression of genes encoding muscle structural proteins such as Myl2, -4, and -6, Tnnt2, Ttn, and Nebl, which appear to be present especially in EpiSC-7 and EpiSC-9 ( Figure 2C, contractile proteins). The highest expression levels of Notch target genes were found in Wt1-negative EpiSC populations, especially EpiSC-10 ( Figure 3A). This is remarkable, since Notch-activated epicardial-derived WT1 + cells were described as a multipotent cell population with the ability to express cardiac genes (Russell et al., 2011). The observation that EpiSC-7 and -9 express cardiac specification markers and sarcomere proteins is suggestive that these populations have cardiomyogenic potential.
A hypoxia-responsive element in the Wt1 promoter was reported to bind HIF-1a, which is required for WT1 induction (Wagner et al., 2003). We found expression of Hif1a and HIF-1-responsive genes, particularly those associated with glycolysis, to be enriched in group I and group II EpiSC populations ( Figure 3C). These include the lactate-generating enzyme lactate dehydrogenase (Ldha) and the lactate-exporting monocarboxylate transporter 4 (Slc16a3) ( Figure 3C and D). Consistently, the genes with significantly higher expression in group I and II populations vs group III populations (Supplementary file 6) showed enrichment of several GO terms related to the hypoxia/HIF-mediated metabolic switch toward oxygen-saving glycolysis, such as glycolytic process (FDR q-value: 3.05E-02), glucose metabolic process (FDR q-value: 3.22E-02), and pyruvate biosynthetic process (FDR q-value: 2.66E-02). Due to the transient nature of the HIF-1 response that strongly depends on local oxygen levels, the enriched expression of HIF-1-responsive genes in group I and group II in comparison to group III suggests that the populations reside in microenvironments with different oxygen levels within the post-MI epicardium.

Post-MI aCSC in comparison to EpiSC
To directly compare the transcriptional profile of EpiSC to myocardial aCSC, we performed canonical correlation analysis (CCA) space alignment of these two post-MI scRNAseq data sets, together with CSC from uninjured hearts. Generated CCA clusters ( Figure 4A, left panel) are also displayed according to original cell IDs ( Figure 4A, right panels). Average RNA expression levels and differentially expressed marker genes of CCA clusters are listed in Supplementary files 7 and 8, respectively. The contribution of EpiSC, aCSC, and CSC fractions to CCA clusters is summarized in Figure 4B. As can be seen, CCA clusters A, D, and F were mainly composed of EpiSC, indicating that the transcriptional profile of EpiSC-1, -4, -7, and -8 is prevalent in the epicardium ( Figure 4B). On the other hand, CCA clusters B, C, H, I, J, and O were dominated by aCSC. Individual assignments showed that Wt1-negative EpiSC-2, -6, and -9 carried to variable degree expression signatures of genuine aCSC-1, -2, -3, -4, -5, and -10. Despite these similarities, there were multiple significant differences in average gene expression levels when comparing Wt1-negative EpiSC and -4 -3 -2 -1 0    . Wt1-expressing CCA clusters A-D were excluded from this comparative analysis, since Wt1 expression was only found at the heart's surface ( Figure 2B), strongly suggesting that aCSC assigned to these clusters probably represent epicardial contamination in the aCSC fraction (see above). As shown in Figure 4C, expression of epithelial/ epicardial genes was generally higher in EpiSC, while mesenchymal/fibroblast genes were more dominant in aCSC. A distinct fraction of cardiogenic factors set I (see Figure 3A) was preferentially expressed in EpiSC ( Figure 4D) while the opposite was true for cardiogenic factors set II ( Figure 4E). Notably, cardiac contractile proteins ( Figure 4F) and HIF-1-responsive glycolytic enzymes ( Figure 4G) were predominantly expressed within EpiSC. Among the paracrine factors ( Figure 4H), we found multiple chemokines highly enriched in EpiSC.
To further compare cell states in our stromal cell fractions, we used SCENIC (Aibar et al., 2017) for gene regulatory network reconstruction. This tool scores the activity of transcription factors by correlating their expression with the expression of their direct-binding target genes (Aibar et al., 2017). As shown in Figure 4I and J, EpiSC and aCSC showed distinct patterns of network activity, again emphasizing their different cellular identities and functions. As to be expected, both post-MI fractions showed major differences in network activity compared to CSC from uninjured hearts ( Figure 4K).  Figure 4-figure supplement 6) and most likely represents cells of the epicardial monolayer of the uninjured heart. Individual single-cell analysis of this monolayer is technically not feasible, because the number of cells liberated by shear forces from uninjured hearts is only rather small (estimated~5600 cells/heart) and there is shear-independent release of cells (background) of unknown origin (Owenier et al., 2020). CSC-12 was the only Wt1-expressing population in the CSC fraction (Figure 4-figure supplement 4E).

Comparison to CSC from uninjured hearts
As shown in Figure 4A-B, CCA space alignment of the CSC and EpiSC data sets revealed that the majority of CSC correlated with EpiSC group III populations EpiSC-3, -5, and -6 (CCA cluster K, L, M). Direct comparison of EpiSC and CSC within individual CCA clusters again showed that expression of epithelial/epicardial genes was enriched in EpiSC, while transcript levels of the conventional fibroblast marker Gsn were higher in CSC (Figure 4-figure supplement 7A). EpiSC showed a generally higher expression of HIF-1-responsive genes (Figure 4-figure supplement 7B) and proliferation-associated genes (Figure 4-figure supplement 7C), consistent with the notion of increased epicardial hypoxia and induction of epicardial cell proliferation in the infarcted heart. Indeed, direct comparison of the epithelial populations EpiSC-7 and CSC-12 in CCA cluster A revealed that HIF-1-       responsive/glycolysis-associated genes (Eno1, Gale, Ldha, Slc16a3) and proliferation-associated genes (Calm,Cenpe,Ran,Ranbp1,Spc25,Ube2c) are upregulated in post-MI EpiSC-7 (Figure 4figure supplement 7D).

Discussion
In this study, we provide a single-cell landscape of the post-MI epicardium with 11 transcriptionally different cell populations. This amazing degree of cellular heterogeneity is similar to that of myocardial cardiac fibroblasts (Farbehi et al., 2019;Forte et al., 2020). The widely used epicardial lineage marker gene Wt1 (Quijada et al., 2020) was selectively expressed in three populations (EpiSC-1, -7, -8), which were localized on the outer surface of the activated epicardium ( Figure 2B). Other commonly used lineage markers and recently identified epicardial population markers of the developing zebrafish heart (Weinberger et al., 2020) were rather heterogeneously distributed (Figure 2A, Figure 2-figure supplement 2A-B), suggesting that they mark different stages of differentiation, commitment, or activity in the adult heart.
The quality of scRNAseq critically depends on the cell isolation technique, which ideally should preserve the native state of cells as close as possible. To minimally perturb the native expression profile of cardiac stromal cells, we used a recently elaborated perfusion protocol for the simultaneous isolation of EpiSC and aCSC from the same infarcted heart, which is short (8 min) and results in a high yield of viable cells (Owenier et al., 2020). The application of mild shear forces on the cardiac surface by a simple motor-driven device permitted the rather selective removal of the EpiSC fraction (Owenier et al., 2020). As compared to a commonly used mincing protocol (30-40 min), we found the yield of aCSC to be considerably higher with our technique, and this was associated with significant lower induction of immediate early response genes (Owenier et al., 2020). Because of these reasons, our data are difficult to compare with published studies at single-cell resolution of the post-MI heart, which relied on a mincing protocol. These studies have identified either no (Skelly et al., 2018;Gladka et al., 2018;Kretzschmar et al., 2018) or only one (Farbehi et al., 2019;Forte et al., 2020;Cui et al., 2019;Asp et al., 2019) epicardial cell population.
RNA velocity analysis ( Figure 1F) and cell cycle scoring ( Figure 1G) suggest that the EpiSC populations can be classified into three different independent population groups, each containing a cluster of proliferating cells. Marker genes of group I, comprising Wt1-expressing EpiSC-1 and -7, have been reported in healthy adult mouse epicardium obtained by laser capture (Bochmann et al., 2010) and have been exclusively used in previous single-cell studies to identify epicardial cells (Wagner et al., 2003;Bujak et al., 2009; Figure 1F). Surprisingly, group I cells accounted for only 26% of the EpiSC fraction from the injured heart. Furthermore, the expression of several paracrine proangiogenic factors that were previously considered specifically derived from WT1 + epicardial cells  was not limited to group I EpiSC, but was even higher in other EpiSC populations ( Figure 3E). Thus, about two-third of all epicardial cells were Wt1-negative, but they are likely to be also involved in the secretion of paracrine factors. Furthermore, aCSC expressed numerous paracrine factors ( Figure 4H). Therefore, secretion of paracrine factors appears to be a general feature of both epicardial and myocardial stromal cells, which all can contribute to cardioprotection. Group II, consisting of four cell populations, represented 40% of the EpiSC fraction and showed enriched expression of chemokines known to be involved in attraction of monocytes and neutrophils ( Figure 3E and F). The expression of these chemokines was a general feature of EpiSC in comparison to aCSC ( Figure 4H). This finding points to a role of epicardial cells in the modulation of the innate immune response post MI, as was already suggested with regard to adaptive immune regulation during the post-MI recovery phase (Ramjee et al., 2017).
Cells with the transcriptional profiles of group I and II populations were generally prevalent in EpiSC in comparison to aCSC ( Figure 4B). Group I and II populations also showed the highest number of potential ligand-receptor interactions ( Figure 1I). Interestingly, EpiSC-8, characterized by expression of Wt1 and high transcript levels of HOX transcription factors, also highly expressed thymosin b4 ( Figure 3A and B), which has tissue-regenerating properties (Smart et al., 2011). The pronounced expression of cardiogenic factors and contractile proteins in group I ( Figure 3A) and the parallel expression of WT1 and thymosin b4 suggest that this cellular network has cardiogenic potential and was involved in the previously reported formation of cardiomyocytes from WT1 + cells after thymosin b4 stimulation (Smart et al., 2011;Wang et al., 2021). Interestingly, group I and II also shared high expression of HIF-1-responsive glycolytic enzymes ( Figure 3C), which again was a general feature of EpiSC compared to aCSC ( Figure 4G). This is in support of the epicardium being a hypoxic niche, which was based on the immunostaining of dispersed HIF-1a-positive cells in the epicardium (Kimura and Sadek, 2012). Since WT1 expression is HIF-1-dependent (Wagner et al., 2003), this is consistent with the view that epicardial HIF-1 signaling is likely to be an important trigger in the ischemic heart to promote cardioprotection.
Group III comprised five cell populations accounting for 34% of all EpiSC. In contrast to group I, group III EpiSC were found localized throughout the activated epicardium ( Figure 2B) and group III population identifiers also labeled stromal cells within the myocardium ( Figure 2B). This finding is consistent with the transcriptional profile of group III EpiSC, which was quite similar to that of major aCSC and CSC populations ( Figure 4A and B). This demonstrated that group III cells exhibit a fibroblast-like phenotype. In addition, group III EpiSC showed the lowest numbers of predicted cell-cell interactions ( Figure 1I) and expressed a lower number of genes ( Figure 1H), suggesting a less active cell state. Another remarkable feature of group III was the expression of a second set of cardiogenic factors ( Figure 3A). However, these cardiogenic factors were prevalently expressed in aCSC ( Figure 4E), which is consistent with the postulated cardiogenic potential of cardiac fibroblasts (Furtado et al., 2014).
That there are different and genuine cellular identities of EpiSC and aCSC is supported by CCA space alignment analysis and gene regulatory network analysis, which revealed distinct expression of epicardial vs mesenchymal/fibroblast genes ( Figure 4C) and individual patterns in transcription factor activity ( Figure 4I and J), respectively. That a direct comparison of cellular expression signatures between EpiSC and aCSC also revealed many transcriptional similarities was not surprising in view of their close developmental relationship (Doppler et al., 2017).
Histological observations by Zhou et al., 2011 andQuijada et al., 2020 found no evidence for migration of WT1 + cells into the infarcted myocardium even after longer periods of time. That WT1 + cells most likely do not migrate into the infarcted heart appears to be in contrast with a study using lentiviral labeling of epicardial cells after pericardial injection of virus expressing fluorescent proteins (Gittenberger-de Groot et al., 2010). Along the same line we have previously shown, that tracking of epicardial cells after labeling with fluorescently marked nanoparticles revealed migration into the injured heart (Ding et al., 2016). Since group II and III EpiSC constitute the majority of post-MI epicardial cells, it is well conceivable that populations of Wt1-negative cells were preferentially marked by the above-mentioned labeling techniques, and this may explain the reported migration into the injured myocardium. Our histological analysis of the epicardial post-MI multi-cell layer revealed that especially cells of group III populations are located in close proximity to the underlying myocardium ( Figure 2B). In contrast, Wt1-expressing group I cells were located in the outmost part of the activated epicardium. Furthermore, our lineage tracing data suggest that until day 5 post MI, WT1 + cells have not contributed to Wt1-negative group II and group III populations ( Figure 2E). Therefore, group I cells appeared to be spatially restricted to the surface layers, while the major epicardial expansion seemed to be driven by Wt1-negative EpiSC of groups II and III with their respective proliferative cell (sub)populations. It is also conceivable that the fibroblast-like cells in group III represent myocardial stromal cells that have migrated into the epicardial multi-cell layer. The notion of a compartmentalization of the epicardial multi-cell layer is supported by the differential expression of hypoxia/HIF-1-responsive glycolytic genes among EpiSC populations, which is suggestive of an exposition to different local oxygen levels.
In summary, our study explored post-MI epicardial cell heterogeneity in the context of all cardiac stromal cells at an unprecedented cellular resolution. Important epicardial properties in post-MI wound healing/regeneration can now be attributed to specific cell populations. A deeper understanding of adult epicardial hierarchy may help decipher the signaling mechanisms by which individual epicardial cell populations interact to specifically stimulate cardiac repair processes originating in the epicardium.

Key resources table
Continued on next page Mice All animal experiments were performed in accordance with the institutional guidelines on animal care and approved by the Animal Experimental Committee of the local government 'Landesamt fü r Natur, Umwelt und Verbraucherschutz Nordrhein-Westfalen' (reference number 81-02.04.2019. A181). The animal procedures conformed to the guidelines from Directive 2010/63/EU of the European Parliament on the protection of animals used for scientific purposes.

Animal procedures
MI followed by reperfusion was performed as previously described (Bö nner et al., 2012). In brief, mice were anesthetized (isoflurane 1.5%) and the left anterior descending coronary artery (LAD) was ligated for 50 min followed by reperfusion. LAD occlusion was controlled by ST-segment elevation in electrocardiography recordings. Sham control animals underwent the surgical procedure without LAD ligation.
Wt1 CreERT2 -mediated lineage tracing was performed as described previously (Smart et al., 2011). In brief, Wt1 CreERT2 Rosa tdTomato mice received tamoxifen injections (2 mg emulsified in sesame oil; i. p.) 5 and 3 days prior to MI to induce CreERT2 activity.

Isolation of EpiSC, aCSC, and CSC
EpiSC and aCSC at day 5 after MI and control CSC from uninjured hearts at day 5 after sham surgery were isolated as previously described (Owenier et al., 2020). In brief, mice were sacrificed by cervical dislocation and hearts were excised for preparation of the aortic trunk in ice-cold phosphate-buffered saline (PBS). Isolated hearts were immediately cannulated and perfused with PBS (3 min), followed by perfusion with collagenase solution (8 min; 1200 U/ml collagenase CLS II (Biochrom, Berlin, Germany)) in PBS at 37˚C.
EpiSC were simultaneously isolated from post-MI hearts by bathing the heart in its collagenasecontaining coronary effluat while applying mild shear forces to the cardiac surface. The effluat was collected, centrifuged (300 g, 7 min), and cells were resuspended in PBS/2% fetal calf serum (FCS)/1 mM ethylenediaminetetraacetic acid (EDTA). Application of shear force and effluat collection were omitted for uninjured sham control hearts due to the absence of an activated, expanded epicardial layer. aCSC and CSC were isolated by mechanical dissociation of the digested myocardial tissue, followed by resuspension in Dulbecco's modified eagle medium (DMEM)/10% FCS. The cell suspensions were meshed through a 100 mm cell strainer and centrifuged at 55 Âg to separate cardiomyocytes from non-cardiomyocytes. The supernatants were again passed through a 40 mm cell strainer, centrifuged (7 min, 300 Âg), and cell pellets resuspended in PBS/2% FCS/1 mM EDTA. Cells were immediately stained for surface markers and applied to fluorescence-activated cell sorting.

scRNAseq
The sorted single-cell suspensions were directly used for the scRNAseq experiments. scRNAseq analysis was performed by using the 10x Genomics Chromium System (10x Genomics, San Francisco, USA). Cell viability and cell number analysis were performed via trypan blue staining in a Neubauer counting chamber. A total of 2000-20,000 cells, depending on cell availability, were used as input for the single-cell droplet library generation on the 10x Chromium Controller system utilizing the Chromium Single Cell 3' Reagent Kit v2 according to manufacturer's instructions. tdTomato lineage tracing experiments were conducted utilizing the Chromium Single Cell 3' Reagent Kit v3 according to manufacturer's instructions. Sequencing was carried out on a HiSeq 3000 system (Illumina, San Diego, USA) according to manufacturer's instructions with a mean sequencing depth of~90,000 reads/cell for EpiSC and~70,000 reads/cell for aCSC. Differences in sequencing depth were necessary in order to achieve a similar sequencing saturation between all samples of~70%.

Processing of scRNAseq data
Raw sequencing data were processed using the 10x Genomics CellRanger software (v3.0.2) provided by 10x Genomics. Raw BCL-files were demultiplexed and processed to Fastq-files using the Cell-Ranger mkfastq pipeline. Alignment of reads to the mm10 genome and unique molecular identifier (UMI) counting were performed via the CellRanger count pipeline to generate a gene-barcode matrix.
The median of detected genes per cell was 3155 for EpiSC, 3265 for aCSC, and 2241 for CSC. The median of UMI counts per cell was 10,689 for EpiSC, 11,110 for aCSC, and 5879 for CSC. Mapping rates (reads mapped to the genome) were about 89% for EpiSC, 90.9% for aCSC, and 86.7% for CSC.
For tdTomato lineage tracing experiments, a custom reference, consisting of the mm10 genome and the full-length sequence of tdTomato, was generated via CellRanger mkref.

Filtering and clustering of scRNAseq data
Further analyses were carried out with the Seurat v3.0 R package (Butler et al., 2018). Initial quality control consisted of removal of cells with fewer than 200 detected genes as well as removal of genes expressed in less than three cells. Furthermore, cells with a disproportionately high mapping rate to the mitochondrial genome (mitochondrial read percentages > 5.0 for EpiSC and aCSC, >7.5 for CSC) have been removed, as they represent dead or damaged cells. Normalization has been carried out utilizing SCTransform. Biological replicates have been integrated into one data set by identifying pairwise anchors between data sets and using the anchors to harmonize the data sets. Dimensional reduction of the data set was achieved by Principal Component Analysis (PCA) based on identified variable genes and subsequent UMAP embedding. The number of meaningful Principal Components (PCs) was selected by ranking them according to the percentage of variance explained by each PC, plotting them in an 'Elbow Plot' and manually determining the number of PCs that represent the majority of variance in the data set. Cells were clustered using the graph-based clustering approach implemented in Seurat v3.0. Doublet identification was achieved by using the tool DoubletFinder (v2.0.2) (McGinnis et al., 2019) by the generation of artificial doublets, using the PC distance to find each cell's proportion of artificial k nearest neighbors (pANN) and ranking them according to the expected number of doublets. Heat maps were generated using Morpheus (https://software.broadinstitute.org/morpheus).

RNA velocity analysis
The Python software velocyto.py (version 0.17) (La Manno et al., 2018) was run on the EpiSC count matrices and BAM files generated by CellRanger (see above) to predict and visualize future cell

RNA in situ hybridization
In situ detection of selected marker gene expression was performed by RNA in situ hybridization using the RNAscope 2.5 HD Reagent Kit-RED (Advanced Cell Diagnostics, Hayward, USA) (Wang et al., 2012) with RNAscope Probes targeting Msln, Wt1, Cd44, Ifit3, Sfrp2, Pcsk6, Top2a, and Dkk3 mRNA sequences. Fresh frozen hearts (5 days after MI) were cut in 10 mm sections. Fixation and pretreatment of the cryosections were performed according to the manufacturer's instructions. The incubation time for the hybridization of RNAscope 2.5 AMP 5-RED was increased to 120 min to enhance the signal. For evaluation of the probe signal, the hybridized heart sections were examined under a standard bright-field microscope (BX61 Olympus, Hamburg, Germany) using a x20 objective. Images were processed for publication using ImageJ/Fiji (Schindelin et al., 2012).

Statistics
Markers defining each cluster, as well as differential gene expression between different clusters, were calculated using a two-sided Wilcoxon Rank-Sum test that is implemented in Seurat.
. Supplementary file 7. Average gene expression levels in clusters from CCA space alignment of EpiSC, aCSC, and CSC with separation in EpiSC, aCSC, and CSC, including p-values of differential expression.
. Supplementary file 8. Genes with significantly enriched expression among clusters from CCA space alignment of EpiSC, aCSC, and CSC.