Coupled myovascular expansion directs cardiac growth and regeneration

In many organs, heterologous cell types establish specialized microenvironments, or niches, that mediate tissue growth and regeneration. Alterations to niche constituents can affect the efficiency of tissue growth and outcomes after injury, making niches a possible target for therapeutic regeneration (Lane et al., 2014; Wagers, 2012). For organs like the intestines, bone marrow, skin and skeletal muscle, niches are classically centered on a stem cell compartment (Fuchs and Blau, 2020). With a few notable exceptions, niches within organs that lack resident stem cells are less defined. For example, in the developing and regenerating neonatal mouse heart, hypoxic niches have been associated with regionalized growth, but the cellular makeup of these niches is unclear (Kimura et al., 2015).

Studies in zebrafish, salamanders and neonatal mice have established a template for innate heart regeneration through proliferation of spared cardiomyocytes (CMs) (Jopling et al., 2010; Kikuchi et al., 2010; Oberpriller and Oberpriller, 1971; Poss et al., 2002). However, innate heart regeneration is a multicellular process with required contributions from epicardial cells, inflammatory cells and nerves (Aurora et al., 2014; Mahmoud et al., 2015; Wang et al., 2015). Recent work has highlighted a crucial role for the vasculature in heart regeneration (Fan et al., 2019; Fernandez et al., 2018; Liu et al., 2020; Marin-Juez et al., 2016). In zebrafish and neonatal mice, cardiac endothelial cells (CECs) rapidly respond to injury, extending nascent vessels into the wound that ultimately guide CM growth (Das et al., 2019; Marin-Juez et al., 2019, 2016). Lineage tracing studies have demonstrated that these new vessels form by proliferation of CECs (Das et al., 2019; Marin-Juez et al., 2019; Zhao et al., 2014). Functional interference with angiogenic responses in the zebrafish or the neonatal mouse heart is associated with defects in CM proliferation (Das et al., 2019; Marin-Juez et al., 2016). Conversely, overexpression of the master angiogenic factor, vegfaa, is sufficient to induce ectopic cardiac growth in zebrafish, suggesting that a stimulated vasculature instructs cardiac growth (Karra et al., 2018). However, a similar role for VEGFA-stimulated CECs in the mammalian heart has yet to be shown.

Here, we spatiotemporally model CM and CEC proliferation in the neonatal mouse heart to investigate the mechanisms underlying cardiac growth and regeneration. We find that CM and CEC proliferation are tightly coupled during postnatal cardiac growth. With cryoinjury (CI), a model of incomplete regeneration, this coupling is spatially impaired. We demonstrate that coupled myovascular expansion is dependent on VEGFA signaling to endothelial VEGFR2 and restoration of coupling after CI through exogenous Vegfa can enhance regeneration. Similarly, we find that genetically-determined levels of circulating VEGFA are associated with higher myocardial mass in humans, suggesting that VEGFA also regulates human myocardial growth. Together, these data demonstrate that coupled expansion of CECs and CMs within a myovascular niche regulates cardiac growth and regeneration.

Neonatal mice are able to regenerate their hearts after injury during the first few days of life (Porrello et al., 2011). The loss of regenerative capacity in the mouse heart is coincident with a developmental decline in the rate of CM proliferation. Although the kinetics of CM proliferation in neonatal mice have been well-documented, dynamics of other cardiac cell types are not well characterized. To better define the relationship of CMs and CECs during the regenerative window of the neonatal mouse heart, we assayed CM and CEC cycling kinetics at various time points following an EdU pulse (Fig. 1A-H). Cardiac sections at different developmental time points were stained for EdU incorporation, along with PCM1 and Erg to specifically mark CM and CEC nuclei, respectively (Alkass et al., 2015; Bergmann et al., 2009; Bergmann et al., 2011; Das et al., 2019). We developed customized image segmentation routines to objectively quantify large numbers of CECs and CMs from cardiac sections (3120±734 CECs per heart and 2328±751 CMs per heart; mean±s.d.). We noted that the relative density of CECs increased between postnatal day (P)1 and P10, whereas the overall density of CMs decreased (Fig. 1E,G). Consistent with previous reports, we found the percentage of EdU⁺ CMs to sharply decline over the first 10 days of life, with a second peak occurring at P5 (Fig. 1C,F). The increase in EdU⁺ CMs at P5 likely coincides with a terminal round of DNA synthesis and binucleation of CMs (Alkass et al., 2015; Soonpaa et al., 1996). EdU incorporation by CECs also declines from P1 to P10 (Fig. 1B,D), with a second peak at P7. Although CM and CEC kinetics differ with regards to the timing of this second peak, the overall trends of their kinetics parallel each other. To determine the strength of this relationship, we compared rates of EdU⁺ CECs and CMs for individual hearts and found them to be correlated (R=0.59, P=0.0007) (Fig. 1H). Above a threshold of ∼8% EdU⁺ CECs, each 1% increase in EdU⁺ CMs is associated with a 1.03±0.27% (P=0.0007) increase in EdU⁺ CECs.

Based on the strong temporal association of CEC and CM proliferation, we next sought to spatially relate CEC and CM cycling. Within thick sections of hearts from P4 mice pulsed with EdU, we observed numerous examples of EdU⁺ CMs adjacent to EdU⁺ CECs (Fig. 1I). To quantify this observation, we assigned coordinates to each CEC and CM nucleus and computed pairwise distances for every CEC and CM. We used this distance information to deconvolve overlapping microenvironments by mapping the position of every CEC relative to an EdU⁺ CM or an EdU⁻ CM as a function of distance, resulting in ‘pseudodistances’ (Fig. 1J). When comparing the density of EdU⁺ CECs relative to an EdU⁺ CM or EdU⁻ CM, we found the percentage of EdU⁺ CECs to be enriched around EdU⁺ CMs compared with EdU⁻ CMs until ∼7 µm (P=0.04 at 7 µm, two-sided Z-test) (Fig. 1K), providing evidence for coupled myovascular expansion during physiologic growth.

Unlike other models of injury to the neonatal mouse heart, CI results in incomplete or inefficient regeneration (Bakovic et al., 2020; Darehzereshki et al., 2015; Polizzotti et al., 2016, 2015). Detailed analyses of CM proliferation kinetics after CI have demonstrated that CM proliferation occurs in this model, but ostensibly at levels lower than that of the uninjured neonatal heart (Darehzereshki et al., 2015; Polizzotti et al., 2015). Thus, we chose to evaluate myovascular expansion following CI at P1 as a model of inefficient regeneration, in a regeneration-competent context (Fig. 2A,B). We first profiled CM and CEC kinetics, considering an average of 2407±891 CECs and 2153±1495 CMs within the border zone for each heart (mean±s.d.). We found that CM and CEC kinetics generally follow the same trend as in the uninjured heart, with an ∼50% decrease in proliferation indices for both cell types over the first 10 days of life (Fig. 2C-E,G). Similarly, CEC and CM density within the border zone had a similar trend to the uninjured heart (Fig. 2F,H). When we compared rates of CM and CEC proliferation in individual hearts, we once again noted that CM and CEC proliferation rates are correlated with ∼1% increase in the percentage of EdU⁺ CECs for every percent increase of EdU⁺ CMs (R=0.54, P=0.008) (Fig. 2I). However, when we performed a spatial analysis of CECs and CMs within the border zone, we did not detect enrichment of EdU⁺ CECs in the immediate 7 µm vicinity of EdU⁺ CMs (Fig. 2J-L). In fact, EdU⁺ CECs may even be depleted around EdU⁺ CMs in the border zone after CI, suggesting that CEC and CM proliferation may not be efficiently coupled after CI.

To better understand the molecular mediators of myovascular coupling during growth and injury, we performed single cell RNA-sequencing (scRNA-seq) using border zones of P7 hearts that were cryoinjured at P1. We profiled 1721 cells and identified nine clusters of cells carrying markers of CECs, CMs, fibroblasts and inflammatory cells (Table S1; Fig. 3A). Clusters 1, 4 and 6 were notable for cells with CEC markers, such as Fabp4, Pecam1, Erg and Vegfr2 (also known as Kdr) (Table S1). Compared with the other CEC clusters, cluster 4 had increased expression of numerous proliferative markers, including Ki67 (Mki67), Prc1, Ccna2 and Ccnb2 (Table S2). A similar cluster of highly proliferative CECs has been also reported by other groups after coronary ligation in neonatal and adult mice (Wang et al., 2020; Wu et al., 2020). In the neonatal mouse heart, VEGFR2 primarily marks cardiac microvascular cells and these cells have recently been implicated in revascularization following injury (Das et al., 2019; Kivelä et al., 2019). Evaluation of sections from neonatal hearts under physiologic growth conditions revealed EdU⁺ CMs within 7 µm of EdU⁺ VEGFR2⁺ CECs that was less common after injury (Fig. 3B,C), thus supporting a role for VEGFR2⁺ microvascular CECs as contributing to a myovascular niche of growth.

For many organs, parenchymal cells secrete angiogenic factors that enable matching of vascular supply to organ size (Rafii et al., 2016). Based on our previous work linking vegfaa overexpression to ectopic cardiomyogenesis in the zebrafish heart, we assayed Vegfa expression during growth and regeneration by quantitative single molecule fluorescent in situ hybridization (Erben and Buonanno, 2019; Karra et al., 2018). We found that, under physiologic growth conditions, Vegfa is expressed in CMs, with a sharp decline in expression from P3 to P10 (Fig. S1; Fig. 3D,F). To determine whether Vegfa levels may be contributing to CEC dynamics after CI, we assayed Vegfa expression in the border zone after CI. Compared with uninjured hearts, Vegfa expression in the border zone at P3 was upregulated by ∼60% (Fig. S1B; Fig. 3D-F). However, following CEC cycling rates after CI (Fig. 2), Vegfa expression markedly declined from P3 to P10 (Fig. S1B; Fig. 3E,F). Together, these data support a dynamic role for myocardial VEGFA to endothelial VEGFR2 signaling as a regulator of the myovascular expansion during growth and regeneration.

Based on our scRNA-seq experiments indicating that many proliferating CECs are Vegfr2⁺, we hypothesized that VEGFR2 signaling is a crucial mediator of myovascular growth in the early neonatal period. To conditionally delete Vegfr2 from CECs, we crossed Vegfr2^flox/flox mice to Cdh5-CreER^T2 mice to generate Cdh5-CreER^T2; Vegfr2^flox/flox (Vegfr2^ΔEC) and Vegfr2^flox/flox (Vegfr2^WT) mice (Hooper et al., 2009; Sorensen et al., 2009). We verified loss of VEGFR2 from CECs by immunostaining for VEGFR2 in Vegfr2^ΔEC mice (Fig. S2A,B). Evaluation of CECs in Vegfr2^ΔEC mice revealed fewer CECs than in Vegfr2^WT animals and ∼80% reduction in CEC proliferation (Fig. 4A-C). To determine how CM growth is affected by the absence of Vegfr2 from CECs, we next assayed EdU incorporation by CMs. Although CM numbers were largely preserved, CM cycling was also attenuated by ∼80% in Vegfr2^ΔEC mice compared with Vegfr2^WT mice (Fig. 4D-F). We next evaluated the effect of Vegfr2 deletion from CECs on myovascular growth after CI. Like our results during early neonatal growth, relative numbers of CECs were decreased whereas relative CM numbers were preserved in Vegfr2^ΔEC hearts (Fig. 4I,L). We identified defects in both CEC and CM incorporation of EdU after injury in Vegfr2^ΔEC mice (Fig. 4G,H,J,K). Specifically, we observed Vegfr2^ΔEC hearts to have ∼80% decrease in cycling CECs after injury and a 40% decrease in cycling CMs. Consistent with our hypothesis, CM expansion during growth and regeneration is dependent on intact Vegfr2 in CECs.

Our expression data suggest a role for myocardial Vegfa in myovascular coupling during growth and after injury (Fig. S1; Fig. 3D-F). To functionally test this concept, we obtained a well-described antibody, B20-4.1.1 (anti-VEGFA), that binds murine VEGFA and prevents interaction with its receptors (Fig. S2C,D) (Liang et al., 2006). Along with decreased CEC density, treatment of neonatal mice with anti-VEGFA decreased CEC proliferation by ∼70% during early neonatal growth and by about ∼90% after injury, indicating a dependency of CEC proliferation on VEGFA (Fig. 5A-C,G-I). Defects in CEC proliferation following anti-VEGFA treatment were accompanied by decreases in the fraction of EdU⁺ CMs during growth and after injury (Fig. 5D-F,J-L), further supporting the need for myovascular coupling during growth and regeneration.

Previous work has described a crucial role for tissue hypoxia as a regulator of CM proliferation in both zebrafish and mice (Jopling et al., 2012; Kimura et al., 2015; Puente et al., 2014). Remarkably, hypoxic pre-conditioning of adult mice results in cardiac growth and enhanced regenerative capacity (Nakada et al., 2017). Mechanistically, hypoxia decreases levels of reactive oxygen species, enabling CM cell cycle-reentry (Puente et al., 2014). As the decreased vascularity in Vegfr2^ΔEC mice and anti-VEGFA-treated mice (Fig. S2B,D) might be expected to result in tissue hypoxia, we examined hypoxia in CMs under these conditions. Pimonidazole is a 2-nitroimidazole used to identify tissue hypoxia based on its formation of stable adducts in the presence of low oxygen tension (Miller et al., 1989; Raleigh et al., 1998). These adducts can be detected by immunofluorescence with the intensity of staining directly proportional to the level of hypoxia. We administered pimonidazole before harvest of hearts from Vegfr2^ΔEC mice, Vegfr2^WT mice and mice injected with anti-VEGFA or vehicle. We noted ∼75% increase in the intensity of pimonidazole uptake by CMs of Vegfr2^ΔEC mice compared with Vegfr2^WT mice and an almost 2-fold increase in mice treated with anti-VEGFA compared with mice treated with vehicle (Fig. S3). In addition, we noted more Hif1α staining in Vegfr2^ΔEC mice and mice treated with anti-VEGFA compared with control animals (Fig. S3) Together, these results demonstrate that hearts of Vegfr2^ΔEC mice and anti-VEGFA-treated mice are hypoxic and suggest that CECs are a required mediator for CM proliferation in response to tissue hypoxia.

Based on the impaired spatial myovascular coupling of CECs and CMs after CI (Fig. 2J-L) and the decrease of Vegfa in the border zone after CI (Fig. 3E,F), we hypothesized that increasing Vegfa levels within the border zone might enhance the efficiency of regenerative growth following CI. To test our hypothesis, we generated adeno-associated virus to overexpress Vegfa (AAV-Vegfa) or GFP (AAV-GFP). We then cryoinjured mice at P1 and injected ∼1.2×10¹⁰ viral genomes of AAV-Vegfa or AAV-GFP into the border zone immediately after injury (Fig. S2E-G). Gross examination of hearts injected with a resin to opacify them at P21 identified increased vascularity at the site of injury of AAV-Vegfa hearts compared with AAV-GFP hearts (Fig. S2H). In addition, we noted that AAV-Vegfa hearts had ∼50% less scarring of the left ventricle and better cardiac function after CI (Fig. 6A-C), all suggestive of enhanced regeneration. Compared with animals treated with AAV-GFP, mice treated with AAV-Vegfa had a significantly higher CEC density in the border zone and a trend towards more Ki67⁺ CECs (8.62±2.59% versus 5.54±1.23%, P=0.31) (Fig. 6D,E). Importantly, this was accompanied by a more than 2-fold increase of Ki67⁺ CMs in the border zone compared with hearts treated with AAV-GFP (Fig. 6F,G). Finally, we evaluated sarcomere morphology and α-SMA expression, two markers of CM maturity (Chen et al., 2021; Porrello et al., 2011). Compared with uninjured hearts and hearts treated with AAV-GFP, AAV-Vegfa-treated hearts were notable for CMs with sarcomeric staining on the periphery of the cell, similar to previous descriptions of sarcomere disassembly (Fig. S4A) (Porrello et al., 2011). Also consistent with a less mature CM phenotype in AAV-Vegfa hearts, CMs were more likely to express the dedifferentiation marker αSMA after AAV-Vegfa treatment (Fig. S4B-D) (Chen et al., 2021). Together, these data indicate that exogenous Vegfa in the border zone can enhance regenerative responses after neonatal CI by increasing CECs, promoting CM dedifferentiation and proliferation, reducing scarring and restoring ventricular function.

Our results after CI in neonatal mice suggest exogenous Vegfa could be used therapeutically to promote cardiac growth and regeneration. Indeed, exogenous Vegfa has been shown to stimulate cardiac repair after experimental infarction many times in adult mammals, even prompting several human clinical trials in patients with ischemic heart disease (Carlsson et al., 2018; Ferrarini et al., 2006; Janavel et al., 2006; Lin et al., 2012; Pearlman et al., 1995; Zangi et al., 2013). However, clinical trials have largely failed to improve outcomes, possibly because of inefficient delivery methods that did not sufficiently raise local Vegfa levels (Oh and Ishikawa, 2019; Taimeh et al., 2013).

Based on its potential as a regenerative factor, we sought to conceptually determine whether VEGFA might regulate cardiac growth in humans using Mendelian randomization (MR). MR is an epidemiologic technique that uses genetic variation to infer causality from observational studies (Emdin et al., 2017; Gray and Wheatley, 1991; Schmidt et al., 2020). Traditional observational studies that associate factors with an outcome are prone to confounding and cannot differentiate causation from reverse causation. However, because innate factors tend to be independent of confounding variables, MR studies are used to point towards causality, and have been used to successfully predict the results of clinical trials (Ference, 2017).

We performed two-sample cis MR analysis with genetically predicted circulating VEGFA concentration as the exposure variable and traits of cardiac structure and function as the outcomes (Fig. 7) (Davey Smith and Hemani, 2014). Under this framework, genetic variants acting in cis that associate with protein abundance (protein quantitative trait loci, pQTL) are used as instrumental variables to estimate the unconfounded causal effects of the protein on the outcomes of interest (Fig. 7A). Genetic variants within the VEGFA gene region were selected from genetic association summary statistics from a genome-wide association study (GWAS) of directly measured circulating VEGFA concentration in over 30,000 individuals (Fig. 7B, Table S3) (Folkersen et al., 2020). Variant association statistics for six parameters of left ventricular (LV) structure and function were then extracted from a GWAS of cardiac magnetic resonance imaging measures in the UK Biobank (Aung et al., 2019). A summary of the GWAS is presented in Table S4. We applied these instruments, using a two-sample MR model that accounts for partial correlation between instruments, to estimate the effect of genetically predicted circulating VEGFA on LV phenotypes (Schmidt et al., 2020). We did not detect a link between genetically predicted VEGFA levels and LV volumes or ejection fraction, but found that higher genetically predicted VEGFA levels resulted in increased LV mass (at type I error rate=0.05/6) (Fig. 7C,D). In total, we found an estimated 1.02 g (95% confidence interval=0.33-1.70 g, P=0.004) increase in LV mass (LVM) per doubling of the circulating VEGFA concentration. Because there is a lack of consensus on how genetic instruments are selected and MR modeling, we performed a sensitivity analysis. Across 120 separate analyses with varying thresholds for instrument selection and four different MR modeling approaches, 119 of these analyses indicate that higher levels of VEGFA are robustly associated with higher LVM (Fig. 7E). These results are consistent with the increased angiogenesis, increased heart size and better cardiac function previously reported in mice that have genetically higher levels of VEFGA (Marneros, 2018). Thus, VEGFA is likely to promote cardiac growth in humans, a finding of therapeutic relevance.

We profiled CEC and CM dynamics to determine that CEC and CM cycling are spatiotemporally coupled in the neonatal mouse heart. Uncoupling CECs from CMs, by deletion of Vegfr2 from CECs or with an anti-VEGFA antibody, dramatically decreases CM proliferation (Figs 4 and 5). Conversely, improving CEC and CM coupling after CI by increasing local Vegfa increases regeneration after CI (Fig. 6). Analogously, myocardial mass of the human heart increases with higher VEGFA levels (Fig. 7). Together, our work demonstrates that the efficiency of cardiac growth is dependent upon myovascular interactions, a finding with strong translational relevance. Approximately 50% of patients with systolic heart failure have flow-limiting coronary artery disease (Benjamin et al., 2018). Even in so-called ‘non-ischemic’ cardiomyopathies, cardiac perfusion is impaired due to microvascular disease and capillary rarefaction (Abraham et al., 2000; Drakos et al., 2010; Mosseri et al., 1991; Parodi et al., 1993; Tsagalou et al., 2008). Thus, as methods to promote proliferation of adult mammalian CMs move towards therapeutics, mechanical or molecular revascularization approaches may be a required adjunct for these approaches to be maximally efficacious.

A key concept suggested by our work is that of growth niches, in which cycling CECs establish inductive microenvironments to promote CM proliferation (Figs 1 and 2). We find that cycling CECs are statistically enriched within 7 µm, or one to two cell lengths, of cycling CMs during physiologic growth. Because we used pairwise comparisons of CEC and CM distances to deconvolve overlapping niches, our spatial analysis likely underestimates the expanse of the niche surrounding cycling CMs. However, our results are remarkably well-aligned with reported increases of CM proliferation within three cell lengths of sprouting angiogenesis in cultured fetal heart sections (Miao et al., 2020). Niches for CM expansion in the postnatal mammalian heart have been previously identified around hypoxic regions, with hypoxia even being able to induce CM proliferation in the adult mouse heart (Kimura et al., 2015; Nakada et al., 2017). Mechanistically, hypoxia exerts its effect on CM proliferation cell autonomously, by modulating reactive oxygen species and the DNA damage response (Puente et al., 2014). Our work suggests another element to the phenomenon of hypoxia-induced CM proliferation. We were able to induce hypoxia in the neonatal heart through two different approaches, genetic deletion of Vegfr2 from CECs and with an anti-VEGFA antibody (Fig. S3). However, even though hypoxia would be predicted to increase CM proliferation, we observed decreased CM cycling in these models. Thus, hypoxia-mediated CM expansion is likely to involve both CM-specific effects but to also depend on CECs within the hypoxic niche. Although hypoxia may be one approach to stimulate heart regeneration, therapeutic hypoxia for patients with heart failure may be challenging. A better understanding of the signaling milieu of the hypoxic niche could lead to alternate approaches for stimulating hypoxia-mediated regeneration.

Because depletion of CECs, either through deletion of Vegfr2 from CECs or by administration of anti-VEGFA, leads to defects in CM proliferation, our work suggests that CECs promote CM cycling (Figs 4 and 5). Mechanistically, CEC expansion may simply lead to new conduits for growth factors or other cell types that promote cardiac growth. However, multiple lines of evidence also point to a direct role for CECs to influence CM proliferation through angiocrines. Indeed, work in zebrafish and neonatal mice support this concept. In zebrafish, inhibition of Notch in CECs results in reduced secretion of Wnt inhibitors that influence CM proliferation (Zhao et al., 2019). More directly, deletion of Igf2 from CECs in developing mice abrogates CM proliferation after neonatal cardiac injury (Shen et al., 2020). Interestingly, deletion of Igf2 from CECs does not affect physiologic CM cycling, raising the intriguing possibility that different angiocrines modulate myocardial expansion during physiologic growth and after injury. Finally, different CEC subsets may contribute different sets of angiocrines, as recent work suggests that Reln from lymphatic endothelial cells has direct effects on CM proliferation (Liu et al., 2020). Future work to identify the complement of angiocrines that instruct CM proliferation in different developmental and injury contexts is needed.

Our work also suggests that restoring myovascular coupling can enhance regenerative capacity. Unlike physiologic growth, CEC cycling is not enriched around cycling CMs in the border zone after CI, a model of incomplete regeneration. Functionally increasing CEC density through overexpression of the master angiokine Vegfa enhances CM cycling and reduces scarring after CI, both signs of enhanced regeneration (Fig. 6). Taken together with our previous work showing that vegfaa overexpression can induce ectopic heart regeneration in zebrafish (Karra et al., 2018), exogenous Vegfa may be one approach to augment regeneration. Indeed, our MR study indicates that VEGFA levels are causally linked to human myocardial mass and predict that more VEGFA can increase human cardiac growth (Fig. 7). Several clinical trials of Vegfa are now underway to revisit whether exogenous VEGFA can be used therapeutically to treat human cardiovascular disease (NCT03409627, NCT03370887 and NCT04125732). Compared with previous clinical trials, these studies use newer delivery methods such as modified RNAs, newer generation viral vectors and more efficient plasmid delivery systems. Based partly on our findings, we would expect that markers of increased CECs, such as improvements in myocardial perfusion, will predict VEGFA-based treatment effects on cardiac mass and function. A second aspect of our work that might be informative to translational work is that the site of Vegfa delivery could be crucial to where growth or recovery can occur. For example, misexpression of vegfaa in the heart of zebrafish after apical resection impairs regeneration, with cardiac growth occurring away from the site of injury (Karra et al., 2018). By contrast, in this work, targeted Vegfa overexpression in the border zone improves regeneration in neonatal mice. Future clinical trials of Vegfa will need to account for the site of Vegfa overexpression when determining effects on cardiac structure and function.

We acknowledge several limitations to our work. First, although we focus on the coupling of CECs and CMs during growth and regeneration, additional cell types that also regulate cardiac growth are likely to be present within niches. Defining niche constituents may be essential to efficiently instruct cardiac growth. Second, although we demonstrate a crucial role of VEGFA-VEGFR2 signaling to coordinate CEC and CM expansion, VEGFR2 broadly marks the cardiac microvasculature. Future work to refine which microvascular cells instruct CM cycling could enable more targeted regenerative approaches via cell transplantation or specific signaling cues. Finally, although our MR studies point to a therapeutic benefit to using VEGFA to stimulate cardiac growth, our MR studies model a lifetime of exposure to VEGFA and do not have the resolution to confirm that increased myocardial mass is the result of CM expansion. Additional work to determine whether VEGFA promotes cardiac growth across developmental stages and to directly link VEGFA to human CM proliferation in vivo are needed to establish VEGFA as a bona fide regenerative factor in humans.

In summary, here, we present work showing that coordinated growth of CECs and CMs guides postnatal cardiac growth and regeneration. We speculate that a better understanding of the signaling milieu within a myovascular niche can inform approaches for heart regeneration.

CD-1 mice (Charles River Labs) were used for profiling myovascular kinetics during growth and regeneration. CD-1 mice were also used for anti-VEGFA antibody and AAV experiments. Cdh5-CreER^T2 and Vegfr2^flox/flox strains have been previously described (Hooper et al., 2009; Sorensen et al., 2009). Cdh5-CreER^T2; Vegfr2^flox/flox and Vegfr2^flox/flox mice were treated with 100 µg of tamoxifen intraperitoneally at P0 and P1. Complete recombination was verified by immunostaining for VEGFR2 (Fig. S2). For EdU incorporation studies, 0.5 mg EdU was given intraperitoneally 4 h before tissue harvest. Cryoinjuries were performed on P1 pups as previously described (Bakovic et al., 2020; Polizzotti et al., 2015).

The number of animals used in these studies is specified in the figure legends. For mouse experiments, we included: (1) animals confirmed to have Vegfr2 deletion by immunostaining for VEGFR2; (2) animals with weight loss after anti-VEGFA treatment; or (3) animals confirmed to have Vegfa overexpression by in situ hybridization after injection with AAV-Vegfa. For assays without computational scoring, readers were blinded to experimental groups. All animal protocols were approved by the Institutional Animal Care and Use Committee at Duke University (NC, USA).

At the time of tissue harvest, hearts were perfused with KCl and then 4% paraformaldehyde. Hearts were then immersed in 30% sucrose overnight and embedded in Tissue Freezing Media for cryosectioning.

Cryosections were blocked with PBST (PBS with 0.1% Tween-20) containing 10% newborn calf serum and 1% DMSO and incubated overnight with primary antibodies at 4°C. Cryosections were then washed with PBST and incubated with secondary antibodies and DAPI (100 ng/ml). For EdU detection, cryosections were incubated in EdU staining solution (100 mM Tris-HCl, 1 mM CuSO₄, 10 mM Azide and 50 mM ascorbic acid in PBS) for 10 min. Primary antibodies used for this study included: anti-PCM1 (Sigma-Aldrich, HPA023370, 1:100); anti-PCM1 (Santa Cruz Biotechnology, sc-398365, 1:100); anti-Erg (Abcam, ab92513, 1:25); anti-CD31 (BD Biosciences, 553370, 1:100), anti-VEGFR2 (R&D Systems, AF644, 1:100), anti-Ki67 (Thermo Fisher Scientific, 4-5698-82, 1:100), anti-actinin (Cell Signaling Technology, 6487P, 1:100), anti-Tnnt (Developmental Studies Hybridoma Bank, CT3, 1:25), anti-αSMA (Abcam, ab5694, 1:200) and anti-HIF1α (Cell Signaling Technology, AF1935, 1:100). Secondary antibodies and azides were conjugated to Alexa-488 (Invitrogen, A21200, A11034 and A21208, 1:200), Alexa-594 (Invitrogen, A11032, A11037 and A11007, 1:200), Alexa-633 (Invitrogen, A21052 and A21072, 1:200) or Alexa-647 (Invitrogen, A21472 and A21469, 1:200). To label CM nuclei, cryosections underwent heat induced epitope retrieval with citrate buffer before immunostaining. Stained slides were imaged on a Zeiss AxioImager M1 epifluorescent microscope, a Zeiss CSU-X1 spinning disk confocal microscope or a Zeiss LSM 510 confocal microscope. For physiologic growth experiments, 9-16 non-overlapping images (40×) of the left ventricle were obtained for each section. For CI experiments, the border zone was defined as the entire region within one to two 40× fields of view along the entirety of the injury plane. For each heart, the three largest sections were imaged.

We developed several customized image segmentation pipelines using CellProfiler and Ilastik for automated scoring EdU⁺ CECs and CMs (Berg et al., 2019; Lamprecht et al., 2007). For CEC quantification, grayscale images for each channel were processed with rolling-ball background subtraction and used as input into a CellProfiler routine that: (1) identified Erg⁺ nuclei based on Erg staining and DAPI intensity; (2) identified EdU⁺ nuclei based on EdU staining and DAPI staining; and (3) identified EdU⁺ Erg⁺ nuclei based on the presence of an EdU object within an Erg object.

For CMs, grayscale images from individual channels were obtained using a Zeiss CSU-X1 spinning disk confocal microscope. Images were pre-processed to generate 32-bit grayscale images and to create a set of images with merged PCM1 and DAPI channels. The DAPI image and the merged PCM1-DAPI images were then input into machine learning routines to generate probability maps for DAPI⁺ nuclei and PCM1⁺/DAPI⁺ nuclei, using Ilastik. Machine learning algorithms were trained using images from our physiologic growth experiments. Probability maps and grayscale images were used as input for a CellProfiler pipeline that: (1) identified nuclei based on a DAPI probability map; (2) filtered nuclei for CMs based on the mean PCM1⁺/DAPI⁺ pixel probability; (3) identified EdU⁺ nuclei based on EdU staining and DAPI staining; and (4) identified EdU+ CM nuclei based on the presence of an EdU object within a CM nucleus. CellProfiler output data was tabulated using the dplyr package in R (https://CRAN.R-project.org/package=dplyr).

For proximity mapping, z-stack images of P4 hearts were obtained using a Zeiss LSM 510 confocal microscope. CM and CEC nuclei were segmented using the Spots tool along with manual refinement, assigned coordinates and categorized as EdU⁺ or EdU⁻ using Imaris (Bitplane). Coordinate files were then used to determine the pairwise distance of each CM and each CEC nucleus. Distances were used to compute the density of EdU⁺ CECs as a function of distance from each PCM1⁺ nucleus using the dplyr and plotly packages in R (https://CRAN.R-project.org/package=plotly; https://CRAN.R-project.org/package=dplyr).

RNAscope Probe Mm-Vegfa-O2 (Lot 16197A), with predicted reactivity against all known murine Vegfa variants, was hybridized against flash-frozen cryosections using the manufacturer's protocol (Advanced Cell Diagnostics). Images of sections were quantified using a CellProfiler routine adapted from Erben and Buonanno (2019). Briefly, the DAPI channel was used to identify nuclei objects and the number of Vegfa spots adjacent to each nucleus was counted.

The left ventricles of cryoinjured hearts were collected at P7. The apical third of the left ventricle that contained the injured area was then dissociated into single cells. Single cells were then captured into droplets using microfluidics followed by preparation of single cell cDNA libraries as previously described (Kobayashi et al., 2020). Samples were sequenced on a single lane of an Illumina HiSeq. Sequencing reads were mapped to Ensembl release NCBIM37.67. Aligned reads were binned and collapsed onto cellular barcodes using the Drop-seq pipeline v1.13.3 (http://mccarrolllab.com/dropseq), resulting in a raw digital expression matrix containing the count of unique UMIs for each gene in each cell (Macosko et al., 2015).

Expression analysis was performed by the Duke Bioinformatics Shared Resource, using the Seurat package (v2.2.0) in R/Bioconductor (Butler et al., 2018). Cells with gene counts over 4000 or less than 200 were removed. Genes expressed in less than three cells were also removed. To avoid over-filtering CMs, we did not filter cells based on high expression of mitochondrial genes. Gene expression for 1721 cells across 12,754 genes was normalized by scaling by the total number of transcripts, multiplying by 10,000, and log transformation. Unwanted sources of variation were adjusted for by regression on the number of detected unique molecular identifiers (UMIs) using the ‘ScaleData’ function. We used the JackStraw method to determine the number of statistically significant principal components (PCs). The ‘FindClusters’ function was used to identify cellular clusters, using the 20 significant PCs and a 0.6 resolution. We used the t-SNE method to visualize cells across 20 PCs in two dimensions. To identify cell types, we used the ‘FindAllMarkers’ function. Parameters were set to test all genes, expressed in at least 10% of the cells in each cluster, for differential expression.

Anti-VEGFA antibody (B20-4.1.1) was a kind gift from Genentech. Neonatal mice were injected with 5 µg/g of anti-VEGFA or vehicle at P0 and P4. Adequate injection of anti-VEGFA was confirmed by significant weight loss compared with control animals at the time of tissue harvest.

Mice were given 60 µg/g pimonidazole (Hypoxyprobe, HP3-100kit) intraperitoneally 1.5 h before tissue harvest. Tissue sections were stained with an anti-pimonidazole antibody (Hypoxyprobe, HP3-100kit, 1:100) and an antibody against α-actinin. For each heart, nine images of the left ventricle for three separate cardiac sections were obtained using an AxioImager M1 microscope. Images were quantified using a CellProfiler routine to create an image mask based on the α-actinin stain and to determine the mean intensity of anti-pimonidazole stain.

AAV encoding Vegfa-164 under the control of a CMV promoter was generated by PCR amplifying the Vegf164 isoform of Vegfa from a mouse endothelial cDNA library and cloning this fragment into the pTR2-eGFP with replacement of EGFP (Hewitt et al., 2009). AAV, serotype 9, was generated through the Duke Viral Vector Core. AAV virus encoding GFP under the control of a CMV promoter was obtained from the Duke Viral Vector Core. For intracardiac injections, 12 µl of AAV at 1×10⁹ vg/µl was injected circumferentially around the injured area immediately after CI (Fig. S2E,F).

The coronary vasculature was visualized by injection of a low viscosity polyurethane resin (PU4ii, VasQTec) into the apex of the left ventricle (Lee et al., 2019). Whole-mount images were obtained using a Zeiss Stemi 508 stereoscope.

P21 hearts were serially sectioned from base to apex and stained with Acid Fuchsin Orange G, as previously described (Karra et al., 2015). Stained slides were scanned using a Leica Aperio Digital Pathology Slide Scanner. Scar area as a percentage of the left ventricle was scored by a blinded reader. The average scar percentage across the five largest sections was determined for each heart.

Fractional shortening was determined by echocardiography in conscious P21 mice by the Duke Cardiovascular Physiology Core, using a Vevo 2100 high resolution imaging system (Visual Sonics) and the previously described protocol (Jean-Charles et al., 2017). M-mode, long-axis images were analyzed by three blinded readers and an averaged reading is reported.

Two-sample MR analysis was performed using summary statistics from GWAS of circulating protein level measured with Olink proximity-extension assay in 30,931 individuals of European ancestries and GWAS of cardiac MRI-derived left ventricular phenotypes in 16,923 UK Biobank participants (Table S4) (Aung et al., 2019; Folkersen et al., 2020). To represent effect in the original unit of measurements, the standardized estimates from both GWAS were back-transformed by multiplying the MR effect estimate with the estimated standard deviation of the traits. MR instruments for VEGFA were selected from variants that: (1) are available in both GWAS; (2) located within the 200 kb-flanking region from the genomic coordinate of VEGFA gene; (3) have a minor allele frequency>0.01; (4) have a P-value for association with circulating VEGFA level<1×10⁻⁴; and (5) are LD-clumped to an r² threshold of 0.4 using ‘plink’ with --clump option (Purcell et al., 2007). Mendelian randomization analysis was performed using an inverse-variance weighted model accounting for correlation between instruments, implemented in the ‘MendelianRandomization’ R package (Burgess et al., 2016; Yavorska and Burgess, 2017). Correlation between instruments was estimated from a random sample of 10,000 UK Biobank participants of European ancestries (Bycroft et al., 2018).

To test the robustness of the cis-MR analysis, a series of sensitivity analyses were performed using varying parameters for instrument selection and MR models. For instrument selection, LD clumping was performed using five different r² thresholds (0.05, 0.1, 0.2, 0.4 and 0.6) and six different P-value thresholds for association with circulating VEGFA level (no threshold), 1e⁻², 1e⁻³, 1e⁻⁴, 1e⁻⁵ and 5e⁻⁸). Four different MR models were also tested for each instrument set (inverse-variance weighted, MR-Egger, principal component MR with 90% variance explained and principal component MR with 99% variance explained) (Bowden et al., 2015; Burgess et al., 2017). To account for correlation between instruments, we implemented the method proposed by Burgess et al. using an instrument correlation model derived from genotype data of a random subset of 10,000 UK Biobank participants of European ancestry (Burgess et al., 2016).

All data are presented as mean±s.e.m. and proportions as proportion±95% confidence interval. Statistical analysis between two groups was performed using a two-tailed unpaired t-test. Levene's test was used to ensure that variances were not unequal between groups. Pairwise comparisons between multiple groups were performed using a two-tailed unpaired t-test test with correction for multiple testing using Holm's method. A P<0.05 was set as an a priori threshold for significance. Statistical analysis and plots were generated in R using the dplyr, ggpubr and ggplot2 packages (Kassambara, 2017; Wickham, 2011; https://CRAN.R-project.org/package=dplyr). For each experiment, individual data points are presented in the plot and the sample size is specified in the figure legend.

We thank Ken Poss, Nenad Bursac, Doug Marchuk and Howard Rockman for comments and discussion. We thank Dr Helene Fradin Kirshner from the Duke Bioinformatics Shared Resource for assistance with the analysis of scRNA-seq data.

The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.200654.reviewer-comments.pdf