Mesenchymal VEGFA induces aberrant differentiation in heterotopic ossification

Charles Hwang, Simone Marini,, Amanda K.Huber, David M. Stepien, Michael Sorkin,Shawn Loder, Chase A. Pagani, John Li,Noelle D. Visser, Kaetlin Vasquez, Mohamed A. Garada, Shuli Li, Jiajia Xu , Aaron W. James,Yuji Mishina, Shailesh Agarwal, Jun Li, Ching-Yun Hsu, Paul B. Yu and Benjamin Levi

INTRODUCTION

Pathology from excess ectopic bone formation following musculoskeletal injury, or trauma-induced heterotopic ossification (HO),presents a substantial barrier to recovery in 20% of patients with hip replacements, musculoskeletal trauma, neurologic injury, and blast/burn injuries.1-2Patients with HO experience substantive morbidity including chronic pain, restricted joint function, and open wounds. Effective therapeutics are lacking in contemporary medicine aside from surgical excision.However,these procedures can fail to completely reverse the joint contractures and restricted range of motion.3Furthermore, even after a successful excision procedure, recurrence may occur via similar cellular mechanisms as observed in HO from the incident trauma.4-6

While the etiologies of HO are diverse, it often follows a predictable progression of endochondral ossification7-10after acute inflammatory insult with proliferation of progenitor populations of multipotent mesenchymal stromal cells (MSCs),11-13chondrogenesis,angiogenesis,and ossification7resulting in highly vascular regions surrounding the HO.14These observations parallel programming seen in normal developing bone and fracture healing, suggesting overlapping pathways between developing bone and posttraumatic HO such as the intimate coupling of vascular invasion and osteogenesis.15Previous studies have further demonstrated MSCs differentiating into a cartilaginous scaffold and supportive matrix,16forming a hypoxic environment responsible for the upregulation of HIF1a and VEGFA.17-19During normal bone formation, VEGFA facilitates coinvasion of this anlagen by neovasculature and Osterix(Osx)-expressing osteoblast precursors,16with increased delivery of oxygen to support growing metabolic activity.20Similar migratory behavior to primary ossification centers by osteoblast precursors was observed in response to VEGFA secreted by chondrocytes.21Disruption of these processes lead to impaired endochondral ossification.8,22Notably, mice with conditional deletion of VEGFA in Col2a expressing cells demonstrate aberrations in cartilage development and increased chondrocyte apoptosis.13,23The role of VEGFA as a pro-angiogenic/osteogenic signal24-26fundamental to bone development and homeostasis27highlights a promising candidate for uncoupling acute injury and ectopic bone formation.Though previous studies have demonstrated the role of endothelial progenitors and VEGFA28in bone development and in a nontraumatic HO model,previous work did not identify the role of mesenchymal-derived VEGFA in traumatic HO.

Acute posttraumatic inflammation stimulates local13,29and recruited MSCs to secrete a collection of soluble factors30-31including VEGFA.32However,MSCs have been well demonstrated to be a heterogeneous population13and thus the true lineage and source of these progenitor cells have remained unclear.Various origins have been proposed for the progenitor cells responsible for HO including Tie2/vWF/VeCadherin-positive endothelium,33Prx1/Mx-1/PDGFRα-positive muscle interstitial fibro-adipogenic progenitors,34-36scleraxis (Scx) lineage tendon derived stem cells37and a myriad of other identifying markers.13,38-41Further illustrating the vast heterogeneity of involved cell lineages, additional peripheral progenitor niches such as those marked by odd-skipped related (Osr1, Osr2)genes42and engrailed1 (En1)43may also be involved although not yet implicated with HO.

Given the role of VEGFA during normal bone development and the infiltration of mesenchymal stem cells throughout the formation of HO, we hypothesized that interruption of progenitor-cell-derived VEGFA would lead to a reduction in posttraumatic HO formation. We demonstrate that mesenchymal progenitors are a significant source of VEGFA important for ectopic bone formation and cell-specific inhibition of this signal significantly reduces posttraumatic HO.

RESULTS

Extremity trauma and inflammation induces a pro-angiogenic environment characterized by an increase in total endothelial structures Previous studies have demonstrated the upregulation of vascular signaling prior to pre-HO chondrogenesis.28However, highresolution imaging of the actual vessels have not been obtained in relation to HO after trauma. To do this, we performed multiple modalities to corroborate the increase in vascular density at multiple timepoints across the progression of endochondral ossification, including the acute phase after injury (20 h) and plateau of ossification (9 weeks). First, we used near infrared in vivo interrogation using Angiosense 750EX to image the mouse hindlimb and demonstrate a robust increase in blood flow (Fig.1a).Next,to obtain a high resolution and anatomic assessment of the vascular network implicated with HO formation,we performed microfil injection and nano-CT imaging.We visualized two distinct regions of increased vascular growth around the calcaneus in the distal hindlimb(Fig.1b,red inset)and proximal to the tibial/fibular fuse point in the injured leg. HO formed in close proximity to corresponding vascular beds at 9 weeks (as demonstrated by microCT analyzed at density thresholds above the intravascular polymer) leading us to conclude that there was indeed involvement of vasculature following the peripheral injury with both increased blood flow and angiogenesis at the injury site. To further characterize the cell surface markers of the infiltrating cells at the HO site at more subacute timepoints, we performed flow cytometry of the injured tissue, specifically interrogating two specific cell types: potentially increasing endothelium per microfil imaging and a progenitor niche capable of undergoing transition to connective tissues44to contribute to the HO anlagen. FACS analysis of the injury site highlighted a significant enrichment of CD31-Tie2-CD34+CD133+45-47(Fig.1c)endothelial progenitor cells (EPCs) in burn/tenotomy mice at 1 and 2 weeks (2.70 and 3.36-fold increases compared with uninjured contralateral leg,P=0.032 and 0.002,respectively).There was not similar enrichment of mature endothelium at the 1 or 2 week timepoints(Fig.1d).Given previous literature demonstrating the autocrine nature of VEGFA in vascular homeostasis,48this cellular enrichment of EPCs is consistent with a robust induction of neovasculature and also suggest these progenitor cells as a possible source for VEGFA signaling involved in the angiogenesis associated with HO development consistent with previously described nontrauma HO models.28Thus, we demonstrate using several different modalities that during HO formation, a rich, patent vascular network intimately surrounds the bone,is more prominent than in the uninjured extremity, and likely a facilitative influence in the progression of ectopic of endochondral ossification.

Fig. 1 Focal injury superimposed upon systemic inflammation induces robust local angiogenesis and produces enrichment of endothelial progenitors at the injury site. a Angiosense IVIS demonstrates increased perfusion to injured hindlimb at 20 h post surgery (representative images from n=3 per group). b MICROFIL perfusion nano-CT images and micro-CT of mature mice show networks of nascent vessel outgrowth in regions of subsequent HO formation following burn tenotomy (green inset) at 9 weeks. Vessels in the distal hindlimb exhibit dense vascular infiltration (red inset). c FACS gating schema of injury site 2 weeks after burn tenotomy. CD31- Tie2- CD34+ CD133+endothelial progenitor cells (EPC, blue) identified via R1+R3. CD31+ Tie2+ CD34- mature endothelial cells (red) identified via R2+R4.Remaining live events were collected by surface marker depletion. d Proportion of gated EPC and mature endothelium from all live cells as assessed by flow cytometry was quantified at 1 week and 2 weeks post burn tenotomy

Fig. 2 Vascular progenitors are not the source for VEGFA in HO formation. a Resultant seven clusters identifying cellular subpopulations within HO anlagen at respective timepoints. b Expression profiles for vascular progenitor genes Prom1 (Cd133), CD34, and Vecad/Cdh5. c Expression profile for Vegfa d Schematic for induction protocol of endogenous reporter mice. VeCadherin-creERT2/tdTomatofl/fl identify cells with endothelial programming within outlined regions.Sections taken from distal hindlimb of injury.Right:Sections of injury site from wildtype 3 week burn/tenotomy mice labeled for VEGFA. e Tiled immunofluorescent micrograph of 3 week burn/tenotomy injury site in endogenous VeCadherin reporters.Quantifications for thresholded total positive area(FIJI/ImageJ)were performed across biological triplicate,with 2 hpf per specimen for both HO anlagen site vs. surrounding stroma (n=6/group)

Vascular progenitors are not the contributing source for VEGFA at the burn/tenotomy injury site

Given the enrichment of EPCs demonstrated by flow cytometry,we performed a directed single-cell RNA sequencing experiment in order to characterize transcriptome-level expression of Vegfa and the cellular subpopulations producing these transcripts at baseline preinjury (day 0), and days 3, 7, and 21 postinductive burn/tenotomy. Clustering analysis was performed by blinded biostatistician. In order to characterize the perturbation from baseline caused by inductive burn/tenotomy, clusters were first preliminarily generated unsupervised, and subsequently rank correlated by centroids.This allowed to consolidate clusters into seven cell families, and subsequently compare them across replicates both within a single timepoint and across all timepoints. Clusters show high intra-timepoint stability, as demonstrated by the red diagonals in rank correlation heatmaps for each timepoint and PCA plots.The only notable exception to correlation with self was cluster A that showed a reciprocal minor correlation with cluster D. Based on this high fidelity, replicates within timepoints were pooled for subsequent analysis, with proportional contribution to each cluster from all replicates confirmed by PCA plots (Supplemental Figs. 1-12). The resulting seven clusters (Fig. 2a) were identified via gene ontology terms and KEGG pathway analyses (Tables 1-7) in addition to characteristic gene expression (Supplemental Figs. 13-15):mesenchymal stem cells (Pdgfra), satellite cells with putative stromal progenitors (Pax7), Schwann cells/neural progenitors(Plp1, Cnp, Dhh), myofibroblasts (MFB) encompassing pericytes and smooth muscle cells (Acta2, Mylk), endothelia/vascular progenitor cells (Cdh5, Pecam, Cd34), lymphocytes (Ms4a1), and myeloid cells including macrophages (Ccr2, Mrc1, Fcgr1), granulocyte (Ccr1, S100a9), dendritic cells (Cd209a). Fold change and adjusted P values for genes of interest are included in Table 8.Notably, the endothelial/vascular progenitors cluster (D) demonstrated expression of both Prom1(Cd133) and Cd34 (Fig. 2b, top/middle), corresponding to the EPC subpopulation identified in flow cytometry (Fig. 1c, d). This cluster also demonstrated highexpression of mature endothelium genes (Supplemental Fig. 1)suggesting a mixed phenotype of the EPC and endothelium identified earlier in the flow cytometry. VeCadherin/Cdh5 (Fig. 2b,bottom)also robustly marked this cellular subpopulation.However,this cluster did not exhibit high Vegfa expression levels (Fig. 2c)across any timepoints,suggesting these cells actually were not the likely source for VEGFA signal. To verify these findings, we looked to reconcile data from previous studies that have identified cells of the endothelial lineage to serve as potential HO progenitor cells using a Tie-Cre33marker. These models, however, did not use inducible systems which allow for adult labeling of cells, which is necessary to assess the true contribution of adult endothelial cells to posttraumatic HO. To address this limitation, inducible tdTomato reporter mice expressing creERT driven by the Cdh5 promoter were treated with tamoxifen with washout prior to burn/tenotomy to interrogate the spatial arrangement of early cells of endothelial lineage. Histological analysis demonstrated Cdh5+cells spatially peripheral to the forming HO lesion but a distinct absence of Cdh5+ cells in the inflamed soft tissue deep to the cutaneous/stromal layers at 1 week (Fig. 2d, green). These cells remained a minor population in the progressing HO anlagen even at the 3 week timepoint (Fig. 2d, yellow), in stark contrast to the copious number of VEGFA-expressing cells seen in the same region(Fig. 2d, right). With the observation of Cdh5 enrichment in the periphery of pre-HO anlagen, we corroborated the significant enrichment in the surrounding stroma compared with pre-HOtissue via thresholded quantification of endogenous reporter micrographs (Fig. 2e). Thus, we verified that progenitors from endothelial lineages were not likely candidate precursors to the formation of HO, nor sources for VEGFA within the HO niche. The relatively few numbers of Cdh5+ cells within the HO analgen is consistent with previous studies showing that mesenchymal cells of nonendothelial lineage (i.e. Prx+, Mx1+, or Scx+) are the traumatic HO progenitors.34,36-37,49Having demonstrated that HO lesions are highly vascular but that cells of endothelial origin do not densely populate the interior of the HO anlagen, we next examined whether VEGFA expression spatially correlated with other cellular subpopulations.

Table 1. Cluster A: Top 10 upregulated terms from KEGG and Gene Ontology Domains

Table 2. Cluster B: Top 10 upregulated terms from KEGG and Gene Ontology Domains

Mesenchymal progenitor cells,instead of endothelial lineage cells,comprise the majority of cells within traumatic HO lesions Unsupervised clustering yielded seven unique clusters (Fig. 2c),representing clusters of infiltrating inflammatory cells including macrophages, dendritic cells, and granulocytes, lymphocytes,endothelial cells/vascular progenitors, satellite cells with putative stromal progenitors, and progenitor subpopulations including mesenchymal stem cells, stromal progenitors, neural progenitors,and vascular progenitors. The mesenchymal stem cell population(cluster A) exhibited a remarkably constrained expression of Pdgfra, Osr2, and En1 (Fig. 3a), all markers for a variety of progenitor populations reported in existing literature. Vegfaexpression was indeed found in the Pdgfra expressing population across all timepoints (Fig. 3b). Interestingly, Vegfa expression was also found in macrophage populations after injury (Fig. 3c).Visualized as a timecourse, the mesenchymal stem cell cluster seems to become more heterogenous as a function of time,differentiating into a population with diluted enrichment in Pdgfra, Prrx1, and Vegfa, along with the expected expression of chondrogenic (Sox9) and preosteoblast (Runx2) markers (Fig. 3d,top). This expression profile is starkly contrasted by robust downregulation in chondrogenic and osteoblastic genes in the vascular cluster (Fig.3d,bottom), further corroborating that these cells are not the progenitors for HO formation, nor do they produce VEGFA. Simultaneously referencing the expression profiles (Fig. 3a, b), Vegfa seems to correspond strongly with a subpopulation of the mesenchymal stem cell cluster strongly expressing Pdgfra.Immunofluorescent histology confirmed robust synthesis of VEGFA at the injury site, colocalizing with PDGFRα positive cells (Fig. 3e), with significant increases in both signals compared with contralateral control and a direct correlation between PDGFRα and VEGFA signal (Fig. 3f).

Table 3. Cluster C: Top 10 upregulated terms from KEGG and Gene Ontology Domains

Vegfa is variably expressed across MSC progenitors and their differentiated progeny

In order to characterize the MSC cluster with more granularity,we performed a subanalysis on the interesting MSC cluster(Fig.4a)in light of the suggestive expression profiles pointing todifferentiating/differentiated cells like chondrocytes (Sox9) and preosteoblasts (Runx2) (Fig. 3d). Transcripts profiling tenocytes(Scx+Col1a1+), MSC (Pdgfra+Sox9-Runx2-), chondrocytes (Sox9+Acan+Col1a1+Col3a1+), and preosteoblasts (Runx2+Col1a1+) were used to stratify subgroups within the MSC cluster ‘A.' These cells were observed to peak In number at day 7, with MSCs (day 7: 1 309 MSCs vs. 269 chondrocytes, 842 preosteoblasts, and 324 tenocytes) the dominant subpopulation comprising this MSC cluster and were subsequently analyzed for their level of Vegfa expression as shown in resulting violin plots (Fig. 4c). Vegfa was found in low levels within residual tendon that was attenuated by day 7. Conversely, at day 3, MSCs along with chondrocytes,and preosteoblasts to a lesser degree were found to be enriched in Vegfa. This transiently dissipated for both MSCs and chondroblasts by day 7. However, more terminally differentiating cells, pre-osteoblasts were found to propagate this Vegfa signal even into day 21.With this sustained signaling from MSCs through more differentiated progeny, we investigated a suitable candidate driver that could simultaneously attenuate Vegfa expression across the predominant MSCs, but also the spectrum of their differentiation following musculoskeletal polytrauma. We observed Prrx1 to overlap robustly with not only Pdgfra cells,but also Sox9 and Runx2 subpopulations as well (Fig. 4d). Given confirmation of transcriptional and translational expression of Vegfa from mesenchymal stem cells and their progeny in the injury niche, we hypothesized that conditional knockout in these cells would produce a reduction in posttraumatic HO.Conditional genetic deletion of VEGFA from Prx1 lineage cells decreases posttraumatic extremity HO Having identified an increase in vascularity at a traumatic HO site and characterized a potential source for vascular signal, we next set out to assess the role of Prx1+lineage cell derived VEGFA.Our single-cell RNA sequencing data reaffirm co-expression of Pdgfra and Prrx1 (Fig. 4d), corroborating our previously demonstrated work showing cells of the Prx1+lineage define the HO progenitor cells.34In addition, VEGFA has been shown to be produced by osteoprogenitor cells and to drive bone development.50Recent studies also demonstrate that VEGFA from cells marked by Col2CreER-mediated recombination, but not vascular cell lineage(Flt1) functions as a survival and differentiation factor for chondrocytes.51Thus, to evaluate the role of VEGFA on ectopic cartilage formation, we deleted VEGFA from cells of the Prx1+lineage.Following VEGF deletion from Prx1+MSCs via Prx1-cre/Vegffl/fl,these mutant mice produced significantly less HO in the distal hindlimb (44% relative reduction, P=0.036 after normalization by cortical thickness due to stunted limb phenotype)when induced by burn/tenotomy (Fig. 5a). Given the density of PDGFRα+cells observed with histology, the reduction in ectopic bone formation following selective gene knockout,implicates VEGFA expression and production by MSCs as a major niche in the acute inflammatory response after burn/tenotomy. Additional histological analyses of injury sites of Prx1-cre/Vegffl/fl9 weeks post injury exhibit a profound absence i.e. significant ablation of VEGFA signaling in the deletion mice compared with littermate controls, despite comparable levelsof PDGFRα+cells(Supplemental Fig.2)validating the degree of gene ablation.

Table 4. Cluster D: Top 10 upregulated terms from KEGG and Gene Ontology Domains

Table 5. Cluster E: Top 10 upregulated terms from KEGG and Gene Ontology Domains

Table 6. Cluster F: Top 10 upregulated terms from KEGG and Gene Ontology Domains

Finally,histology was obtained from human HO specimens with diverse morphologies that reflect the full histologic spectrum of human disease (Fig. 5b, n=10 human specimens). VEGFA immunostaining was performed, which paralleled our findings from mouse tissues. In early spindle cell lesions before frank ossification (Fig. 5b1,2), a high density of VEGF immunoreactive cells were present. In other areas with early and inconscipuous amounts of woven bone, noticeably stronger VEGF immunoreactivity was present within bone-lining cells (Fig. 5b3,4). These findings were consistent across specimens and in areas of more conspicuous lamellar bone (Fig. 5b5,6), in which bone-lining cells were the again prominently VEGF immunoreactive. In samples with HO of a thicker,more‘cortical'appearance(Fig.5b7,8)bonelining cells continued to show VEGF immunoreactivity while osteocytes within thicker bone fragments showed little staining.Thus, and across human samples, VEGF highlights similar populations of HO stromal progenitor cells and osteoblasts as in our mouse models, suggesting a potentially similar pathophysiology and disease course.

DISCUSSION

HO is the formation of ectopic bone highly reliant on angiogenesis to progress through endochondral ossification to mature into its mineralized, final form. VEGFA is regulated by many transcription factors including HIF1a which is thought to be a major stimulus incancer, bone development, and bone pathology.19,34,50,52Having identified the central role of HIF1a in traumatic HO, we set out to determine whether VEGFA, a downstream signaling mediator of HIF1a, is involved in HO formation.

Table 7. Cluster G: Top 10 upregulated terms from KEGG and Gene Ontology Domains

Previous studies have shown that VEGFA orchestrates ossification during normal bone development.53Furthermore, elegant studies have demonstrated that VEGFA produced by cells of Osx-Cre lineage regulate blood vessel recruitment and early osteoblast differentiation21and that VEGFA produced by Col2-Cre-marked cells is necessary for chondrocyte survival.51We chose to target VEGFA in cells of the Prx1-Cre lineage as we were the first to identify these cells as the HO progenitor cells,34a decision further supported by our transcriptomic and histologic data correlating Prrx1 and Pdgfra expression. Notably, in our model, chondrocytes provided limited contribution to Vegfa expression levels.However,preosteoblasts were found to be the responsible cell population for Vegfa expression at later timepoints(D21),demonstrating that in addition to recruited MSCs, their differentiating progeny propagated the robust expression Vegfa. While MSCs were the dominant source for Vegfa expression both in magnitude and volume (Fig. 4), the heterogeneity of Vegfa also deriving from chondroblasts and preosteoblasts further reinforced our strategy in producing VEGFA knockout via a Prx1 driver.Of particular note from a bioinformatics perspective, while the overwhelming majority of clusters demonstrated high specificity in alignment,two clusters A and D showed moderate levels of correspondence to each other following the initial pseudo clustering (Fig. 2a).These clusters correspond to related progenitor cells, bothexpressing a common ancestry given their prrx1 enrichment.However, the musculoskeletal mesenchymal progenitor population (D) exhibited sparse VEGFA signaling, essentially rendering this an inert subpopulation. Here we show that genetic loss of VEGFA from the mesenchymal stem cell population (along with some chondrocyte and preosteoblasts)is sufficient to significantly reduce HO volume.VEGFA inhibition in cells of the Prx-Cre lineage did not, however, completely mitigate all HO formation likely due to alternative sources for VEGFA as demonstrated in our single-cell RNA sequencing, namely inflammatory cells.51,54-55Despite previous literature implicating vascular sources of VEGFA secreted in a cell autonomous manner,48our sequencing data concretely demonstrated a notable absence of Vegfa transcripts.In addition, these cells did not express common markers of chondrogenic and osteoblastic differentiation, cementing their supportive role in the periphery of the developing ectopic bone.

Table 8. Cluster identities-characteristic genes

Histological analysis of our acute injury model demonstrated a lack of spatial proximity between Cdh5-creERT2 cells and VEGFA within the HO region when activated just prior to injury. Cells positive for Cdh5, however, were present in a large amount surrounding the mesenchymal condensation that underwent endochondral ossification corroborating our sequencing results and the Cdh5+ cell's important supportive role to the HO niche.This presence of Cdh5+ cells also corresponded with an increase in blood flow as well as blood vessels likely derived from the identified EPCs, which is consistent with the rubor and swelling seen at a joint predisposed to HO prior to radiographic evidence.

Two pathways could explain our observed phenotype of HO inhibition following loss of mesenchymal-derived VEGFA:decreased angiogenesis secondary to loss of the trophic signal leading to an indirect deleterious effect on endochondral ossification or direct phenotypic changes to the mesenchymal cells themselves. Interestingly, our mesenchymal stem cell population did not exhibit upregulation of any VEGF receptors,including Flt1/VEGFR1, Kdr/VEGFR2, or Flt4/VEGFR3 (data not shown) rendering them unresponsive to VEGFA. This is largely concordant with previous literature characterizing VEGFA activity in fetal tissues.23Given the large influx of inflammatory cells including VEGFA expressing macrophages, one could presume that there should be enough residual VEGFA at the locus of injury to rescue the loss of mesenchymal-derived VEGFA,especially with the former pathway.Given our observed reduction in HO,it seems the latter pathway is more likely. Indeed, previous work by Zelzer et al.23demonstrated that VEGFA loss from mesenchymal stem cells leads to aberrant growth plate development despite no changes to surrounding vasculature,instead pointing to transcriptional changes within the differentiated mesenchymal cells and chondroblasts themselves. Furthermore, in this paradigm, the partial bone phenotype is unsurprising as Pdgfra/Prrx1 cells highly expressing VEGFA were a subpopulation of a larger pool of progenitors. Nevertheless, the targeted ablation of VEGFA was observed to produce marked, clinically-relevant reductions of posttraumatic HO.As such,future work could involve utilization of selective, deliverable next-generational inhibitors such as siRNAs for cell-specific knockdowns to attenuate HO formation.

This work is also the first to our knowledge to further characterize the role for novel progenitor-cell populations rapidly expanding in the literature. Subpopulations characterized as dermal lineage progenitors like engrailed1 (En1)56and fibroadipogenic progenitors like odd-skipped related-1/2 (Osr1/2)57were found in the larger Pdgfra cluster,implicating En1 and Osr1/2 expressing cells as subsets of a larger progenitor population involved with HO formation. These cells appear transcriptionally distinct from Vegfa expressing progenitors, suggesting next steps that will be necessary to further characterize the differentiation capacity of these cells in the HO niche and to what degree they differ in the disease course. In addition, given previous work demonstrating the spatial variation of different lineages i.e.concentration of Prx-positive cells in the distal hindlimb, these separate cells with progenitor phenotypes may perhaps mark more proximal contributing cells to HO formation. Downstream high-throughput analyses facilitated by Next-Gen sequencing modalities as performed in this work will offer additional insights into important phenotypes and candidate targets.

This study has several notable limitations. Our RNA-sequencing results did not offer substantial resolution between mature endothelium (Supplemental Fig. 1) and EPCs. However, the flow cytometry did mark this separation well. Furthermore, this data would suggest that analysis of VeCadherin endogenous reporters would be likely indicating both mature endothelium in addition to the EPCs on histology. Thus, a complete absence of signal necessarily still validates an absence of EPCs.We chose to use Prx-Cre which will mark all cells from embryonic and postnatal development rather than an inducible Cre given our previous studies that demonstrate their central role in HO using this model.34In addition, we do not delineate whether VEGF stimulates MSC differentiation58or just promotes chondrocyte survival. In addition, because of our approach to characterize perturbations of cellular populations found at baseline, these constraints limited are ability to characterize the cell populations with further granularity i.e. myeloid cells vs. individual clusters for macrophages, granulocytes, and dendritic cells. This leaves us potentially susceptible to artifacts in results specifically due to dilution of clusters(as seen in the Vegfa expression patterns in the mesenchymal stem cells.

Notwithstanding these limitations, our findings demonstrate that mesenchymal-derived VEGFA is a crucially expressed signal in the genesis of HO,suggesting that therapies targeting VEGFA may present a reasonable option for reducing this debilitating disease and expanding the repertoire for its prevention. Subsequent studies are required to determine the efficacy and safety profile of targeted siRNA delivery systems and existing systemic anti-VEGFA antibodies such as bevacizumab in more directed wound healing models for HO patients.

Fig. 3 Mesenchymal progenitor cells comprise the majority of VEGFA expressing cells within HO lesions. Clusters from 10× RNA sequencing across days 0,3,7,and 21 timepoints.Expression profiles corresponding to a mesenchymal progenitor cells:Pdgfra,Osr2,and En1；b Vegfa；c and macrophages: Mrc1 and Adgre1. d Expression heatmap delineating timecourse for log-fold changes in expression for genes of interest within(Top)mesenchymal stem cell cluster A and(Bottom)endothelial cells/vascular progenitor cells cluster E.Red and blue signify upregulation and downregulation respectively. e Representative immunolabeling from burn/tenotomy on wild-type mice 1 week after injury. n=3/group. Color dotted squares designate magnified regions.White dotted lines circumscribe residual bellies of Achilles tendon.Dash dot lines outline reference native bone(calcaneus).Scalebars represent 100 μm.f Quantifications for relative signal in 1-week injury site(left leg)vs.uninjured contralateral limb (right leg) following automated thresholding (Yen). Linear regression of PDGFRα vs. VEGFA demonstrates significant correlation

Fig.4 Vegfa is variably expressed amongst MSC progenitors and their differentiated progeny.a Schematic of cluster subanalysis of MSC cells from single-cell RNA sequencing.b Discrete counts of cell types and classification criteria.c Expression levels of Vegfa across all indicated cell types. d Expression profiles of characteristic markers Prrx1； mesenchymal progenitors: Pdgfra； chondrocytes: Sox9； and preosteoblasts: Runx2

MATERIALS AND METHODS

Ethics statement All animal experiments described were approved by the University Committee on Use and Care of Animals at the University of Michigan, Ann Arbor (Protocols: #05909 and 07930). This study was carried out in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health.59All animal procedures were carried out in accordance with the guidelines provided in the Guide for the Use and Care of Laboratory Animals: Eighth Edition from the Institute for Laboratory Animal Research.60

Animals

Mice evaluated for ectopic bone were wild-type C57BL/6J(Jackson Laboratories), Cdh5-Cre/tdTomatofl/fl61, VeCadherin(Cdh5)-creERT2/tdTomatofl/fl62(a gift from Dr Luisa Iruela Arispe,UCLA),Prx1-cre34/Vegfafl/fl, or Prx1-cre/Vegfafl/fllittermate controls. All breeding was performed at the University of Michigan in facilities managed by the Unit for Laboratory Animal Medicine. CreERT reporters were induced with two doses of 1 mg tamoxifen in mice one week prior to burn/tenotomy for adequate washout by the time of injury.

Fig.5 VEGFA knockout and molecular inhibition attenuates formation of HO in the distal hindlimb.Following 9 weeks burn tenotomy,a Prx-Cre；Vegffl/fl； treated WT mice exhibit lower volumes of heterotopic ossification than their corresponding controls at 800 HU. Prx-Cre HO volumes normalized by mean cortical thickness(mm)due to congenital phenotype.n=4-9/group.Bars represent means and SD.*P ＜0.05 vs.Control by Student's t test.b Magnified images of human HO samples immunolabeled for VEGFA:(1,2)pre-HO stroma,(3,4)early HO formation with surrounding stromal fibroblasts,(5,6)intermediate HO,and(7,8)mature HO with terminal ossification.VEGFA indicated by red,scalebar indicates 100 μm

Injury models

All mice received presurgical analgesia consisting of 0.06 mg·kg-1buprenorphine for 48 h, followed by anesthesia with inhaled isoflurane, and close postoperative monitoring with analgesic administration.

Mice received 30% total body surface area partial-thickness burn on a shaved dorsum followed by transection of the left Achilles tendon. Dorsal burn was induced using a metal block heated to 60°C in a water bath and applied to the dorsum for 18 s continuously. HO anlagen was observed by week 3 with mature bone formation visible by microCT by 9 weeks.

In vivo imaging system

Mice were anesthetized via inhaled isoflurane and depilated with Nair. Specimens were imaged before injection for baseline fluorescence levels at 770-nm wavelength using an IVIS Spectrum system(PerkinElmer 124262).Per manufacturer's instructions,100 μL of 20 μmol·L-1AngioSense 750EX (PerkinElmer NEV10011EX)solution was injected into the tail vein followed by re-imaging 24 h post injection at excitation/emission of 750/800 nm. Angio-Sense protocol was performed 20 h post burn/tenotomy (n=3/group).

10× Single-cell RNA sequencing

Baseline uninjured tendon (day 0) and post surgery day 3, 7, and 21 harvested tissue samples were digested for 45 min in 0.3%Type 1 Collagenase and 0.4% Dispase II (Gibco) in Roswell Park Memorial Institute (RPMI) medium at 37°C under constant agitation at 120 r·min-1. Digestions were subsequently quenched with 10% FBS RPMI and filtered through 40 μm sterile strainers.Cells were then washed in PBS with 0.04% BSA, counted and resuspended at a concentration of ～1 000 cells per μL. Cell viability was assessed with Trypan blue exclusion on a Countess II(Thermo Fisher Scientific) automated counter and only samples with ＞85% viability were processed for further sequencing.

Single-cell 3′ library generation was performed on the 10×Genomics Chromium Controller following the manufacturer's protocol for the v2 reagent kit (10× Genomics, Pleasanton, CA,USA). Following generation of single-cell gel bead-in-emulsions(GEMs), reverse transcription was performed and the resulting Post GEM-RT product was cleaned up using DynaBeads MyOne Silane beads (ThermoFisher Scientific, Waltham, MA, USA). The cDNA was amplified, SPRIselect (Beckman Coulter, Brea, CA, USA)cleaned and quantified then enzymatically fragmented and size selected using SPRIselect beads to optimize the cDNA amplicon size prior to library construction. An additional round of doublesided SPRI bead cleanup is performed after end repair and Atailing.Another single-sided cleanup is done after adapter ligation.Indexes were added during PCR amplification and a final doublesided SPRI cleanup was performed. Libraries were quantified by Kapa qPCR for Illumina Adapters(Roche)and size was determined by Agilent tapestation D1000 tapes. Read 1 primer sequence are added to the molecules during GEM incubation. P5, P7 and sample index and read 2 primer sequence are added during library construction via end repair, A-tailing, adapter ligation and PCR.Libraries were generated with unique sample indices for each sample. Libraries were sequenced on a HiSeq 4000, (Illumina, San Diego, CA, USA) using a HiSeq 4000 PE Cluster Kit (PN PE-410-1001) with HiSeq 4000 SBS Kit (100 cycles, PN FC-410-1002)reagents, loaded at 200 pmol·L-1following Illumina's denaturing and dilution recommendations. The run configuration was 26×8×98 cycles for Read 1, Index, and Read 2, respectively. Cell Ranger Single Cell Software Suite 1.3 was used to perform sample de-multiplexing, barcode processing, and single-cell gene counting (Alignment, Barcoding and UMI Count) at the University of Michigan Biomedical Core Facilities DNA Sequencing Core.

A total of ～200-500 million reads were generated from the 10×Genomics sequencing analysis for each replicate. The sequencing data were first preprocessed using the 10× Genomics software Cell Ranger (10× Genomics Inc., Pleasanton, CA, USA) and aligned to mm10 genome (deposited to the Gene Expression Omnibus database: GSE126060). Downstream analysis steps were performed using Seurat. Cells with fewer than 500 genes, more than 10% of reads mapping to the mitochondrial genome, or more than 60 000 UMIs；and genes expressed in fewer than 10 cells per replicate,were filtered for quality control. The downstream analysis steps for each sample type include normalization,scaling,dimensionality reduction(PCA,t-SNE,and UMAP),unsupervised clustering,and the discovery of differentially expressed cell-type specific markers.We verified the absence of the within-time-point batch effect by visual inspection of the per-sample contribution in the per-day merged set, once we pooled together samples from each timepoint (Fig. 2b). In order to characterize the induced changes in cellular populations following burn/tenotomy injury, provisional clusters were identified unsupervised clustering followed by centroid rank correlation comparison within and between timepoints (days 0, 3, 7, and 21). Provisional clusters that failed to align across these comparisons were discarded,leaving behind cells and clusters that were identified in the steady state(day 0)timepoint and subsequently perturbed across days 3,7,and 21. Final clusters were refined by pooling cells from all the timepoints to create a merged set；correcting for batch effect of the merged set； and performing a supplemental unsupervised clustering. Expression profiles of known markers, GO term and KEGG pathway analysis, along with characteristic expression profiles were utilized to identify the cell types of each cluster. FDR adjustments were performed to determine statistical significance of gene expression fold changes within a cluster vs. all other clusters within the merged set.Enrichment was performed with LRPath.The whole Bioinformatics analysis summarized here is detailed in the supplementary material (Figs. S3-15).

Histology and immunofluorescence

Histologic evaluation was performed at indicated timepoints in wild-type C57BL/6J and tdTomato endogenous reporter mutants following burn/tenotomy. TdTomato and wild-type hindlimbs were fixed at 4°C in 4% paraformaldehyde or 10% neutralbuffered formalin. Specimens were subsequently decalcified in 14% (m/v) EDTA solution for 4-6 weeks at 4°C until deformable manually. Hindlimbs were cryo- or paraffin-embedded and 8 μm-10 μm sections were cut and mounted on Superfrost plus slides (Fisher) and stored at -20°C.

Paraffin sections were consistently selected from similar regions defined by talus as the anatomical landmark, rehydrated,processed for antigen retrieval in citrate buffer for 20 min at 95°C-100°C and blocked with 2% serum, 1% BSA, and 0.1% fish skin gelatin. Human sections were separately treated with trypsin enzymatic antigen retrieval solution (Abcam) for 10 min at 37°C and blocked with 3%hydrogen peroxide(20 min,RT)and 5%goat serum(1 h,RT).Immunostaining of HO anlagen was performed on paraffin sections with the following antibodies:goat anti-mPDGFRalpha (R&D Systems, AF1062) and rabbit anti-mouse VEGFA (EMD Millipore, ABS82). Secondary antibodies consisted of donkey antirabbit or anti-goat Alexafluor-488 or -594 (Invitrogen). Primary and fluorophore conjugated secondary antibodies were diluted 1:50-100 (2-4 μg·mL-1) and 1:200, respectively. Nuclear counterstain was performed with Hoechst 33342 (Life Technologies).Immunostaining of human HO was performed on paraffin sections with polyclonal rabbit anti-mouse VEGFA (Santa Cruz, Cat No. sc-152). Slides were subsequently blocked with BLOXALL solution(Vector)and 2.5%normal horse serum,and incubated with second primary antibody,IMPRESS-AP reagent(Vector)with development and counterstain via Alkaline Phosphatase Substrate and hematoxylin (Vector). Appropriate primary antibody negative controls were run simultaneously with each tested sample.

Microscopy

Cryosections of tdTomato endogenous reporter mice were imaged with a Zeiss LSM 510-META Laser Scanning Confocal Microscope equipped with a helium neon 1 laser (AF594, RFP).Cryosections for immunolabeled samples were performed on a Leica TCS SP8 Laser Scanning Confocal Microscope with tunable laser whose excitation and emission parameters were set to fluorophore manufacturer's instructions.Antibody labeled paraffin sections were imaged with epifluorescent upright scope(Olympus BX-51) equipped with DAPI and dual-cube 488 nm/TRITC filters attached to an Olympus DP-70 high-resolution digital camera.Each site was imaged in all channels and were overlaid in Olympus DPViewer or Photoshop before examination and quantifications in ImageJ. Images were adjusted linearly only in brightness and contrast identically across comparison groups for clarity where necessary.

For quantification of positively fluorescence, images were automatically thresholded via Yen method, and then measured across 1-3 separate hpf per biologic specimen (ImageJ, n=1-3/specimen, N=3 specimens/treatment group). Images were then analyzed for area via the standard ImageJ measurement panel.

Flow cytometry

FACS analysis was performed on the injury site after harvesting the soft tissue of the posterior compartment between the muscular attachment of the Achilles tendon and the calcaneal insertion. Correlating soft tissue was also harvested from the contralateral uninjured leg.Tissue was digested for 60 min in 0.3%Type I Collagenase and 0.4% Dispase II (Gibco) in RPMI at 37°C under constant agitation. Digestions were quenched with RPMI supplemented with 10% FBS and filtered using a 70-μm sterile strainer, transferred to 15 mL conical tubes, and centrifuged at 300×g for 5 min. Tubes were decanted and washed in PBS three times. Digested specimens were blocked for Fc receptors with anti-mouse CD16/32 for 15 min at 4°C, and subsequently stained with the following antibodies: Tie2-PE (12-5987-83, eBioscience),CD31-APC (551262, BD Pharmingen), CD133-AF488 (11-1331-82,eBioscience), and CD34-eFluor450 (48-0341-82, eBioscience).Following 30 min of incubation at 4°C, samples were washed,filtered through a 35 μm mesh strainer,incubated with propidium iodide,and run on a FACSAria II(BD Biosciences)Cell Sorter at the University of Michigan Flow Cytometry Core in the Biomedical Science Research Center. For sorting, after appropriate FSC/SSC,singlets, and viable cells were gated, the CD31 and Tie2 doublenegative subpopulation was enriched for CD34 positive EPCs.Tie2 and CD31 double positive subpopulations were further refined for CD133 and CD34 double-negative mature endothelium. Event counts were normalized to total viable cells to control for possible differences in the amount of tissue harvested. Data were then analyzed using FlowJo software.

MICROFIL imaging

In vivo collateral vessel formation was assessed using MICROFIL.At 9 weeks after burn/tenotomy, wild-type mice were euthanized with cannulation of descending thoracic aorta and venting of inferior vena cava. Mice were perfused with MICROFIL (MV120-blue； Flow Tech Inc.) which was allowed to polymerize overnight at 4°C and subsequently imaged with μCT (n=2/group).

MicroCT analysis

Mouse hind limbs were imaged with μCT at 9 weeks post injury(Bruker, 35 μm resolution and 357 μA/70 kV beam). Images were reconstructed and ectopic bone volume formation calculated with Microview (GE Healthcare v2.2, Parallax rc16). All scans were analyzed with threshold Hounsfield units (HU) of 800, 1 250, or 1 800 to determine the gross volume of mineralized tissues.Additional analyses delineated ectopic bone proximal to the tibial-fibular confluence from bone distal to this landmark and were defined as proximal and distal HO respectively.

Statistical analysis

Primary outcome of interest for a priori power analysis is volume of mature ectopic bone formation. To detect 50% decrease from 7.5 mm3in untreated mice at significance level of 0.05 and power of 0.80, three mice were required per group. Means and SD were calculated from numerical data, as presented in text, figures, and figure legends. Bar graphs represent means with error bars specifying one SD. Sample groups were confirmed to be homoscedastic via Levene's test. Pairwise comparisons were conducted with Student's t test after confirmation for homoscedasticity via F-test. P values are included in figure legends.

ACKNOWLEDGEMENTS

We would like to thank Christopher Breuler, Caitlin Priest, Serra Ucer, and the University of Michigan Orthopaedic Research Laboratories for their technical assistance.We also thank the Department of Radiology at The University of Michigan for the use of The Center for Molecular Imaging and The Preclinical Imaging &Computational Analysis Shared Resource which are supported in part by Comprehensive Cancer Center NIH grant P30 CA046592. We also thank to the University of Michigan Biomedical Research Core Facilities and the DNA sequencing core.

FUNDING

B.L.: Supported by funding from NIH/National Institute of Arthritis and Musculoskeletal and Skin Diseases NIH1R01AR071379 and American College of Surgeons Clowes Award.D.M.S.:Supported by Plastic Surgery Foundation Resident Research Award.M.S.: Supported by Plastic Surgery Foundation National Endowment Award. C.H.:Supported by Howard Hughes Medical Institute Medical Research Fellowship. J.L.:Supported by Vascular Surgery T32 5-T32-HL-076123-14. A.W.J.: Supported by the NIH/NIAMS (R01 AR070773, K08 AR068316, S10OD016374), the Orthopedic Research and Education Foundation with funding provided by the Maryland Stem Cell Research Foundation, and the Musculoskeletal Transplant Foundation. P.B.Y.:Supported by funding from NIH/NIAMS R01 AR057374 and NHLBI R01 HL131910.Y.M.: Supported by funding from NIH/NIDCR R01 DE020843 and DE027662

AUTHOR CONTRIBUTIONS

C.H.,S.M.,A.K.H.,D.M.S.,M.S.,S.Loder,C.A.P.,John Li,N.D.V.,K.V.,M.A.G.,S.Li.,C.Y.H.,J.X., A.W.J., S.A., Jun Li, and B.L. conceived and/or designed the work that led to the submission, acquired data, and/or played an important role in interpreting the results.C.H.,S.M.,M.S.,P.B.Y.,A.W.J.,Jun Li,Y.M.,S.A.,and B.L.drafted or revised paper.All listed authors approved the final version.

ADDITIONAL INFORMATION

The online version of this article (https://doi.org/10.1038/s41413-019-0075-6)contains supplementary material, which is available to authorized users.

Competing interests:A.W.J. serves on the scientific advisory board for Novadip LLC and also receives laboratory financial support from MTF Biologics for research unrelated to the current project. All other authors declare no competing interests..