Pigs underwent AR_P1, MI_P28, AR_P1MI_P28, or neither surgical procedure (CTL). Subsequent experiments confirmed that pigs in the AR_P1MI_P28 group, but not MI_P28 animals, completely regenerated the myocardial tissue that was damaged by MI (Zhao et al., 2020; Nakada et al., 2022). snRNAseq was performed with myocardial tissue collected from MI_P28 and AR_P1MI_P28 animals on P30, P35, P42, and P56; from AR_P1 and CTL animals on P28 and P56; from CTL animals on P1, and from fetal pigs.

The cell-cycle-specific Autoencoder (CSA) and clustering pipeline visualized five EC clusters (Figure 1A), each of which expressed elevated levels of the EC-specific surface markers Platelet And Endothelial Cell Adhesion Molecule 1 (PECAM1) (Newman, 1994), Kinase Insert Domain Receptor (KDR) (Terman et al., 1992), and vascular endothelial cadherin (CDH5) (Wu et al., 2021) (Figures 1B–D; Supplementary Figure S1). Three of the clusters were enriched for expression of the vascular-EC (VEC) markers CD34 (Fina et al., 1990) and Plasmalemma Vesicle Associated Protein (PLVAP) (Denzer et al., 2023) (Figures 1E, F), while the other two expressed high levels of the lymphatic-EC (LEC) markers C-C Motif Chemokine Ligand 21 (CCL21) (Johnson and Jackson, 2010), Prospero Homeobox 1 (PROX1) (Wigle et al., 2002), and Lymphatic Vessel Endothelial Hyaluronan Receptor 1 (LYVE1) (Banerji et al., 1999) (Figures 1G–I). Thus, the CSA appeared to separate the global cardiac EC population into three clusters of vascular ECs (VEC1-3) and two clusters of lymphatic ECs (LEC1-2).

Each EC cluster contained cells from multiple experimental groups and time points (Supplementary Figure S2), confirming that our results were not compromised by sampling bias. Further, we verified that the results were not compromised by batch effects. VEC1-3 collectively comprised >95% of ECs in fetal hearts and 59%–81% of ECs in hearts from all nonfetal animal groups (Figure 1J), whereas LEC1-2 included the remaining <5% of ECs in fetal hearts and 19%–41% of ECs in nonfetal hearts. The majority of VECs were clustered in VEC2, regardless of animal groups or time points, with the highest proportions of VEC1 and VEC3 ECs observed in CTL hearts on P1 and in fetal hearts, respectively. Nearly all LECs were clustered in LEC2.

VEC1 ECs comprised 8.7% of ECs in fetal hearts and 14.8% of ECs in CTL hearts on P1, they declined to nearly undetectable levels in CTL hearts on P28, and then expanded to include 3.5%–4.2% of ECs in MI_P28 hearts from P30-P42 and in AR_P1MI_P28 hearts on P30-P35 (Figure 2A). VEC1 was also the only VEC cluster to express elevated levels of each of five markers for mitosis and cytokinesis (Aurora Kinase B [AURKB], marker of proliferation Ki-67 [MKI67], Inner Centromere Protein [INCENP], Survivin [BIRC5], and Borealin [CDCA8] (Li et al., 1998; Vader et al., 2006); Figures 2B–G; Supplementary Figure S3) and was enriched for the expression of genes associated with DNA synthesis and repair (Gene Ontology Annotations, 2023b), chromosome condensation (Gene Ontology Annotations, 2023c), and the gene ontology terms DNA replication (p = 3.86 × 10⁻⁸), mitotic spindle organization (p = 8.71 × 10⁻⁷), and mitotic cytokinesis (p = 2.42 × 10⁻⁶) (Figure 2G; Supplementary Tables S1A, B), as well as numerous transcription factors that have been linked to increases in EC proliferation and angiogenesis (Figure 2G; Table 1). Collectively, these observations are consistent with previous reports that the coronary vasculature undergoes a period of rapid growth and remodeling shortly after birth (Tian et al., 2014; Tan and Lewandowski, 2020) and suggest that a small proportion of cardiac ECs revert to a more proliferative phenotype in response to MI on P28.

Notably, although LEC1 ECs were exceptionally rare (<1.5%) in hearts from all groups and time points, they were most abundant in CTL hearts on P1 and in MI_P28 and AR_P1MI_P28 hearts on P30 (Figure 1I; Supplementary Figure S4), and were enriched for many of the same mitosis/cytokinesis markers, gene ontology terms, and transcription factors that were upregulated in VEC1 ECs (Figure 2G), which suggests that expansion of the LEC1 cluster may have a role in the inflammatory response (Prabhu and Frangogiannis, 2016) and/or contribute to the re-uptake of interstitial fluid after MI. Also, Immunohistological staining of VEC clusters’ markers (PLVAP + PECAM1) (Figure 3) showed that on P30, the percentage of PLVAP + ECs significantly decreased, compared to CTL-P28 level as the baseline, following MI_P28 in both regenerative (p = 3.382 × 10⁻⁴) and non-regenerative heart (p = 3.4995 × 10⁻⁶); then, these percentage returned to the baseline level on P42. Since the percentages of VEC and LEC clusters were complemented, this result also suggested an increase of LEC between 2 and 7 days after MI_P28. This LEC expansion was also suggested by an increase of cycling LEC1 cluster found in the snRNAseq analysis. Besides, trajectory analysis of LEC clusters (Supplementary Figure S4) also showed continuous paths from LEC2 to LEC1; the ordered was from cell-cycle G1/S phase marker Cylin D1 (CCND1) through S-phase marker CCNE2, M-phase marker CCNB2, to proliferation marker MKI67.

Although VEC1 ECs became more common in pig hearts after MI was induced on P28, they tended to be less common (though not significantly) in AR_P1MI_P28 hearts than in MI_P28 hearts from P30-P42 and, consequently, are unlikely to contribute to the enhanced regenerative potential observed in the hearts of AR_P1MI_P28 pigs. Conversely, the VEC3 cluster was significantly more prominent in fetal (30.46%) than CTL hearts (3.9%–2.0% from P1-P56; p < 0.05 at each time point) and in AR_P1MI_P28 hearts than in MI_P28 hearts on P30 (AR_P1MI_P28: 6.33%, MI_P28: 1.6%; p < 0.05) and P35 (AR_P1MI_P28: 6.86%, MI_P28: 1.8%; p < 0.05) (Figure 4A), but was not enriched for the expression of mitosis/cytokinesis markers (Figure 2G; Supplementary Table S2A). Instead, VEC3 ECs expressed elevated levels of genes associated with the EndMT, with the Transforming Growth Factor β (TGFβ) and Bone Morphogenetic Protein (BMP) (Figure 4B; Supplementary Table S2B) pathways (Huang et al., 2022; Gene Ontology Annotations, 2023a; KEGG, 2023), and with ontology terms that are associated with the upregulation of blood-vessel development (Figure 4C; Supplementary Table S2B). Notably, both TGFβ and BMP signaling regulate the EndMT, and although the EndMT appears to contribute to the formation of a number of pathological vascular abnormalities (Huang et al., 2022), all three processes are crucial for embryonic cardiovascular development (Gordon and Blobe, 2008; Morrell et al., 2016; Kovacic et al., 2019). Furthermore, VEC3 ECs expressed elevated levels of Plakophilin 2 (PKP2, fold-change: 24.62, p-value: 5.74 × 10⁻²⁴⁷), which regulates the assembly of intercellular junctions (Bass-Zubek et al., 2008), as well as the adhesion molecules Cadherin 11 (CDH11) and Cadherin 18 (CDH18) (Figures 4D, E; Supplementary Figure S5), while expression of the junctional proteins KDR and CDH5 (Giannotta et al., 2013) was somewhat lower in VEC3 than in VEC1 ECs (Supplementary Figure S1).

The uniqueness of the VEC1 and VEC3 clusters was also evidenced by the results from Monocle pseudotime trajectory analysis (Qiu et al., 2017) of cell-cycle snRNAseq data for cells in all three VEC clusters. The analysis produced a branched structure with the tips of two branches each composed almost exclusively of VEC1 or VEC3 ECs and the third branch, as well as the region where all three branches converged, occupied by primarily VEC2 ECs (Figure 5A). Furthermore, the distribution of cells expressing Cyclin D1 (CCND1), which regulates the G1/S phase transition, and the proliferation markers AURKB and MKI67 traced a pathway between VEC2 and VEC1, but not between VEC2 and VEC3 or VEC3 and VEC1 (Figures 5B–D). Collectively, these observations indicated that VEC2 ECs could transform into either VEC1 or VEC3 ECs, but there was no direct transformational pathway between the VEC1 and VEC3 clusters, and identified a trajectory from VEC2 to VEC1 that recapitulated the onset and progression of cell-cycle activity and proliferation.

Following MI_P28, the vessel density, quantified by the percentage of image PECAM1 area and by the manual counting of vessels in each image (Figure 6), significantly decreased on P30 in non-regenerative hearts; then, on P35 and P42, their vessel densities recovered to the level of CTL-P28 hearts. Meanwhile, in regenerative hearts, the vessel densities did not significantly change, remaining similar to CTL-P28. On the other hand, Marker of Proliferation Ki67 immunofluorescence analysis (Figure 7) showed that seven to 14 days following the MI_P28 injury, endothelial cells in both regenerative and non-regenerative hearts actively proliferate. Before the injury, less than 1% of EC expressed Ki67. Then, on day 7 and 14 after MI_P28, the percentages of Ki67+ EC increased to, and remained between 6% and 8% in the non-regenerative hearts (MI_P28P35 and MI_P28P42, respectively); meanwhile, in the regenerative hearts, % Ki67+ EC increased to 9.88% on P35, then decreased to 5.20% on P42. The difference of Ki67+ EC between the two groups was insignificant at these time points (p = 0.062 and p = 0.467, respectively), and was consistent with the snRNAseq analysis. Then, we compared the expression of Vascular Endothelial Growth Factor A (VEGFA) and Vascular Endothelial Growth Factor B (VEGFB), which are known as the driver genes in angiogenesis involved in wound healing (Johnson and Wilgus, 2014), in the whole-heart snRNAseq data. Supplementary Figure S6 showed that it was unclear whether these two genes expressions were consistently higher in the regenerative (or on-regenerative) groups. This observation may explain why we did not observe significant angiogenesis between the regenerative and non-regenerative hearts.

In Figure 8, the expression of EndMT marker PDGFRA in EC was represented by the intensity of red fluorescence overlapping with green (EC marker PECAM1) fluorescence. A small number (∼1.5%) of EC in the CTL-P28 hearts were VEC3; therefore, the PDGFRA expression in CTL-P28 EC was used as the baseline. Then, following MI_P28 injury, PDGFRA expression in the non-regenerative hearts EC sharply decreased on P30, then gradually increased on P35 and P42; however, the increased expression did not reach beyond the baseline level. Meanwhile, PDGFRA expression sharply increased in the regenerative heart on P30 (p = 7.69 × 10⁻⁴) and P35 (p = 4.71 × 10⁻⁵), then reverted to the baseline level on P42. These trends were consistent with the snRNAseq analysis. Also, surprisingly, PDGFRA immunofluorescence analysis showed that non-endothelial cells surrounding PDGFRA + EC also expressed PDGFRA; these cells could be cardiac fibroblasts and cardiomyocytes, as reported in (Bloomekatz et al., 2017; Yao et al., 2022).

Our previous works demonstrated the upregulation of Hippo/YAP signaling pathway in regenerative pig hearts that underwent myocardial infarction on P1 (Zhang et al., 2020), as well as ERBB4 in AR_P1P28 hearts (Nguyen et al., 2022) at the protein level. Also, Pi3k-Akt signaling pathway was found enriched in AR_P1MI_P28 cardiomyocytes (Nakada et al., 2022; Nguyen et al., 2023b). Therefore, we examined whether the VEC3 cluster may have genes that interact with these pathways (Figure 9A). Thirty genes encoding secretable proteins were found upregulated in VEC3; among them, Fibronectin 1 [FN1] (Wang et al., 2013), Neuregulin 1 [NRG1] (Bersell et al., 2009; Gemberling et al., 2015), Versican [VCAM] (Feng et al., 2023), and Periostin [POSTN] (Kuhn et al., 2007) were known to promote cardiomyocyte proliferation in zebrafish, mice, and rat. We also queried the protein-protein interactions between these 30 genes and their receptors linked to Pi3k-Akt, Hippo/YAP, and ERBB4 signaling pathways in the Biological General Repository for Interaction Datasets (BioGrid) (Oughtred et al., 2021). Figure 9B shows that POSTN, FN1, and VCAM interact with Integrin Alpha 1 [ITGA1] and Integrin Beta 1 [ITGB1], which are receptors of the Pi3k-Akt signaling pathway (Feng et al., 2023). Also, Fibrillin 2 [FBN2], which was upregulated in VEC3, interact with Insulin Receptor [INSR] of Pi3k-Akt signaling pathway. Furthermore, NRG1 interacts with ERBB4, and Bone Morphogenetic Protein 6 [BMP6] interacts with receptors Bone Morphogenetic Protein Receptor Type 1A/B and Type 2 [BMPR1A/B, BMPR2] of the Hippo/Yap signaling pathway. Overall, these secret protein-receptor interactions suggested that VEC3 cluster may promote cardiomyocyte proliferation and overall regeneration.

The pig snRNAseq dataset was obtained from the Gene Expression Omnibus database (accession number GSE185289). The complete dataset included 250,700 cells and was processed as described previously (Nakada et al., 2022; Nguyen et al., 2022; Nguyen et al., 2023b). The number of pig hearts for each group was summarized in Supplementary Table S3. The data were normalized with Seurat v.4 (Hao et al., 2021); then, expression data for 14,753 genes with at least 1,000 unique molecular identifiers were embedded into 10 dimensions with an Autoencoder. The results were visualized via Uniform Manifold Approximation (UMAP), and clustering was performed with the UMAP density-based clustering (dbscan) algorithm (McInnes et al., 2018; Meehan, 2021). Cells in clusters that specifically expressed the EC markers PECAM1 and KDR were considered ECs. A total of 29,826 ECs were identified.

Frozen tissues, from which the snRNAseq datasets (Nakada et al., 2022; Nguyen et al., 2022; Nguyen et al., 2023b) were generated, were cut into 1 cm thick sections, embedded in optimal cutting temperature compound, and preserved in −80°C freezers. All experimental protocols were approved by The Institutional Animal Care and Use Committee (IACUC) of the University of Alabama, Birmingham, and performed in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals (NIH publication No 85-23).

Upregulated genes for each EC cluster were identified according to the following criteria: 1) a cluster p-value (Fisher’s Exact test (Fisher, 1922)) of less than 10⁻⁶, 2) expression by more than 25% of cells in the cluster, and 3) mean expression at least 2-fold greater than among ECs in other clusters. For transcription factors, which are typically expressed in low abundance, criteria for upregulation were modified to 1) a cluster p-value of <10⁻⁶, 2) expression by >10% of cluster cells, and 3) mean expression >2-fold greater than among cells in other clusters. The list of transcription factors was obtained from the Transcriptional Regulatory Relationships Unraveled by Sentence-based Text mining database (TRRUST) (Han et al., 2018) version 2, and the set of upregulated genes for each cluster was analyzed with the DAVID functional annotation tool (Huang da et al., 2009); only terms present in the Gene Ontology (Huntley et al., 2015) and KEGG (Kanehisa et al., 2017) categories were selected. To avoid false positives caused by multihypothesis testing, the selected results were required to have a Benjamini-adjusted (Li and Barber, 2019) p < 0.05.

The immunofluorescence analyses were conducted similar to our previous work in (Nakada et al., 2022). Briefly, samples were cut into transverse sections (thickness: 10 μm), embedded into slides, washed in PBS, and fixed in 4% paraformaldehyde for 20 min at room temperature. Then, the slides were permeabilized at 4°C, 0.1% triton-x for 10 min, and then blocked with Ultra-V block for 7 min. Samples were incubated with mouse anti- PECAM1 at 1:100 (MCA1738, BioRad, United States) or rabbit anti Ki67 (1:100 dilution, ab15580, Abcam, United States), or rabbit antiPDGFRA (1:100 dilution, AF0241, Affinity, United States), or rabbit anti-PLVAP (1:100 dilution, AF0241, Affinity), in PBS containing 10% goat serum (10% PBS)for overnight at 4°C. On the next day, donkey anti-mouse IgG conjugated with FITC at 1:200 dilution and donkey anti-rabbit IgG conjugated with Cy3 at 1; 200 dilution were added for 1 h at room temperature Then, samples were washed after 1:1000 DAPI solution was added for 1 min, mounted in vectorshieldand cover slide for imaging.

The slides were put into an Olympus IX83 microscope to capture images. The image frame size was 1360-width × 1024-height pixels. Exposure time for PECAM1 fluorescence (green) was initially set to 2.8 s, and could be case-by-case adjusted between 2.6 and 3.2 s. Exposure time for Ki67 and PLVAP fluorescence (red) was initially set to 1.2s, and could be adjusted between 0.9 and 1.4 s. Exposure time for PDGFRA fluorescence (red) was always set to 1.2s. DAPI (blue) exposure time was between 12 and 18 μs. Z-stack parameter was set to 0 for all images and channels. For each slide, ten image sets were captured; each image set include three channels: red, green, and blue; then, the three-channels were put together by Matlab to create a ‘merge’ RGB-color image.

Three-channel images (green: PECAM1; red: either Ki67, PDGFRA, or PLVAP; blue: DAPI) were processed using Matlab as in (Nguyen et al., 2022). Multi-thresholding segmentation method (Otsu, 1979) (Matlab multithresh function) was applied to extract image fluorescence regions. The vessel density was computed as the percentage of image area with green fluorescence. EC nuclei were recognized as the DAPI (blue) blocks overlapping with green regions; EC proliferating cells were recognized as the Ki67 red block overlapping with green regions; and EC proliferation marker was computed by the ratio between EC proliferating cells and EC DAPI. EC PDGFRA expression was assessed by the average red intensity overlapping with green fluorescence. For PLVAP + EC analysis, the PECAM1 and PLVAP fluorescence intensity was amplified by 2x before merging into one 3D image channel, PLVAP signal that were not overlapping with PECAM1 was removed for enhanced visualization (PLVAP/DAPI/PECAM1(e) images); then, the number of PLVAP + ECs was manually counted. Statistical comparison among the groups was completed using non-parametric Wilcoxon’s Ranksum test. A p-value < 0.05 indicates statistical significance.