Neurovascular ingrowth after soft tissue trauma accompanies transcriptional changes within vascular cells

In prior work, we observed neuronal sprouting to occur during the development of extremity trauma associated heterotopic bone.6 Using the same model of HO induction, we sought to confirm this finding and relate axonal ingrowth to vascular proliferation and maturation (Fig. 1a–d, Fig. S1). Confirming prior findings,6 innervation of the uninjured tendon was limited and confined to the peritenon as demonstrated by Beta III Tubulin immunofluorescent staining (TUBB3, Fig. 1a, left). Upon HO inducing tenotomy, brisk ingrowth of tendon-associated nerves was found at day 7 (previously characterized fibroproliferative phase of injury,6 Fig. 1b) and progressed at day 21 (previously characterized cartilaginous phase of injury,6 (Fig. 1c). Quantification of TUBB3+ axons confirmed an increase in overall staining intensity, but also an increasing complexity of neuronal projections with increasing junctions and branchpoints overtime (Fig. 1d). CD31 immunofluorescent staining was performed to highlight angiogenesis over the same time course (Fig. 1e–g). A similar spatial pattern of vascular growth was observed, in which endothelium was principally confined to the peritenon among uninjured tendons (Fig. 1e), and a progressive increase in CD31+ vascular channels was found at 7 and 21 days after HO induction (Fig. 1f, g). In parallel to changes in neuronal histomorphometry, quantification of vascularity demonstrated similar kinetics of progressive increase in overall CD31+ endothelium, as well as an increase in the mean diameter of vascular channels, presumed to reflect a more mature vascular network (Fig. 1h).

Fig. 1

Neurovascular ingrowth and transcriptional profile after soft tissue trauma. Traumatic injury (complete Achilles tenotomy with dorsal burn to induce heterotopic ossification) was performed on 6–8 weeks old, male and female C57BL/6 J mice, with analyses at 7 and 21 days. Day 0 indicates uninjured control. a–c Representative pan-neuronal Beta III Tubulin (TUBB3) immunofluorescent staining within sagittal sections of tendon injury site, day 0, 7 and 21 post-tenotomy. d Quantification of TUBB3+ nerves within the injury site over time, including overall TUBB3 staining intensity, number of nerve junctions, and number of nerve branches. e–g Representative immunofluorescent staining of the endothelial marker CD31 over the 21-day time course post-injury. h Quantification of CD31+ blood vessels within the injury site over time, including overall CD31 staining intensity and vessel diameter. i–m Single cell RNA sequencing analysis of pooled cells from the injury site, obtained at day 0, 7 and 21 post-injury. i UMAP visualization of endothelial and pericyte clusters at all timepoints, and (j) UMAP visualization indicating cells at day 0, 7 and 21 post-injury. k Percentage of endothelial cells and pericytes across time as a function of both cell types combined. l, m GO term analysis of endothelial cells (red) and pericytes (blue) at day 7 and day 21 post injury in relation to uninjured control. n Modular index scoring of cell proliferation, vessel growth, and vessel maturation in endothelial cells across time points. Gene lists for modular scores are shown in File S1. Dashed white lines indicate edges of Achilles tendon. Scale bars: 100 μm. For histology, n = 4 animals per timepoint. For scRNA-Seq, n = 3 animals per timepoint. Individual dots in scatterplots represent values from single animal measurements, while mean and one SD are indicated by crosshairs and whiskers. *P < 0.05; **P < 0.01 in relation to D 0 timepoint. Two-tailed student’s t test was used for two-group comparisons (b, d). For more than two groups, one-way ANOVA followed by post hoc analysis was performed (n). Non-parametric data was analyzed using the Kruskal–wallis with Dunn’s multiple comparison test (n)

Single cell RNA Sequencing (scRNA-Seq) profiling of the HO induction site was next examined within this time frame (Fig. 1i–m). Of 22,902 total cells, 4043 cells were identified as endothelial cells (Pecam1high, Cdh5high) or pericytes (Acta2high, Pdgfrbhigh) for further analysis (Fig. 1i, Fig. S2). Both cell clusters were represented across all timepoints (Fig. 1j), with a relative enrichment in pericyte frequency after injury (Fig. 1k). To characterize the vascular response to traumatic injury, GO term analysis of both endothelial cells and pericytes was performed at 7 and 21 days in comparison to uninjured control (Fig. 1l, m). Seven days after injury, both vascular cell types showed enrichment in GO terms related to cell proliferation, angiogenesis, and the regulation of inflammation (Fig. 1l). After 21 days, vascular cell types showed enrichment for GO terms related to smooth muscle cell proliferation, cell adhesion and migration, matrix organization, and TGFβ signaling (Fig. 1m). Having observed these dynamic changes over time, we next examined shifts in endothelial cells using gene modules related to cell proliferation, vessel growth and maturation (Fig. 1n). Confirming prior findings, module scorings of proliferation and vessel growth acutely increased at 1 week post-injury, and after 3 weeks re-approximated the state of uninjured endothelium. In contrast, a module score of ‘vessel maturation’ demonstrated the inverse trend, in which an acute reduction was observed at 1 week post-injury, which re-approximated the score of uninjured endothelium at 3 weeks post-injury. Thus, neurovascular growth and network maturation accompanies trauma-induced HO, and the period between 1 and 3 weeks after injury is a dynamic time frame to examine this process on a transcriptional and histologic level.

Endothelial cell subclustering highlights an increase in type H vessels after HO induction

Having demonstrated shifts in endothelial cell phenotype after injury, endothelial cell subclustering was next performed (Fig. 2). Upon re-clustering, three subclusters of Pecam1high, Cdh5high cells were identified: Endothelial 1, endothelial 2, and a small ‘transitioning’ cell population (Fig. 2a), the latter of which expressed a combination of endothelial and pericyte markers (e.g., Mcam, Pdgfrb, Acta2). All populations were represented before and after HO induction, however endothelial 2 was greatly overrepresented at 7 and 21 days in comparison to baseline uninjured samples (Fig. 2b, c, 58.0% and 35.4% of endothelial cells at 7 and 21 days in comparison to 6.4% at baseline). Characterization of the endothelial 2 cell subcluster revealed a high proliferation phenotype, as determined by module score in relation to other cell clusters (Fig. 2d). Tip cells and stalk cells represent two distinct endothelial cell phenotypes in the nascent sprout of vessels, one guides the growing vessel (tip) and the other that proliferates and follows (stalk).34 Further transcriptomic characterization revealed a stalk-like cell phenotype in the endothelial 2 subcluster, while endothelial 1 demonstrated enrichment for a tip cell-like phenotype, using module scoring (Fig. 2e). GO term analysis of these clusters again demonstrated relative enrichment in terms related to proliferation, angiogenesis and response to extremity injury among the endothelial 2 subcluster, whereas endothelial cluster 1 shows a regulatory phenotype (Fig. 2f) giving further evidence of tip and stalk cell phenotypes.

Fig. 2

Endothelial cell subclustering reveals distinct cell populations after soft tissue trauma. a–h ScRNA-seq analysis of pooled endothelial cells from the HO induction site at 0, 7 and 21 days post-injury. a, b UMAP visualization of all endothelial cells demonstrating (a) endothelial subclusters, and (b) cells across time (day 0, day 7, and day 21). c Percentage of endothelial cell subclusters across time. d Cell proliferation index across endothelial cell subclusters as determined by modular scoring, day 7. e Modular index scoring of stalk cell-like and tip cell-like phenotypes in endothelial subclusters. Gene lists for modular scores are shown in File S1. f GO term analysis of endothelial 1 and endothelial 2 subsets, based on differentially expressed genes between the groups. g Emcn expression among endothelial sub-clusters across time. h Hif1α (Hypoxia-inducible factor 1-alpha) expression among endothelial sub-clusters across time. i Representative immunofluorescent staining of EMCN (Endomucin) and CD146 within the tendon injury site at day 7 and 21 in relation to uninjured control (D0). Dashed white lines indicate edges of Achilles tendon. Scale bars: 100 μm. For histology, n = 4 animals per timepoint. For scRNA-Seq, n = 3 animals per timepoint. *P < 0.05; **P < 0.01 in relation to endothelial subcluster 1. For two groups, non-parametric data was analyzed using the Kolmogorov–Smirnov test (e, g, h). For more than two groups, Non-parametric data was analyzed using the Kruskal–wallis with Dunn’s multiple comparison test (d)

Much attention has been paid to endothelial subtypes involved in osteo-angiogenic coupling, including type H vessels characterized by high expression of endomucin (Emcn).35 To our knowledge, the documentation of type H vessels in HO has not been performed. Indeed, we observed that the endothelial 2 subcluster expressed higher levels of Emcn in comparison to the endothelial 1 subcluster (Fig. 2g). Moreover, the endothelial 2 subcluster demonstrated higher expression of Hif1α (Hypoxia inducible factor 1 subunit alpha) at baseline, and upon HO induction - high Hif1α expression has been previously linked to a type H vessel phenotype36 (Fig. 2h). Results thus far by scRNA-Seq suggested early expansion of a proliferative endothelial cell subset with markers suggestive of type H vessel identity. To further examine this, immunofluorescent staining for Emcn was performed of the injury site in relation to uninjured control, using co-immunohistochemistry for the pan-vascular marker CD146 (Fig. 2i). Results showed infrequent Emcnhigh vessels within uninjured tendon, primarily restricted to the peritenon. A progressive increase in Emcnhigh vessels was observed at 7 and 21 days after HO inducing injury. Together, these data suggest an increase of endothelial cells with a type H phenotype at early timepoints in the evolution of trauma-induced HO. To our knowledge this is the first time that type H vessels have been reported in the context of HO which contributes to our understanding of this pathology and may help better develop new strategies for treatment.

Denervation prevents vascular ingrowth and shifts endothelial phenotype after soft tissue trauma

Our previous studies demonstrated that surgical denervation via the use of sciatic neurectomy largely prevented the formation of HO in our tendon injury model (97.0% reduction in bone formation as determined by microCT imaging analysis).6 However, the extent to which vascularity was altered with surgical denervation (i.e., the dependence of an HO-associated angiogenic response on intact innervation) had not been assessed. Here, ipsilateral sciatic neurectomy (or sham surgery with visualization only of the nerve) was performed at the time of HO induction (Fig. 3). As expected and confirming prior results,6 sciatic neurectomy led to a 90.1% reduction in nerve fiber frequency at the HO site (Fig. 3a, 21 d post-injury, assessed by TUBB3 immunofluorescent staining). Next, CD31 immunofluorescent staining showed a similarly dramatic reduction in vascularity after sciatic neurectomy (Fig. 3b, 92.4% reduction in CD31 immunostaining in comparison to sham control).

Fig. 3

Neurectomy impairs blood vessel infiltration after soft tissue trauma. Animals underwent proximal sciatic neurectomy or sham surgery at the time of HO induction. a TUBB3 (Beta III Tubulin) immunostaining and quantification within the tendon injury site among sham or neurectomized mice, 21 d post-injury. b CD31 immunostaining and quantification within the tendon injury site among sham or neurectomized mice, 21 d post-injury. Dashed white lines indicate edges of Achilles tendon. c–l ScRNA-seq analysis of the HO site at 7 d post injury, among sham or neurectomized animals. c UMAP visualization of subclustering of endothelial and pericyte populations. d Violin plots of characteristic gene markers identifying four vascular subclusters. Pecam1 (Platelet endothelial cell adhesion molecule 1), Cdh5 (Cadherin 5), Pdgfrb (Platelet-derived growth factor receptor beta), Acta2 (Actin alpha 2), Mcam (Melanoma cell adhesion molecule). e Modular index scoring for cell proliferation across subclusters. f Module index scoring for ‘stalk cell-like’ and ‘tip cell-like’ phenotypes across Endothelial 1 and Endothelial 2 subclusters. g Emcn expression across Endothelial 1 and Endothelial 2 subclusters, by violin plot. h GO terms enriched in Endothelial 2 subcluster as compared to Endothelial 1. i UMAP visualization of vascular cell subclusters with or without sciatic neurectomy. j Cellular percentage of Endothelial 2 subcluster with or without sciatic neurectomy. k Emcn expression in Endothelial 2 subcluster in the context of sham or neurectomy conditions, as shown by violin plot. l GO term analysis of endothelial 2 subcluster in the context of sham and neurectomy surgery. m Module index scoring of proliferation, VEGF (Vascular endothelial growth factor) signaling, FGF (Fibroblast growth factor) signaling, and Wnt signaling under neurectomy conditions in comparison to sham control. Values greater than 1 indicates increased signaling activation in neurectomy group while values less than 1 indicates reduced expression among the neurectomy group. Scale bars: 100 μm. For histology, n = 4 animals per condition. For scRNA-Seq, n = 3 animals per condition. Individual dots in scatterplots represent values from single animal measurements, while mean and one SD are indicated by crosshairs and whiskers. *P < 0.05; **P < 0.01 in relation to endothelial subcluster 1 (e–g) or sham control (a, b, k, l). A two-tailed student’s t test was used for two-group comparisons (a, b, f, g, k, l). For more than two groups, non-parametric data was analyzed using the Kruskal–wallis with Dunn’s multiple comparison test (e). Non-parametric data was analyzed using the Kolmogorov–Smirnov test (f, g, m)

To further understand the transcriptomic changes in vascular cells with or without intact innervation, re-analysis of our prior scRNA-Seq data was performed in which HO induction was performed combined with or without sciatic neurectomy, and cells of the HO site were analyzed at 7 days after injury (GSE163446).6 Of a total of 16 823 cells, 3 616 vascular cells were captured and distributed across 4 subclusters with strong similarity to the previously analyzed dataset, including endothelial cell 1, endothelial 2, pericytes, and again a ‘transitioning’ cluster with markers shared between endothelial and pericytes (Fig. 3c, d). Consistent with prior observations, the endothelial 2 subcluster showed the highest proliferative index by module scoring (Fig. 3e), expression more consistent with a stalk cell-like phenotype (Fig. 3f), and high Emcn expression (Fig. 3g). GO term analysis for the endothelial 2 subcluster again showed relative enrichment for biological processes such as endothelial proliferation, migration and angiogenesis (Fig. 3h). Thus, this scRNA-Seq dataset demonstrated overall similar vascular cell populations for analysis, and endothelial subcluster 2 again appeared to represent a proliferative, type H-like endothelium appropriate for focused comparison.

Next, comparisons were performed between sham and neurectomy conditions (Fig. 3i–m). Each vascular subcluster was present under both sham and neurectomy conditions (Fig. 3i), and showed overall a slight increase in frequency across sham as compared to neurectomy samples (Fig. 3j). The endothelial 2 subcluster showed a shift in phenotype in response to neurectomy as observed by a reduction in Emcn expression (Fig. 3k), as well as GO term analysis (Fig. 3l). For example, terms related to angiogenesis, tissue remodeling and MAPK signaling were enriched among control (sham operated) conditions. Focused analysis of specific pathways was further undertaken using gene module scores (Fig. 3m). Proliferative indices were reduced under neurectomy conditions. Signaling pathway changes were also observed including reductions in VEGF, FGF, and WNT signaling pathway scoring among endothelial cells of neurectomized mice (Fig. 3m). Thus, surgical denervation blunts HO-associated angiogenesis accompanied by a shift in endothelial cell transcriptional phenotypes, which may be requisite features of HO formation and provides new insights into molecular pathways that could be targeted.

Pharmaco-genetic or pharmacological inhibition of TrkA signaling impairs angiogenesis at sites of trauma

Surgical denervation using the sciatic nerve represents a mixed motor-sensory deficit which may have multiple parallel effects on angiogenesis and HO formation within the extremity. Next, two established models were assessed to specifically inhibit TrkA-expressing nerves within the context of an HO inducing injury (Fig. 4). Here, TrkA-expressing nerve ingrowth was first inhibited using a chemical-genetic approach, in which a point mutation in TrkA (TrkAF592A) renders mutant mice susceptible to small molecule mediated inhibition of TrkA catalytic activity using the small molecule 1NMPP1 (Fig. 4a). Previously, this animal showed 74.1% reduction in heterotopic bone by microCT imaging in comparison to wildtype control animals.6 In similarity to past observations,6 TrkAF592A mutant animals showed a prominent reduction in HO site innervation (Fig. 4b, c), quantified as a 74.5% reduction (Fig. 4d). Injury site vascularity was next assessed in TrkAF592A mutant animals using CD31 immunofluorescent staining. Results showed a significant impairment in HO site angiogenesis in the context of TrkA inhibition (Fig. 4e, f), quantified as a 65.1% reduction in overall staining intensity (Fig. 4g). Emcn and CD146 co-immunofluorescent staining and quantification was next performed (Fig. 4h–j). Here, a similar reduction in injury site vascularity was observed among TrkA mutant animals, including a notable reduction in the numbers of Emcnhigh blood vessels (84.8% reduction at 21 d post injury), as well as CD146 (85.9% reduction at 21 d post injury). Next, TrkA inhibition was achieved in a parallel method by oral administration of the small molecule inhibitor AR786 to wildtype mice during the formation HO (Fig. 4k). In our past work, AR786 led to a prominent reduction in the amount of cartilage after HO induction (69.3% reduction at 21 d after injury6). As in prior experiments, TrkA inhibition via AR786 treatment led to a significant reduction in injury site innervation (Fig. 4l–n, 88.1% reduction), but also led to a prominent reduction in vascular density (Fig. 4o–q, 76.6% reduction). Thus, three independent methods to inhibit TrkA-expressing sensory nerve ingrowth led to significant deficits in extremity injury associated angiogenesis, further supporting neuron-to-endothelial signaling after musculoskeletal trauma.

Fig. 4

Genetic and pharmacological inhibition of TrkA signaling impairs vascular ingrowth after soft tissue trauma. a Schematic representation of experiment. (Created in Biorender.com). TrkAF592A mice have a point mutation which renders them susceptible to temporally controllable inhibition of TrkA catalytic activity, via administration of the small molecule 1NMPP1. Either TrkAWT or TrkAF592A animals (6–8 weeks old, male and female mice) were subjected to HO induction and both administered 1NMPP1 throughout the study period. Analyses performed at 21 d post-injury. b–d Representative immunohistochemical staining for TUBB3 (Beta III Tubulin, appearing green) and quantification among TrkAWT or TrkAF592A animals, as assessed at the distal tenotomy site as visualized using sagittal cross-sections. e–g Representative immunostaining for CD31 (appearing red) and quantification among TrkAWT or TrkAF592A animals. h–j EMCN (Endomucin) and CD146 co-immunohistochemical staining (appearing red and white) and quantification among TrkAWT or TrkAF592A animals. k Schematic representation of experiment. C57BL/6 J mice (6–8 weeks old, male and female mice) were administered AR786 or vehicle control throughout the study period. HO induction was performed (burn/tenotomy), with analysis at day 21. l–n Representative immunohistochemical staining for TUBB3 and quantification among control or AR786 treated animals, as assessed at the distal tenotomy site. o–q Representative immunostaining for CD31 (appearing red) and quantification among control or AR786 treated animals. Dashed white lines indicate edges of Achilles tendon. Scale bars: 100 μm. n = 4–5 animals per group. Individual dots in scatterplots represent values from single animal measurements, while mean and one SD are indicated by crosshairs and whiskers. * P < 0.05; ** P < 0.01. A two-tailed student’s t test was used for all comparisons

Traumatic HO incites a pro-angiogenic gene profile among sensory neurons

Having shown the dependence of neural ingrowth for angiogenesis at sites of HO predisposing trauma, we next set out to identify potential vasculogenic growth factors enriched within sensory neural populations. First, re-analysis of a scRNA-Seq dataset of mouse dorsal root ganglia (DRG) cells was performed (8–12 weeks old C57BL/6 J mice).37 In agreement with prior reports, four main cell clusters were identified, including: Neurofilament (NF), peptidergic nociceptors (PEP), non-peptidergic (NP) and tyrosine hydroxylase (TH) based on the expression of specific markers (Fig. 5a, b). Ntrk1 (encoding TrkA) was found in peptidergic cells and to a lesser extent in NF cells (Fig. 5b). Next, candidate secreted pro-vasculogenic factors were examined in this dataset. Heatmap demonstration showed unique angiogenic gene profiles across cell clusters (Fig. 5c). For instance, among the Ntrk1low NF neurons Sema5a, Fgf18, Fgf1, Vegfc and Wnt5a demonstrated enrichment. Among the Ntrk1high PEP neurons, Vegfa and Pdgfa among others showed enrichment. The impression that Ntrk1-expressing neurons showed high expression of vasculogenic regulatory genes was confirmed via module score analysis, in which NF and PEP neuron clusters showed high relative angiogenesis related genes in comparison to Ntrk1-non-expressing cell clusters (Fig. 5d).

Fig. 5

TrkA+ peripheral nerves express key angiogenic genes after soft tissue trauma. a–d Re-analysis of scRNA-Seq of mouse lumbar DRG neurons (8–12 weeks old male C57BL/6 J mice, n = 3 mice, 11 191 neuronal DRG cells were examined) (a) UMAP visualization of 4 cell clusters including Neurofilament (NF), Peptidergic nociceptors (PEP), Non-peptidergic (NP) and Tyrosine hydroxylase (TH) in pooled cells. b Violin plots demonstrating characteristic markers of each cell cluster. Nefh (Neurofilament Heavy Chain), Calca (Calcitonin Related Polypeptide Alpha), Tac1 (Tachykinin Precursor 1), Mrgprd (MAS Related GPR Family Member D), Ntrk1 (Neurotrophic Receptor Tyrosine Kinase 1), Th (Tyrosine Hydroxylase). c Heatmap of angiogenic gene expression among each cell cluster. Sema5a (Semaphorin 5 A), Fgf18 (Fibroblast Growth Factor 18), Fgf1 (Fibroblast Growth Factor 1), Vegfc (Vascular Endothelial Growth Factor C), Wnt5a (Wnt Family Member 5 A), Vegfa (Vascular Endothelial Growth Factor A), Pdgfa (Platelet Derived Growth Factor Subunit A), Vegfb (Vascular Endothelial Growth Factor B). d Module score analysis of angiogenic genes among each cell cluster. e–g qPCR analysis of angiogenic gene expression among ipsilateral lumbar DRGs after HO induction, in comparison to uninjured control. Ipsilateral L4-L6 pooled for analysis. e Schematic of qPCR analysis of angiogenic gene expression among ipsilateral lumbar DRGs after HO induction. Created in Biorender.com. f Heatmap of 15 pro-angiogenic genes enriched in Ntrk1-expressing neurons after HO induction, in comparison with uninjured control. Analysis performed at 3 d post-injury. g Representative most up-regulated genes upon HO induction, including Fgf18, Pdgfc, Vegfa, and Vegfd. h Representative immunofluorescent staining of ipsilateral L4 DRG after HO induction in relation to uninjured control, including FGF18 (Fibroblast Growth Factor 18), PDGFC (Platelet Derived Growth Factor Subunit C), VEGFA (Vascular Endothelial Growth Factor A), and VEGFD (Vascular Endothelial Growth Factor D). Analysis performed at 3 day post-injury, and quantification performed in relation to uninjured control. Scale bars: 100 μm. n = 3 or 4 animals per group. Individual dots in scatterplots represent values from single animal measurements, while mean and one SD are indicated by crosshairs and whiskers. * P < 0.05; ** P < 0.01. A two-tailed student’s t test was used (g, h). Non-parametric data was analyzed using the Kolmogorov–Smirnov test (d)

Having established a set of secreted pro-angiogenic factors expressed in situ within mouse sensory neurons, we next examined the extent to which these factors would be over-expressed by innervating sensory nerves after HO induction (Fig. 5e–h). Here, ipsilateral lumbar DRGs were harvested 3 days after HO induction in comparison to uninjured control animals (Fig. 5e). Analysis of 15 pro-angiogenic genes by qRT-PCR showed dynamic modulation of most candidates, with significant change in 7 of 15 genes (Fig. 5f). Those genes with the highest increase in relative expression after HO induction included Fgf18, Pdgfc, Vegfa, and Vegfd (Fig. 5g). To confirm these transcriptional changes, immunofluorescent staining of these factors was pursued among innervating DRGs from the ipsilateral lumbar levels (Fig. 5h). Results confirmed a relative increase in immunostaining for all four pro-angiogenic proteins among HO site innervating DRGs (Fig. 5h).

HO induction results in loss of an anti-angiogenic environment

In contrast to the highly vascular environment of heterotopic bone, tendons are, for the most part, a non-vascularized tissue. Maintaining this state under homeostatic conditions may rely on the release of soluble factors by local cells. To better understand the potential balance between neural-derived angiogenic factors and potential anti-angiogenic factors present within the local tissues, we returned to our scRNA-Seq dataset of the HO induction site, analyzed at 0, 7 and 21 days after injury (Fig. 6). Seven basic cell clusters were identified by characteristic gene markers, including mesenchymal progenitor cells (MPCs), endothelial cells, pericytes, Schwann cells, immune cells, skeletal muscle cells and lymphatic endothelium (Fig. 6a, b). A literature review identified a number of secreted molecules related to angiogenic inhibition, (including Lect1, Ifng, Il18, Il1b, Serpinf1, Pf4, Timp2 and Tnni3),38 which identified overall transcriptional enrichment in MPCs, as visualized by violin plot and UMAP projection (Fig. 6c, d). Other well characterized anti-angiogenic factors showed little or no expression in our dataset, including Plg, and Il12. Having demonstrated that tendon-associated MPCs were most enriched for anti-angiogenic factors in a quiescent state, we next examined changes in factor expression across time after injury (Fig. 6e, f). Further analysis of the MPC cluster showed that overall anti-angiogenic molecule expression was significantly reduced after injury at day 7 and starting to return to basal levels at day 21 (Fig. 6e). Of the 8 genes analyzed, Ifng and Il18 were the most upregulated among MPCs within the uninjured tissues, whereas Serpinf1 and Timp2 have higher expression at both baseline and day 21. Lastly, Lect1 and Tnni3 were upregulated at day 21 which may indicate a specific role in the resolution of healing and tendon re-growth (Fig. 6f). These transcriptional results preliminary suggest a balance of angiogenic and anti-angiogenic factors which collectively may play a permissive role in vascular invasion at sites of trauma—what appears to be a necessary feature for HO formation.

Fig. 6

HO induction prohibits anti-angiogenesis phenotype. Single cell RNA sequencing analysis of pooled cells from the injury site, obtained at 0, 7 and 21 days post-injury. a UMAP visualization of seven cell clusters including mesenchymal progenitor cells (MPCs), pericytes, endothelial cells, pericytes, immune cells, Schwann cells, myoblast and lymphatic cells. b Violin plots demonstrating characteristic markers of each cell cluster. Pdgfra (Platelet-derived growth factor receptor alpha), Pdgfrb (Platelet-derived growth factor receptor beta), Acta2 (Actin alpha 2), Pecam1 (Platelet endothelial cell adhesion molecule 1), Ptprc (Protein tyrosine phosphatase receptor type c), Sox2 (SRY-box transcription factor 2), Sox10 (SRY-box transcription factor 10), Pax7 (Paired box 7), Lyve1 (Lymphatic vessel endothelial hyaluronan receptor 1). c Module score analysis of angiogenic inhibition across seven cell clusters. Gene lists for modular scores are shown in File S1. d UMAP visualization of expression of anti-angiogenic factors across 7 cell clusters. e Module score analysis of angiogenic inhibition across time. f Heatmap of 8 anti-angiogenic genes across time, including Ifng (Interferon gamma), Il18 (Interleukin 18), Serpinf1 (Serpin family F member 1), Timp2 (Timp metallopeptidase inhibitor 2), Il1b (Interleukin 1 beta), Pf4 (Platelet factor 4), Lect1 (Chondromodulin), and Tnni3 (Troponin I3). ** P < 0.01, non-parametric data was analyzed using the Kruskal–wallis with Dunn’s multiple comparison test (c, e)