Identification of distinct non-myogenic skeletal-muscle-resident mesenchymal cell populations.

Mesenchymal progenitors of the lateral plate mesoderm give rise to various cell fates within limbs, including a heterogeneous group of muscle-resident mesenchymal cells. Often described as fibro-adipogenic progenitors, these cells are key players in muscle development, disease, and regeneration. To further define this cell population(s), we perform lineage/reporter analysis, flow cytometry, single-cell RNA sequencing, immunofluorescent staining, and differentiation assays on normal and injured murine muscles. Here we identify six distinct Pdgfra+ non-myogenic muscle-resident mesenchymal cell populations that fit within a bipartite differentiation trajectory from a common progenitor. One branch of the trajectory gives rise to two populations of immune-responsive mesenchymal cells with strong adipogenic potential and the capability to respond to acute and chronic muscle injury, whereas the alternative branch contains two cell populations with limited adipogenic capacity and inherent mineralizing capabilities; one of the populations displays a unique neuromuscular junction association and an ability to respond to nerve injury.

INTRODUCTION

Limb muscles are derived from early myogenic progenitors of the somatic mesoderm. These progenitors initially delaminate from the hypaxial dermomyotome, proliferate, and migrate into the limb field to form dorsal and ventral masses as a response to patterning cues (Buckingham et al., 2003; Deries and Thorsteinsdottir, 2016). The masses of myogenic progenitors undergo further patterning, followed by their differentiation into myoblasts that ultimately fuse into myofibers forming all limb skeletal muscles (Bismuth and Relaix, 2010; Braun and Gautel, 2011). Concomitant with this process, cells of the lateral plate mesoderm (LPM) proliferate and migrate along with the early myogenic progenitors to generate the developing limb field. During limb growth, the LPM-derived mesenchymal progenitors continue to proliferate and become patterned along the proximal-distal, dorsal-ventral, and anterior-posterior axes. They ultimately undergo a series of cell fate choices as they condense and differentiate into chondrocytes (cartilage cells), osteoblasts (bone cells), ligament fibroblasts, and tenocytes (tendon cells) of the limb skeleton, as well as into fibroblastic cells within skeletal muscles that contribute to the formation of muscle-associated fascia and connective tissues (Nassari et al., 2017). Expression and genetic reporter studies have identified a number of transcription factors (Prrx1, Osr1, Tcf7l2/Tcf4) expressed within LPM-derived mesenchymal cells that give rise to some or all of the limb skeletal and connective tissue lineages (Kardon et al., 2003; Logan et al., 2002; Nassari et al., 2017; Stricker et al., 2006). Developmental genetic approaches to eliminate some of these cells, or remove specific signaling (Shh and Wnt) pathway components within these cells, have underscored their requirement during early limb muscle patterning and proper myogenesis (Helmbacher and Stricker, 2020; Hu et al., 2012; Kardon et al., 2002, 2003; Mathew et al., 2011; Stricker et al., 2012; Vallecillo-Garcia et al., 2017).

In postnatal and adult muscles, LPM-derived mesenchymal cells localize to the interstitial regions between muscle fibers and have been further characterized via the utilization of genetic and flow-cytometric approaches exploiting their expression of Platelet-derived growth factor receptor α (PDGFRα), Stem cells antigen 1 (SCA1; Ly6a), CD34, and Hypermethylated in cancer 1 (HIC1) (Joe et al., 2010; Kardon et al., 2003; Mathew et al., 2011; Scott et al., 2019; Uezumi et al., 2010; Uezumi et al., 2014; Uezumi et al., 2011; Vallecillo-Garcia et al., 2017; Stumm et al., 2018), with PDGFRα and/or SCA1 being the most prominent and well-recognized markers for the identification and isolation of these cells. While CD34 is expressed in muscle-resident mesenchymal cells, it is also expressed in muscle stem cells (satellite cells) and therefore makes it a less useful tool in trying to understand their function(s) (Alfaro et al., 2011; Beauchamp et al., 2000). In vitro and in vivo studies indicate that some or all of these muscle-resident mesenchymal cells are multipotent and capable of differentiating into fibroblasts, adipocytes, osteoblasts, and chondrocytes (Joe et al., 2010; Oishi et al., 2013; Uezumi et al., 2010, 2011). Their propensity for either fibrogenic and/or adipogenic differentiation in vitro and in vivo earned them the name fibro-adipogenic progenitors (FAPs). While these fates may be their most common differentiated state, some FAPs have also demonstrated a potential to contribute directly to the myogenic lineage. Specifically, a Twist2/PDGFRα/PDGFRβ-expressing cell population appears to be capable of directly contributing to the myogenic lineage and forming progressively more muscle fibers with age (Liu et al., 2017). Other work has also suggested that the myogenic potential may be regulated by epigenetic reprogramming (Biferali et al., 2021). More recently, a PDGFRα-expressing group of FAPs has also demonstrated an ability to form muscle fibers localized specifically to the myotendinous junction (Esteves de Lima et al., 2021; Yaseen et al., 2021). Collectively, these studies highlight the potential heterogeneity and/or multifunctionality of muscle-resident mesenchymal cells or FAPs during postnatal skeletal muscle development and homeostasis.

In the context of skeletal muscle injury, muscle-resident mesenchymal cells or FAPs establish a defined extracellular matrix, secrete numerous growth factors and cytokines, aid in the clearance of tissue debris (Chapman et al., 2017; Lukjanenko et al., 2019; Schuler et al., 2021; Zou et al., 2008), and communicate with other cell types such as satellite cells and immune cells to facilitate muscle regeneration (Heredia et al., 2013; Kang et al., 2018; Kastenschmidt et al., 2021). Their genetic removal from skeletal muscles or alterations to their basic cellular functions following muscle injury leads to impaired muscle regeneration (Murphy et al., 2011; Uezumi et al., 2014; Wosczyna et al., 2019). Experimentally, these cells generate the intramuscular fat following glycerol injuries and alter normal muscle regeneration (Kopinke et al., 2017). Similarly, FAPs have been implicated in the fatty atrophy of the supraspinatus and infraspinatus muscles following massive rotator cuff tears, at least in murine models (Liu et al., 2016).

FAPs are also significant contributors to the disease pathologies of several skeletal muscle disorders. Dysregulated differentiation and function of FAPs leads to the fibrosis and fatty deposition associated with various muscular dystrophy pathologies (Contreras et al., 2016; Giuliani et al., 2021; Hogarth et al., 2019; Kopinke et al., 2017; Malecova et al., 2018). Further highlighting their multipotentiality in skeletal muscle diseases, FAPs have been identified as a pathologic cell type in heterotopic ossification (HO). In HO, these cells undergo aberrant differentiation into an osteogenic-like cell resulting in mineral deposition and/or bone formation within skeletal muscles, often following muscle injury or damage (Eisner et al., 2020; Leblanc et al., 2011; Oishi et al., 2013; Wosczyna et al., 2012). Specifically, in fibrodysplasia ossificans progressiva (FOP), a genetic form of HO, PDGFRα-expressing FAPs are the pathogenic cell type when an activating mutation in the bone morphogenetic protein (BMP) type I receptor, ACVR1 (ALK2) occurs and leads to robust mineralization of skeletal muscle following muscle damage (Lees-Shepard et al., 2018). The diverse abilities of muscle-resident mesenchymal cells or FAPs to potentially contribute directly to the muscle lineage, generate fatty, fibrotic, and/or osteogenic matrices within diseased or damaged skeletal muscles, and generally be required for normal skeletal muscle patterning, development, and regeneration raises the question as to whether these cells are a single class of multipotent progenitors or whether they are transcriptionally, physically, and/or functionally distinct cell populations that share the expression of a few common markers.

Here we aim to answer some of these questions and begin to parse the potential heterogeneity of muscle-resident mesenchymal cell populations. To better understand FAP populations in the context of developing muscle, we chose to work with postnatal mice at 21 days of age (P21). This juvenile time point sits at the transition from development to adulthood, which allows for insight into both adult and developmental cell populations. Utilizing a Prrx1Cre; R26-tdTomato reporter mouse line, fluorescent activated cell sorting (FACS), single-cell RNA sequencing (scRNA-seq), immunofluorescence (IF), differentiation assays, and muscle injury models, we identify six distinct Pdgfra-expressing non-myogenic muscle-resident mesenchymal cell populations within juvenile skeletal muscles. Our study provides unique insights into these muscle-resident non-myogenic mesenchymal cell populations while providing a means for both their identification and isolation, which will have broad implications for future skeletal muscle research.

RESULTS

Muscle-resident mesenchymal cells are distinct from myogenic and other cell lineages during muscle development and regeneration

The Prrx1Cre (also known as Prx1Cre) mouse line has been widely used to study LPM-derived mesenchymal lineages (Logan et al., 2002), often in bone and connective tissues of the limb skeleton. Therefore, we sought to confirm its utility in identifying mesenchymal lineages within skeletal muscle tissue by utilizing a Prrx1Cre; R26-tdTomato reporter mouse line, hereafter referred to as Prrx1;R26-tdT. We first examined the distribution of fluorescent tdTomato-positive (tdT+) cells at various time points throughout musculoskeletal development. At embryonic day 11.5 (E11.5), tdT+ cells are largely restricted to the LPM and limb bud mesenchyme within whole embryos (Figure 1A). By E16.5, tdT+ cells are observed broadly within the muscle interstitium and are associated with limb skeletal elements (Figures 1B1 and B2). As the skeletal muscle approaches homeostasis following early postnatal development (P21), tdT+ cells can be identified within the muscle interstitium and surrounding the blood vessels (Figure 1C). The tdT+ cells persist in these locations as the mice reach adult stages (4 months of age, 4M) (Figure 1D). When the skeletal muscle was challenged with barium chloride (BaCl2) injury or was aged (12 months of age, 12M), the tdT+ mesenchymal cells remained within the skeletal muscle interstitium and did not contribute to the myofibers within the medial aspects of the muscles studied (Figures 1E and 1F). Flow cytometry confirmed that tdT+ cells are not muscle stem cells (Pax7+ satellite cells), endothelial cells (CD31+), or immune cells (CD45+, F4/80+) (Figures 1G and S1). Taken together, these results demonstrate that the Prrx1Cre-expressing cells and descendants represent a class of skeletal-muscle-resident mesenchymal cells that are in the interstitial space of muscle fibers.

Figure 1.

Muscle-resident mesenchymal cells are distinct from myogenic and other cell lineages during muscle development and regeneration

(A) Whole embryo fluorescent/bright-field imaging of Prrx1Cre;TdT at E11.5.

(B1) Whole hindlimb cross-section of Prrx1Cre;TdT at E16.5.

(B2) Muscle interstitial area from (B1) at increased magnification.

(C–F) Tibialis anterior muscle cross-sections of Prrx1Cre;TdT at (C) P21, (D) 4 months, (E) 3 weeks post BaCl2 injury, performed at 4 months of age, and (F) 1 year of age. Immunofluorescence for tdTOMATO (tdT), DAPI, and PHALLOIDIN marking F-ACTIN myofiber bundles. n ≥ 3 per time point for all image assessments.

(G) Flow cytometry of Pax7GFP and TdT from Prrx1Cre;TdT;Pax7 hindlimb muscle cells performed at P21. n = 4.

Single-cell RNA sequencing exposes the heterogeneity of muscle-resident mesenchymal cells

To investigate the heterogeneity of muscle-resident mesenchymal cells, we performed scRNA-seq analysis on sorted tdT+ cells isolated from P21 Prrx1;R26-tdT hindlimb muscles that were liberated of tendon tissue (Figure 2A). Utilizing the 10x Genomics platform and Seurat-R workflow, 2,480 cells were analyzed via unsupervised clustering, revealing nine independent cell clusters (Stuart et al., 2019) (Figure 2B). All clusters expressed high levels of tdTomato and some level of Prrx1 expression, demonstrating the purity of the sorted cells and the continued expression of Prrx1 in many muscle-resident mesenchymal cells (Figure S2). Each cluster can be identified by unique markers, exemplified by the expression of some of the most highly and uniquely expressed genes in each cluster (Kcnj8, Cnn1, Scx, Osr1, Adam12, Gap43, Clu, Hsd11b1, Gli1) (Figures 2C and 2D; Table S1). Below we describe the individual clusters in greater detail regarding each of their unique gene signatures and functions.

Figure 2.

scRNA-seq reveals nine clusters of tdTomato+ muscle-resident mesenchymal cells

(A) Graphical representation of experimental workflow for scRNA-seq. Hindlimb muscles liberated from tendon tissue were harvested from P21 Prrx1Cre;TdT mice, digested to single cells, and isolated via FACS. 2,480 tdT+ cells were sequenced using the 10x Genomics scRNA-seq platform.

(B) Unsupervised clustering in Seurat returns nine clusters from single-cell transcriptomic analysis.

(C) Violin plots for genes uniquely expressed in each cluster.

(D) Feature plots for genes uniquely expressed in each cluster.

Clusters 1–3 exhibit a gene expression signature consistent with known muscle-resident mesenchymal cell populations. The transcriptional profile of cluster 1 includes high levels of Regulator of G-protein signaling 5 (Rgs5), Nestin (Nes), and Potassium inwardly rectifying channel subfamily j member 8 (Kcnj8) (Figures S3A and S3C), which are established pericyte-associated genes. The cells of cluster 2 are vascular smooth muscle cells and express genes associated with mature vascular smooth muscle, including: Smooth muscle actin alpha2 (Acta2), Transgelin (Tagln), and Calponin1 (Cnn1) (Figures S3B and S3C). The identity of both clusters was validated by IF, where some of the tdT+ cells co-express α-smooth muscle actin and localize in close proximity to CD31+ endothelial cells (Figures S3D and S3E). Previously, these perivascular cells in muscle have been grouped together and were described as smooth muscle mesenchymal cells (Giordani et al., 2019); however, our analysis indicates that these are distinct groups. Cluster 3 is highly enriched in the expression of tenocyte-associated genes including Scleraxis (Scx), Tenomodulin (Tnmd), and Mohawk (Mkx) as well as WNT1-inducible signaling pathway protein 1 (Wisp1) (Figures S4A and S4B). We validated by IF the co-expression of TENOMODULIN and tdT within the skeletal muscle interstitium at sites distant to the tendon (Figure S4C). This cell type has previously been identified as a unique muscle interstitial tenocyte-like population (Giordani et al., 2019).

Additional muscle-resident mesenchymal cells are currently described as PDGFRα+, SCA-1+, and/or HIC1+ FAPs or TCF4+ muscle connective tissue fibroblasts (CTFs) (Joe et al., 2010; Mathew et al., 2011; Scott et al., 2019; Uezumi et al., 2010, 2011). When assessing FAP gene expression in our scRNA-seq dataset, we identified high expression of Pdgfra in clusters 4–8, high expression of Lymphocyte antigen 6 (Ly6a; which encodes SCA1) in clusters 5–7 with more restricted expression in clusters 4 and 8, and high Hic1 expression in clusters 1 and 4–9 (Figures 3A and 3B). We also identified the CTF-associated gene, Transcription factor 7 like 2 (Tcf7l2; which encodes TCF4), to be predominantly expressed within clusters 4–9 (Figures 3A and 3B). We validated these results via IF showing co-localization of PDGFRα with tdT (Figure 3C). When examined by flow cytometry, nearly all (94%) PDGFRα+/SCA1+ cells were also positive for tdT (Figures 3D and 3E). However, just over 30% of the PDGFRα+ cells are both lineage negative (lin−; CD31−, CD45−, and F4/80−) and tdT−. Given the strong recombination efficiency of the Prrx1Cre, these results open the possibility that some of the PDGFRα+/tdT− cells within the muscle interstitium may arise from outside the LPM (Seo and Serra, 2007) (Figures 3D and 3F). While these results demonstrate that Prrx1Cre-expressing cells and descendants include FAP and CTF populations, our data suggest that these broad classes of cells may represent distinct subpopulations of non-myogenic muscle-resident mesenchymal cells.

Figure 3.

scRNA-seq exposes the heterogeneity of muscle-resident mesenchymal cells

(A and B) Single-cell (A) feature and (B) violin plots for fibro-adipogenic genes Pdgfra, Ly6a, Hic1, and Tcf7l2.

(C) Immunofluorescence staining for PDGFRα, TdT, and F-ACTIN (stained by PHALLOIDIN) in a cross-section from P21 Prrx1Cre;TdTf/+ tibialis anterior muscles. n ≥ 3.

(D–F) Flow cytometry of Prrx1Cre;TdT postnatal hindlimb muscle cells for lineage markers (DAPI/CD31/CD45/F480), SCA1, PDGFRα, and tdT. n = 5.

Osr1 expression identifies a primitive muscle-resident mesenchymal cell population

Cluster 4 has high and concentrated expression of Odd-skipped related transcription factor 1 (Osr1) (Figures 4A and 4B). Osr1 is expressed in developing muscle-resident mesenchymal cell populations and is thought to mark a precursor population for adult FAPs (Stricker et al., 2012; Stumm et al., 2018; Vallecillo-Garcia et al., 2017). Previous utilization of the Osr1-CreERT2 has shown that OSR1+ populations give rise to PDGFRα+ cells in muscle formation, and that OSR1+ cells increase and give rise to PDGFRα+ FAPs after injury (Stricker et al., 2012; Stumm et al., 2018; Vallecillo-Garcia et al., 2017). Our scRNA-seq analysis found that many genes associated with stemness and/or developmental processes were specifically enriched within cluster 4, including Gata6, Stc1, and Nog (Chaturvedi et al., 2009; Fisher et al., 2017; Li et al., 2018; Rifas, 2007; Whissell et al., 2014; Yoon et al., 2018). Taken in the context of the established understanding of OSR1 in FAP development, this expression profile further supports the role of cluster 4 as a precursor population (Figures 4A and 4B). To examine the trajectories of our cell clusters we performed RNA velocity, which computationally assesses the differentiation trajectory of cells over developmental time (La Manno et al., 2018). Consistent with these cells being identified as progenitors, cluster 4 cells localize to the base of our RNA velocity analysis, indicating that it is a primitive FAP precursor population (Figures 4C and 4D). Interestingly, the two observed trajectories parallel a recently suggested, but not functionally verified, FAP trajectory which diverges into separate Dipeptidylpeptidase 4 (Dpp4+) and C-X-C motif chemokine ligand 14 (Cxcl14+) arms (Oprescu et al., 2020) (Figures 4E and 4F). In our RNA velocity analysis, cluster 4 gives rise to the Dpp4+ and Cxcl14+ trajectories that can each be further divided into more discrete and distinct cell populations. The Dpp4+ arm includes cells of clusters 6 and most of cluster 5, and the Cxcl14+ arm includes cells of clusters 7 and 8/9 (Figures 4C–4F).

Figure 4.

Osr1 expression identifies a primitive muscle-resident mesenchymal cell population

(A and B) Single-cell (A) feature and (B) violin plots for Osr1, Gata6, Stc1, and Nog.

(C) RNA velocity analysis of Pdgfra-expressing cells.

(D) Key to RNA velocity analysis corresponding to Pdgfra+ uniform manifold approximation and projection (UMAP) clusters.

(E and F) (E) Feature and (F) violin plots for Dpp4 and Cxcl14.

Adam12 and Gap43 expression identifies immune-responsive muscle-resident mesenchymal cell populations

Clusters 5 and 6 demonstrate unique and high expression of ADAM metallopeptidase domain 12 (Adam12) and Growth associated protein 43 (Gap43), respectively (Figures 5A and 5B). Beyond these unique markers, clusters 5 and 6 have enriched expression of genes associated with fibrosis and transforming growth factor β (TGF-β) signaling including Transforming growth factor β receptors 2 (Tgfbr2) and Fibronectin 1 (Fn1), with less enriched expression of Latent transforming growth factor β binding protein 4 (Ltbp4), a negative regulator of TGF-β signaling whose expression is associated with decreased fibrosis (Figures S5A and S5B) (Biernacka et al., 2011; Ceco et al., 2014; Delaney et al., 2017; To and Midwood, 2011). These clusters also express adipogenic progenitor and precursor genes, including Dpp4, Peptidase inhibitor 16 (Pi16), Annexin a3 (Anxa3), Intercellular adhesion molecule 1 (Icam1), Peroxisome proliferator activated receptor γ (Pparg), CCAAT enhancer binding protein α (Cebpa), and Fatty acid binding protein 4 (Fabp4) (Figures S6A and S6B). These gene expression profiles highlight the fibrogenic and adipogenic potential of these cells.

Figure 5.

Adam12 and Gap43 expression identify immune-responsive muscle-resident mesenchymal cell populations

(A and B) Single-cell (A) feature and (B) violin plots for cluster 5 (Adam12) and cluster 6 (Gap43) genes and type 2 immune-receptor genes (Il4ra, Il13ra1).

(C) LipidTOX/DAPI-stained PDGFRα+ sorted (C1) ADAM12+ cluster 5, (C2) GAP43+ cluster 6, (C3) CLU+ cluster 7, (C4) HSD11B1 clusters 8/9 cells cultured in adipogenic medium for 5 days.

(D) LipidTOX/DAPI-stained PDGFRα+ sorted (D1) ADAM12+ cluster 5, (D2) GAP43+ cluster 6, (D3) CLU+ cluster 7, (D4) HSD11B1 clusters 8/9 cells cultured in adipogenic medium supplemented with recombinant IL-4 for 5 days.

(E) LipidTOX/DAPI-stained PDGFRα+ sorted (E1) ADAM12+ cluster 5, (E2) GAP43+ cluster 6, (E3) CLU+ cluster 7, (E4) HSD11B1 clusters 8/9 cells cultured in adipogenic medium supplemented with recombinant IL-13 for 5 days.

(F) Quantification of adipogenic differentiation as visualized by %LipidTOX/DAPI. Data are represented as mean ± SEM. Intergroup analysis performed with two-way ANOVA and post hoc Tukey’s test. Intragroup analysis performed with paired Student’s t test. n ≥ 9 per cluster and condition. Significance denoted by an asterisk.

(G) Flow-cytometry contour plots for PDGFRα and either IL-4Ra (cluster 5/6), CLU (cluster 7), or HSD11B1 (cluster 8/9) gated on tdT+ live cells obtained from BaCl2-injured muscle and NaCl-injected contralateral controls at 4 days post injury as well as glycerol-injured muscle and NaCl-injected contralateral controls 3 days post injury.

(H) Percent quantification of tdT+/PDGFRα+/IL4Ra (cluster 5/6), tdT+/PDGFRα+/CLU+ (cluster 7), and tdT+/PDGFRα+/HSD11B1 (cluster 8/9) cells in BaCl2-injured muscle or NaCl-injected contralateral controls. n ≥ 6 per cluster and condition. Data are represented as mean ± SEM. Significance determined using two-way ANOVA and post hoc Tukey’s test. Significance denoted by an asterisk. NS, not significant.

(I) Percent quantification of tdT+/PDGFRα+/IL4-Ra (cluster 5/6), tdT+/PDGFRα+/CLU+ (cluster 7), and tdT+/PDGFRα+/HSD11B1 (cluster 8/9) cells in glycerol-injured muscle or NaCl-injected contralateral controls. Data are presented as mean ± SEM. n ≥ 5 per cluster and condition. Significance determined using two-way ANOVA and post hoc Tukey’s test. Significance denoted by an asterisk. NS, not significant.

All scale bars represent 200 μm.

To examine the fibrogenic capacity of these cells, we sorted clusters 5–9 based on their expression of PDGFRα and specific cluster type markers (ADAM12, GAP43, CLU, or HSD11B1) and plated the resulting cells. Cells were treated with fibrogenic media supplemented with recombinant TGF-β for 5 days. Despite differences in expression patterns in Fn1, Tgfb2, and Ltbp4 within clusters 5–9 there was no significant difference in the expression of Fn1 or Col1a1 in the sorted clusters following 5 days of fibrogenic medium (Figure S5C). This result suggests that cells from each cluster have comparable fibrotic capacity. However, it does not determine when different clusters may be triggered into a more fibrotic condition. To explore this further, we delve deeper into the unique attributes of each cluster.

Our scRNA-seq shows that clusters 5 and 6 exhibit high expression of the genes encoding the interleukin-4 (IL-4) (Il4ra) and interleukin 13 (IL-13) (Il13ra1) receptors (Figures 5A and 5B). The IL-4 and IL-13 receptors are associated with the type 2 immune response in skeletal muscle, facilitating FAP responses during acute muscle injury and preventing FAP- mediated intramuscular adipogenesis (Heredia et al., 2013). The expression of immune receptors in specific muscle-resident mesenchymal cell populations at this stage suggests that there is already a pre-existing group of cells primed to respond to a regenerative muscle injury. To determine the response of clusters 5 and 6 to type 2 immune cytokines, we sorted clusters 5–9 based on their expression of PDGFRα and specific cluster type markers (ADAM12, GAP43, CLU, or HSD11B1) and plated the resulting cells. Cells were treated with adipogenic media ± recombinant IL-4 or IL-13 (rIL-4, rIL-13) cytokines for 5 days and the amount of lipid formed was examined using LipidTOX, which stains for the accumulation of neutral lipid. The ADAM12+ cluster 5 and GAP43+ cluster 6 accumulated significantly more lipid (21.40% and 16.55%) than either CLU+ cluster 7 or HSD11B1+ cluster 8/9 (11.66% and 7.35%) in control conditions (Figures 5C1–C4 and 5F). Following treatment with rIL-4 or rIL-13 cytokines, the ADAM12+ cluster 5 and GAP43+ cluster 6 are the only cells that respond to type 2 immune cytokines, as evidenced by a decrease in adipogenesis of cluster 5 and 6 cells treated with type 2 immune cytokines (Figures 5D–5F). Interestingly, the ADAM12+ cluster 5 cells preferentially respond to rIL-4 and not rIL-13. However, GAP43+ cluster 6 cells respond to both rIL-4 and rIL-13 (Figures 5C–5F), which is likely a reflection of the different expression levels of Il4ra and Il13ra1 in clusters 5 and 6 (Figures 5A and 5B).

Given the importance of the type 2 immune signaling in acute muscle injury, we next sought to determine whether IL-4/IL-13 responsive cells (clusters 5/6) react differently to acute muscle injury in comparison with clusters 7 and 8/9. Four days following acute muscle damage via intramuscular BaCl2 injection, we performed flow cytometry for tdT+/PDGFRα+/IL4Ra+ cells (clusters 5/6), tdT+/PDGFRα+/CLU+ cells (cluster 7), and tdT+/PDGFRα+/HSD11B1+ cells (clusters 8/9) to examine their response to injury. We observed a significant increase in TdT+/PDGFRα+/IL4Ra+ cells (cluster 5/6), in contrast to the minimal change observed from cells in clusters 7–9 (Figures 5G and 5H). To further examine the ability of clusters 5 and 6 to respond in acute injury contexts that elicit intramuscular fat accumulation, we performed intramuscular glycerol injections. Glycerol injuries result in intramuscular adipogenic deposition as opposed to the more fibrotic damage observed following BaCl2 injury (Mahdy et al., 2015). Three days following acute muscle damage by intramuscular glycerol injection, we performed flow cytometry for tdT+/PDGFRα+/IL-4Ra+ cells (clusters 5/6), tdT+/PDGFRα+/CLU+ cells (cluster 7), and tdT+/PDGFRα+/HSD11B1+ cells (clusters 8/9) to examine their response to injury (Figures 5G and 5I). Once again, we observed that clusters 5/6 show a significant increase in TdT+/PDGFRα+/IL4-Ra+ cells, contrasting with little to no increase in clusters 7 (CLU+) or 8/9 (HSD11B1+). Collectively, these data suggest that the ADAM12+ and GAP43+ cells (clusters 5/6; IL-4Ra+ cells) are primed and preferentially respond to acute skeletal muscle injury as compared with other PDGFRα+ muscle-resident mesenchymal cells.

Beyond the increased response of cluster 5/6 cells to acute muscle injury, these cells also demonstrate an increased presence in situations of chronic muscle injury, specifically a mouse model of muscular dystrophy. To illustrate this, we performed flow cytometry on D2-Mdx mice and D2 littermate controls. For this experiment, we utilized DPP4 to label cluster 6 and part of cluster 5, ADAM12 for cluster 5 more exclusively, CLU for cluster 7, and HSD11B1 for cluster 8/9 (Figure S7A). Due to the enhanced adipogenic and fibrotic responses that occur in D2-Mdx muscles, our flow-cytometry analysis demonstrates an increase in CD31−/CD45−/PDGFRα+/ADAM12+ (cluster 5) and CD31−/CD45−/PDGFRα+/DPP4+ (cluster 5/6) cells within D2-Mdx mice and compared with controls. By contrast, we also observed a decrease in CD31−/CD45−/PDGFRα+/CLU+ (cluster 7) and CD31−/CD45−/PDGFRα+/HSD11B1+ (clusters 8/9) when compared with D2 controls (Figure S7B). In sum, these results demonstrate that clusters 5 and 6 preferentially expand in chronic injury settings such as muscular dystrophy, while cells from clusters 7 and 8/9 decrease.

Clu expression identifies a mineralizing muscle-resident mesenchymal cell population

Cluster 7 exhibits unique expression of Clusterin (Clu) and Hemicentin 1 (Hmcn1), among other genes. These cells are also enriched in the expression of genes known to regulate mineralization and osteogenesis, including Ectonucleotide pyrophosphatase/phosphodiesterase 1 (Enpp1) and Fam20C (Figures 6A and 6B) (Liu et al., 2018; Nam et al., 2011; Vogel et al., 2012). Previous studies have demonstrated that muscle-resident mesenchymal cells, specifically PDGFRα+ cells, are capable of mineralizing and do so in diseases of HO and certain models of neuromuscular disease (Eisner et al., 2020; Leblanc et al., 2011; Lees-Shepard et al., 2018; Oishi et al., 2013; Wosczyna et al., 2012). To assess the mineralization capacity of the cells in clusters 5–9, we sorted tdT+/PDGFRα+ cells (control), tdT+/PDGFRα+/ADAM12+ cells (cluster 5), tdT+/PDGFRα+/GAP43+ cells (cluster 6), tdT+/PDGFRα+/CLU+ cells (cluster 7), and tdT+/PDGFRα+/HSD11B1+ cells (clusters 8 and 9). All cell populations were treated with osteogenic medium for 7 days and then stained for Alizarin red (Figure 6C). The CLU+ cluster 7 and to a minor extent HSD11B1+ cluster 8/9 are capable of mineralizing. Only the mineralization of cluster 7 cells is significant and comparable with tdT+/PDGFRα+ (clusters 4–9) control. Neither the ADAM12+ cluster 5 nor the GAP43+ cluster 6 cells are capable of mineralizing (Figures 6C and 6D).

Figure 6.

Clu expression identifies a mineralizing muscle-resident mesenchymal cell population

(A and B) Single-cell (A) feature and (B) violin plots for cluster 7-specific gene markers (Clu, Hmcn1) and mineralizing genes (Enpp1, Fam20C).

(C) Sorted clusters (C1) tdT+/PDGFRα+ clusters 4–9, (C2) tdT+/PDGFRα+/ADAM12+ cluster 5, (C3) tdT+/PDGFRα+/GAP43+ cluster 6, (C4) tdT+/PDGFRα+/CLU+ cluster 7, (C5) tdT+/PDGFRα+/HSD11B1 clusters 8/9 cells cultured in osteogenic medium for 7 days and then stained with Alizarin red.

(D) Quantification of Alizarin red staining relative to tdT+/PDGFRα+ sorted cells. Data are presented as mean ± SEM. n ≥ 4 for each cell type. Significance denoted by an asterisk as calculated by one-way ANOVA.

PDGFRα+ cells have been implicated in the disease FOP, a rare genetic condition caused by activating mutations in the BMP type I receptor gene, ACVR1 (Lees-Shepard et al., 2018; Shore et al., 2006). FOP is a debilitating disorder wherein muscle is transformed to heterotopic bone after even the slightest of soft tissue injuries (Pignolo et al., 2011). Therefore, we assessed the expression of BMP receptors (Acvra1 and Bmpr1a), signaling components (Smad 1/5), and downstream target genes (Id1–Id4) in our scRNA-seq analysis. While BMP receptors are widely expressed across all PDGFRα clusters, the downstream BMP genes show an enrichment of expression within the CLU+ cluster 7 (Figures S8A and S8B). Additionally, gene ontology (GO) analysis of clusters 4–9 demonstrates that terms indicating mineralization and bone formation are more commonly associated with cluster 7, while terms indicative of an immune response and fat are more enriched in the GO analysis of clusters 5 and 6 (Figure S9). To determine the ability of PDGFRα+ cell clusters to respond to BMP signaling and subsequently mineralize, we sorted tdT+/PDGFRα+ cells (control), tdT+/PDGFRα+/ADAM12+ cells (cluster 5), tdT+/PDGFRα+/GAP43+ cells (cluster 6), tdT+/PDGFRα+/CLU+ cells (cluster 7), and tdT+/PDGFRα+/HSD11B1+ cells (clusters 8 and 9). All cell populations were treated with osteogenic media ± recombinant BMP2 for 7 days and then stained for Alizarin red (Figure S8C). Similar to Figure 6, only the cells of CLU+ cluster 7 and HSD11B1+ cluster 8/9 are capable of mineralizing. ADAM12+ cluster 5 and GAP43+ cluster 6 remain incapable of mineralizing even after supplementation with recombinant BMP2 (Figures S8C and S8D). Collectively, these data indicate that only the CLU+ cluster 7 and HSD11B1+ cluster 8/9 cells contribute in a direct manner to the mineralization and bone formation that occurs within PDGFRα+ cells.

Hsd11b1 expression identifies a neuromuscular junction-associated muscle-resident mesenchymal cell population

Clusters 8 and 9 are uniquely identified by the expression of Hydroxysteroid 11β dehydrogenase 1 (Hsd11b1) and many genes that are associated with neuromuscular junctions (NMJs) and neuromuscular diseases, including GDNF family receptor alpha 1 (Gfra1) and Ret and membrane metalloendopeptidase (Mme), among others (Figures 7A and 7B) (Auer-Grumbach et al., 2016; Baudet et al., 2008; Nguyen et al., 1998; Paratcha and Ledda, 2008). IF imaging at P21 for MME confirms their co-localization with tdT+ cells and the α-bungarotoxin+ (αBTX+) NMJ (Figures 7C–7E). To investigate the response of clusters 8 and 9 to neuromuscular injury, we performed sciatic nerve transection followed by IF and flow cytometry at 1, 2, and 4 weeks post injury (WPI). IF staining demonstrates that MME+/tdT+ cells, surrounding αBTX+ NMJ, expand following nerve transection at 1, 2, and 4 WPI (Figures 7F, 7G, S10A, and S10B). Flow cytometry further revealed that the TdT+/PDGFRα+/HSD11B1+ cell population (cluster 8/9) is significantly expanded at 1, 2, and 4 weeks post sciatic nerve injury as compared with contralateral controls (Figures 7H4, 7I4, 7J, 7K, S10C, and S10D). The ability of cluster 8/9 cells to respond to nerve injury was also compared with the immune-responsive clusters 5 and 6 and the mineralizing cluster 7 populations. Flow cytometry revealed that none of the TdT+/PDGFRα+/ADAM12+ (cluster 5), TdT+/PDGFRα+/DPP4+ (cluster 5/6), or TdT+/PDGFRα+/CLU+ (cluster 7) cell populations were significantly expanded at 1 or 2 weeks following sciatic nerve injury (Figures 7H1–7H3, 7J, S10C, and S10D). However, by 4 WPI nearly all populations were significantly expanded as compared with contralateral controls, revealing a wider fibrotic or mesenchymal response (Figures 7I1–7I4 and 7K). Overall, these results establish that the HSD11B1+ cluster 8/9 cells are a population of NMJ-associated muscle-resident non-myogenic mesenchymal cells with an early response to peripheral nerve injury.

Figure 7.

Hsd11b1 expression identifies an NMJ-associated muscle-resident mesenchymal cell population

(A and B) Single-cell (A) feature and (B) violin plots for cluster 8/9-specific (Hsd11b1) and neuromuscular (Gfra1, Ret, Mme) gene markers.

(C–E) Immunofluorescence on P21 Prrx1Cre;tdT tibialis anterior muscle cross-sections for: (C) αBUNGAROTOXIN (αBTX) to mark NMJs, tdT, and DAPI; (D) MME, tdT, DAPI, and F-ACTIN; (E) αBTX, MME, and DAPI.

(F and G) Immunofluorescence of adult Prrx1Cre;tdT tibialis anterior muscle cross-sections following sciatic nerve transection for αBTX, tdT, MME, and DAPI in contralateral control and sciatic nerve injury muscle at (F1, F2) 1 week post injury (WPI) and (G1, G2) 4 WPI.

(H) Histogram from flow cytometry at 1 WPI showing (H1) ADAM12+, (H2) DPP4+, (H3) CLU+, and (H4) HSD11B1+ cells within PDGFRα+/tdT+ live cells isolated from nerve-injured and contralateral control muscles.

(I) Histogram from flow cytometry at 4 WPI showing (I1) ADAM12+, (I2) DPP4+, (I3) CLU+, and (I4) HSD11B1+ cells within PDGFRα+/tdT+ live cells isolated from nerve-injured and contralateral control muscles.

(J and K) Percent quantifications of HSD11B1+ cells within the PDGFRα+/tdT+ population from nerve-injured and contralateral control muscles at (J) 1 WPI and (K) 4 WPI. Data are presented as mean ± SEM. n ≥ 5. Significance denoted by an asterisk as calculated by one-way ANOVA with post hoc Tukey’s test. NS, not significant. All scale bars represent 25 μm.

Of note, cluster 9 is a small cell population that has a largely overlapping gene signature compared with cluster 8, which includes enriched expression of NMJ-related genes such as Gfra1 and Mme. Therefore, these two clusters have been linked in our functional analysis. However, differences are observed between these two cell populations/states. Cluster 9 exhibits relatively unique expression of Interleukin 15 (Il15), while both clusters 9 and 4 show enriched expression of Hedgehog (Hh) signaling genes that include Gli1, Patched1 (Ptch1), and Patched2 (Ptch2) (Figures S11A and S11B). Notably, cluster 9 has very low levels of FAP marker expression, including Pdgfra and Ly6a (Figure S11B). These data suggest that cluster 9 may be a more differentiated subset of the NMJ-associated cluster 8 cells, which are Hh-responsive and have generally downregulated the expression of traditional FAP markers.

Single-cell RNA sequencing and functional analyses identify novel non-myogenic muscle-resident mesenchymal cell populations that differentially respond to muscle injuries

Our scRNA-seq and functional analyses, described above, have identified nine distinct non-myogenic muscle-resident mesenchymal cell populations within juvenile skeletal muscle. These include three previously identified interstitial niche cell populations (pericytes, vascular smooth muscle cells [SMCs], and tenocyte-like cells), as well as six Pdgfra-expressing cell populations that fit into a bipartite differentiation trajectory originating from a common Osr1+ progenitor. This differentiation trajectory branches toward either immune-responsive, fibrogenic, and highly adipogenic FAP populations or NMJ-associated and/or mineralizing fibroblast populations with limited adipogenic capacity and unique responsiveness to peripheral nerve injury and Hh signaling (Figure S12).

To further validate the existence of these non-myogenic muscle-resident mesenchymal cell populations in other contexts, we compared our data with published scRNA-seq datasets of muscle-associated cell populations in adult skeletal muscle at homeostasis and following muscle injury. To ensure that we were investigating comparable “FAP-like” cell populations across datasets, we first subset the published results based on Pdgfra or Gelsolin (Gsn) expression due to the fact that these datasets were obtained from the entirety of skeletal muscle cell populations. We next integrated the Pdgfra+ populations of our scRNA-seq data with existing FAP datasets from uninjured adult muscle (De Micheli et al., 2020; Oprescu et al., 2020). Our Pdgfra+ subpopulations of cells integrated well within these previous works (Figures S13A and S13B). When multiple cluster-associated genes (clusters 4–9) were utilized to detect these populations, we were able to identify the presence of nearly all of our defined subpopulations (Figure S13C). Of note, cluster 5 (ADAM12+) cells from our study of juvenile skeletal-muscle-associated cells were not identified to a significant degree in either published scRNA-seq dataset obtained from adult muscle under homeostatic conditions.

To further orient our work in the context of muscle injury, we first examined how the acute and chronic muscle injury and FAP-associated genes, Vcam1, Lox, Tek, and Myoc, identified previously (Malecova et al., 2018), were expressed among our non-myogenic muscle-resident mesenchymal cell populations. We observe that Vcam1 and Lox are more enriched in our immune-responsive clusters 5 and 6, with some expression in the primitive Osr1+ (cluster 4) cell population. Tek primarily shows balanced expression between cluster 6 and our NMJ-associated cell population (cluster 7). Myoc, on the other hand, exhibits enrichment primarily in mineralizing and NMJ-associated fibroblasts (clusters 7–9) (Figures S14A and S14B). We next verified the dynamic changes in the expression of these FAP-related genes and their associated populations in the scRNA-seq dataset of Oprescu et al. (2020) over the course of 21 days following a cardiotoxin muscle injury. Our analyses indicate that Vcam1- and Lox-expressing cells (mostly clusters 4–6) increase within 5 days following cardiotoxin injury and ultimately return to baseline numbers by 21 days post injury, while MyoC- and Tek-expressing cells (mostly clusters 7–9) decrease and return to baseline over the same time course (Figure S14C). Consistent with these trends, we observe that our more cluster-specific genes Adam12, Il4ra, and Il13ra (clusters 5/6; immune-responsive) demonstrate an increase in frequency or expression by 3.5–5 days post injury that return to baseline by 21 days post injury. Meanwhile, the mineralizing and NMJ-associated fibroblasts (clusters 7–9), marked mostly by Hsd11b1 here, show a decrease in frequency or expression, which ultimately returns closer to baseline levels by 21 days post injury (Figure S14D). This increase in immune-responsive FAPs following cardiotoxin muscle injury strongly resembles the cellular response we observed following both BaCl2 and glycerol acute muscle injuries (Figures 5H and 5I); however, the decreasing populations of NMJ-associated or mineralizing fibroblasts (clusters 7–9) are more reminiscent of what we observed in the D2-Mdx model of chronic injury and muscle degeneration/regeneration (Figures S7A and S7B). These various muscle injury models are known to have distinct modes of action and cellular responses that likely account for these differential changes in our non-myogenic muscle-resident mesenchymal cell populations (Hardy et al., 2016). Finally, it is also of note that Adam12 expression returns to this immune-responsive population of adult non-myogenic muscle-resident mesenchymal cells only following muscle injury (Figures S13C and S14D), indicating that it may be developmentally regulated or represents a transient state that ultimately gives rise to a more permanent or long-lived Dpp4-, Gap43-, Il4ra-, and Il13ra-expressing immune-responsive FAP population.

DISCUSSION

In this study we highlight the heterogeneity of LPM-derived muscle-resident mesenchymal cell populations utilizing a Prrx1Cre;R26-tdTomato reporter mouse line, FACS, scRNA-seq, IF, differentiation assays, and muscle injury models. Our genetic reporter studies spanning embryonic and postnatal development, as well as adult muscle in both homeostatic and acute muscle injury settings, demonstrate that the LPM-derived mesenchymal cells and other Prrx1Cre-expressing cells and their descendants do not give rise to the myogenic lineage to form myofibers within skeletal muscles. These data were further confirmed via flow cytometry. Our scRNA-seq analysis identified non-myogenic muscle-resident mesenchymal cell populations that included both non-PDGFRα+ cells (pericytes, vSMCs, and interstitial tenocyte-like cells) and PDGFRα+ cells. Based on our flow-cytometry studies, our data suggest the presence of a PDGFRα-expressing cell population(s) within skeletal muscles that may not originate from the LPM, or at least is not marked by the expression of the Prrx1Cre transgene. Liu et al. (2017) identified a TWIST2+/PDGFRα+/PDGFRβ+ muscle-resident mesenchymal cell population capable of contributing directly to muscle, which they also suggested may not be derived from the LPM. Our data open the possibility that this population is represented within our PDGFRα+/tdT-population of muscle-resident mesenchymal cells, which we did not otherwise directly study (Liu et al., 2017). Other groups have also identified muscle-resident mesenchymal cells of an LPM origin that contribute to muscle fibers specifically at the myotendinous junctions (Esteves de Lima et al., 2021; Yaseen et al., 2021). While our work used similar genetic approaches, we did not observe this population because we removed tendon and myotendinous tissues and exclusively examined muscle-resident mesenchymal populations at muscle regions distant to the tendons.

At present, PDGFRα-expressing muscle-resident cells are classified as FAPs, implying both a fibrogenic and adipogenic fate capacity. This study advances our understanding of FAP identity, demonstrating unique and distinct populations that comprise the LPM-derived PDGFRα+ cells within skeletal muscles. Consistent with prior work, we find that PDGFRα+ cells diverge into two trajectories that can be identified by the expression of either Dpp4 or Cxcl14, among other genes (Oprescu et al., 2020). Unlike their analysis, we were able to place a primitive Osr1+ population of cells at the root of these trajectories while also determining specific aspects of the functionality, identity, and localization of cell populations within each trajectory. One trajectory includes two populations of fibro-adipogenic cells which are immune-responsive and highly adipogenic, and which expand following acute or chronic muscle injury. These populations of immune-responsive fibro-adipogenic cells are also incapable of mineralizing even in the presence of exogenous BMP signaling. In comparison, a second trajectory includes cells which are similarly capable of fibrosis, but with limited adipogenic capacity. This trajectory includes two cell populations prone to mineralization, with at least one of these populations being NMJ-associated and responding quickly to neuromuscular injury that also contains a subset of Hh-responsive cells (Figure S12).

The identification of an immune-responsive subset of muscle-resident mesenchymal cells in an uninjured state enhances our understanding of their many roles in muscle regeneration. Earlier studies have indicated that PDGFRα+ and/or SCA1+ cells are vital in muscle regeneration by regulating key immune responses (Heredia et al., 2013; Kang et al., 2018). Our work indicates that IL-4 and IL-13 signaling only within these cell populations likely co-ordinates the proliferation and differentiation injury response and prevents their aberrant intramuscular adipogenesis following muscle damage (Heredia et al., 2013; Schiaffino et al., 2017). Here, we advance previous work by identifying specific type 2 immune-responsive populations. We demonstrate that not all muscle-resident mesenchymal PDGFRα+ cells express IL-4Ra, nor do all PDGFRα+ cells respond to IL-4 or IL-13 signaling. Instead, there exist only two populations of PDGFRα+ muscle-resident mesenchymal cells that express Il4ra and Il13ra1 (ADAM12+ cluster 5 and GAP43+ cluster 6). These groups are both competent to respond to IL-4 signaling but differentially respond to IL-13 signaling. In addition to demonstrating differences in IL-4/IL-13 signaling, this work demonstrates the existence of PDGFRα+ subgroups that do not robustly respond to acute muscle injury (BaCl2 or glycerol injury) and decrease in chronic muscle injury (D2-Mdx). This significant finding supports the importance of immune-responsive subgroups of non-myogenic muscle-resident mesenchymal cells in muscle injury, since the IL-4Ra+ cells are the predominant PDGFRα+ cell type in regenerative acute and chronic muscle injuries.

In addition to the identification of immune-responsive FAPs, our work identifies a distinct group of muscle-resident mesenchymal cells that are prone to mineralization (primarily CLU+ cluster 7). PDGFRα+ cells are the pathologic cell type associated with FOP and non-genetic traumatic HO (Eisner et al., 2020; Lees-Shepard et al., 2018; Wosczyna et al., 2012). As such, these findings may have significant implications for our understanding of these diseases. Beyond pathologic mineralization, other work has suggested a role for muscle-resident mesenchymal cells in bone repair, showing their ability to respond during bone healing (Julien et al., 2021). The existence of a muscle-resident mesenchymal population prone to mineralization identified by our study may have broad implications for the isolation and application of such cells in bone repair.

Cluster 8/9 cells, the trajectory partner to cluster 7, are physically identified as an NMJ-associated cell population. The presence of an NMJ-associated cell population has been hypothesized previously, sometimes referred to as kranocytes (Court et al., 2008). Kranocytes were primarily characterized via electron microscopy as an NMJ-capping cell with a theorized function of regulating neuromuscular activity. While some have shown that NMJ-associated fibroblasts are beneficial (Uezumi et al., 2021), others suggested a pathogenic role for NMJ-associated fibroblasts in the context of neuromuscular degenerative diseases, such as amyotrophic lateral sclerosis (Gatchalian et al., 1989; Gonzalez et al., 2017; Madaro et al., 2018). The origin of such a cell type and the precise molecular makeup are presently not well understood. Our study begins to bridge this gap by identifying, and defining at the transcriptional level, an NMJ-associated muscle-resident mesenchymal cell population that is likely derived from the LPM. Furthermore, we have demonstrated that this cell type responds to neuromuscular injury, expanding following sciatic nerve injury. The specific increase of the NMJ-associated cell population in response to nerve injury, but not BaCl2 or glycerol injuries, and their propensity to decrease in settings of chronic injury and neuromuscular disease (muscular dystrophy; D2-Mdx model) demonstrates the functional differences observed between distinct non-myogenic muscle-resident mesenchymal cell populations.

Muscle-resident mesenchymal cells are an ill-defined heterogeneous group of cells; however, they are a dynamic and critical member of the muscular niche. In this study, we identify distinct cell populations and differentiation trajectories that emanate from a common Osr1+ progenitor within PDGFRα-expressing muscle-resident mesenchymal cells. While one trajectory contains two distinct populations that robustly form adipocytes and respond to specific immune cytokine signaling, the second trajectory contains two cell populations demonstrating limited adipogenic potential, with one population demonstrating strong functional capacity to mineralize and another population that localizes to NMJs and responds to neuromuscular damage. Given the current understanding of FAPs as “fibro-adipogenic” progenitor cells encompassing all PDGFRα+ muscle-resident mesenchymal cells, our work described here comes into conflict with this current definition, and therefore suggests that it may be time to re-examine how FAPs are classified to account for the existence of these distinct PDGFRα+ non-myogenic muscle-resident mesenchymal cell populations.

Limitations of the study

This work relies on the use of the Prrx1Cre mouse line, which exhibits Cre expression in most cells of the LPM and many of its descendants. While this model has a high recombination efficiency, it may not capture all descendants of the LPM. Therefore, it is difficult to determine whether the non-myogenic muscle-resident mesenchymal cells that have escaped our analysis (tdTomato−; PDGFRα+) are descended from the LPM and can be categorized into similar cell populations or not.

Another important but necessary technical limitation of our work is the use of singular genes to label or identify each cell population. All scRNA-seq clusters are not defined by the expression of single genes, but rather from a more complex summation of gene expression that creates a unique genetic signature per cell. To capture each population in the most direct way possible, we chose markers which were most specifically restricted to, and enriched within, the populations of interest. This allows us to develop and test tools (antibodies or genetic reporters) for the isolation and characterization of these largely unique populations. We recognize that this technique may not capture the full complement of cells within each cluster; however, utilizing population markers to study each group of non-myogenic muscle-resident mesenchymal cells is powerful and the best available methodology for the validation and functional analyses of these cells.

STAR★METHODS

RESOURCE AVAILABILITY

Lead contact

Further questions regarding information and requests for resources should be directed to the lead contact, Matthew Hilton (matthew.hilton@duke.edu).

Materials availability

No new materials were generated during this study. All resources and reagents are commercially available and the source and accession or catalogue numbers can be found in the STAR resource table.

Data and code availability

Raw data and the code utilized in this work has been deposited into the GEO server with the identifier GEO: GSE200234. These files will be publicly available as of the date of publication. Processed data files have been included as Table S1. No original code was written to analyze the single-cell RNA-seq data, instead the Seurat platform was utilized (as described in greater detail in the method details section). Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS

Mouse lines

The Prrx1Cre mice (Stock No. 005584), Rosa26-tdTomato mice (Ai9 Stock No. 007909), and D2-Mdx mice (Stock No. 013141) were obtained from Jackson Laboratories. The Pax7-eGFP mice were graciously provided by Dr. David Kirsch. Both sexes of mice were used in this study and they were housed at 23°C on a 12-h light/dark cycle and maintained on a PicoLab Rodent Diet. All animal work was approved by the Duke University Institutional Animal Care and Use Committees (IACUC). No sex differences were noted. scRNA-sequencing was performed on a female mice. All other experiments included both male and female animals in sex matched cohorts.

Barium chloride injury

Solid barium chloride (BaCl2) was dissolved in sodium chloride (NaCl) at a concentration of 1.2% and then filter sterilized. Animals were put under anesthesia using isoflurane before a 25g transdermal needle and syringe pre-loaded 50uL of 1.2% BaCl2 was used to administer BaCl2 throughout the tibalias anterior muscle. Animals recovered from anesthesia under observation.

Glycerol injury

Since a 50% glycerol solution produces robust fat deposition (Mahdy et al., 2015), we performed acute muscle injuries using 50uL of 50% glycerol. Prior to injury, 50% glycerol by volume in HBSS was prepared and 50uL was pre-loaded into a syringe affixed with a 25G transdermal needle. Animals were put under anesthesia using isoflurane and then intramuscular injections were performed. Animals recovered from anesthesia under observation.

Sciatic nerve injury

Prior to injury animals were administered Buprenorphine HCl for pain and were placed under anesthesia using isoflurane. A trained veterinarian (D.R.) prepared the animal, opening the skin of the upper hindlimb, and severing the sciatic nerve. The contralateral side was used as an uninjured control.

METHOD DETAILS

Tissue preparation and immunostaining

Harvest and preparation

Freshly isolated Tibalias Anterior muscles were washed in 1X PBS and then incubated overnight in 30% sucrose before imbedding in Shandon Cryomatrix (ThermoFisher Scientific). Samples were flash frozen in isopentane, supercooled in a liquid nitrogen bath. Cryosections were cut at 10μM or thick 50μM and stored at −20°C until staining. Sections of the muscle were taken from the middle of the muscle, sectioned after the end of the myotendinous junction.

Immunostaining

Slides were brought to room temperature and dehydrated in PBS. Following a 5-min post fixation with 4% PFA, antigen retrieval was performed as shown in Table 1. For unconjugated antibodies, slides were blocked with 1% goat serum for 1 h before incubating with primary antibodies at 4°C in a humified chamber. Slides were washed and incubated with secondary antibody for 45 min at room temperature. Sections were mounted with ProLong Gold Mountant with DAPI and imaged on a Leica DMI3000B microscope using a Leica DFC3000G camera or a Leica SP5 inverted confocal microscope with the latest Leica imaging suite. Images were reconstructed and pseudocolored on Fiji.

Muscle resident cell isolations

Muscle interstitial cells were released from various hindlimb muscles using a previously described protocol (Latroche et al., 2018). Hindlimb muscles were isolated and mechanically disrupted, then enzymatically digested in a volume of 1mL per mouse containing 10mg/mL Collagenase B and 2.4U/mL Dispase for 30 min at 37°C in a rotating oven. Samples were shaken vigorously every ten minutes during the 30-min digestion.

Flow cytometry and Fluorescence-activated cell sorting (FACS)

Following digestion described above, the cell suspension was passed through a 70μM filer and centrifuged at 350 x g for 7 min. After red blood cell lysis cells, were stained on ice in the dark. If staining with unconjugated primary antibodies, cells were first stained with primary antibodies on ice for 30 min in PBS +10% FBS. Staining was quenched with 5–10mL of PBS+10% FBS, depending on staining volume, and cells were spun for at 350 x g for 7 min. Cells were then stained for 30 min with secondary and conjugated antibodies in a mastermix. Staining was quenched with 5–10mL PBS+10% FBS and spun for 350 x g for 7 min. Cells were resuspended in PBS+2% PBS and taken for FACS, or flow analysis. Cells sorted from FACS were sorted into 0.5mL DMEM+10%FBS and 1% PenStrep.

Cell culture and cellular staining

In vitro cell culture of sorted clusters

Following sorting, cells recovered in growth media (DMEM+10%FBS, and 1%PenStrep) for multiple days until they reached full plate coverage. Once cells reached full plate coverage, but were still slightly sub-confluent, they were switched to differentiation media.

Adipogenic culture

Fresh adipogenic media was made for each group. Adipogenic media was composed of growth media as previously described supplemented with adipogenic supplement (STEMCELL technologies Cat#05507) then filter sterilized. For cytokine treatment assays, adipogenic media was supplemented with 1 ng/uL of recombinant IL4 or 1 ng/mL of recombinant IL13, reconstituted in PBS, and filter sterilized before use. Cells were treated with adipogenic media, or IL4/IL13 supplemented adipogenic media for 5 days before analysis.

Fibrogenic culture

Fresh fibrogenic media was made for each group. Fibrogenic media was composed of growth media supplemented with recombinant TGF β at a concentration of 0.25ng/mL. Cells were treated with fibrogenic media for 5 days prior to analysis.

LipidTOX staining

HCS LipidTOX was used to assess the adipogenic potential of plated clusters. Cells were washed and then fixed with 4% paraformaldehyde for 10 min. After 3 washes following fixation, cells were stained with DAPI for 15 min to allow for cellular visualization. Following DAPI, cells were stained with LipidTOX for at least 30 min. Plates are imaged on a ZEISS 880 Airyscan Inverted Confocal using a plate adaptor.

Osteogenic culture

Osteogenic media was generated from growth media supplemented with 50ug/mL Ascorbic Acid, 100nM Dexamethosone, and 10mM Beta-glycerophosphate. BMP2 supplemented osteogenic media was made using osteogenic media with recombinant BMP2 at a concentration 100ng/mL then filter sterilized. Cells were treated with osteogenic media, or osteogenic media + rBMP2 for 7 days before analysis.

Alizarin Red Staining and quantification

Osteogenic cultured cells were washed and fixed with formalin for 10 min before staining with 2% Alizarin red. Following a 15-min incubation, Alizarin red staining was washed with distilled water. Alizarin red stain was measured as absorbance after 10% Cetylpryridinium Chloride (CPC) destaining.

Single-cell isolation and scRNA-seq

10x Transcriptome library prep

Single muscle interstitial cells were isolated as described above. Following red blood cell lysis, cells were stained with DAPI for 30 min on ice. Cells were washed with PBS+10%FBS, spun down and resuspended in PBS +2% FBS for fluorescent activated cell sorting (FACS). Cells were gated on forward side scatter, side scatter, DAPI and TOMATO, to yield single live tomato cells as shown in Figure 2. As with similar work, previously described Long et al. (Long et al., 2022), cells were taken to the Duke Molecular Physiology Institute Molecular Genomics core where they were combined with a master mix that contained reverse transcription reagents. The gel beads carrying the Illumina TruSeq Read 1 sequencing primer, a 16bp 10x barcode, a 12bp unique molecular identifier (UMI) and a poly-dT primer were loaded onto the chip, together with oil for the emulsion reaction. The Chromium Controller partitions the cells into nanoliter-scale gel beads in emulsion (GEMS) within which reverse-transcription occurs. All cDNAs within a GEM, i.e. from one cell, share a common barcode. After the RT reaction, the GEMs were broken and the full-length cDNAs cleaned with both Silane Dynabeads and SPRI beads. After purification, the cDNAs were assayed on an Agilent 4200 TapeStation High Sensitivity D5000 ScreenTape for qualitative and quantitative analysis.

Enzymatic fragmentation and size selection were used to optimize the cDNA amplicon size. Illumina P5 and P7 sequences, a sample index, and TruSeq read 2 primer sequence were added via End Repair, A-tailing, Adaptor Ligation, and PCR. The final libraries contained P5 and P7 primers used in Illumina bridge amplification. Sequence was generated using paired end sequencing (one end to generate cell specific, barcoded sequence and the other to generate sequence of the expressed poly-A tailed mRNA) on an Illumina NextSeq500 with 136M read per run.

Computational analysis of scRNA-seq data

Analysis

The primary analytical pipeline for the SC analysis followed the recommended protocols from 10X Genomics. Briefly, we demultiplexed raw base call (BCL) files generated by Illumina sequencers into FASTQ files, upon which alignment to the mouse reference transcriptome, filtering, barcode counting, and UMI counting were performed using the most current version of 10X’s Cell Ranger software. We used the Chromium cell barcode to generate feature-barcode matrices encompassing all cells captured in each library. The secondary statistical analysis was performed using the last R package of Seurat (Stuart et al., 2019). In Seurat, data were first normalized and scaled after basic filtering for minimum gene and cell observance frequency cut-offs (200–4000 genes, 7% mitochondrial genes, genes expressed in >5 cells). We then closely examine the data and performed further filtering based a range of metrics in attempt to identify and exclude possible multiplets (i.e. instances where more than one cell was present and sequenced in a single emulsified gel bead). We thresholded tomato at 2 expression or higher to exclude a extremely minor amount of cells that were FACS sorted but were tomato negative by gene expression.

After quality control procedures were complete, we performed linear dimensional reduction calculating principal components using the most variably expressed genes in our dataset (2000 variable genes, dims = 30). The genes underlying the resulting principal components are examined in order to confirm they are not enriched in genes involved in cell division or other standard cellular processes (subsetting out percent.mito). Significant principal components for downstream analyses are determined through methods mirroring those implemented by Macosko et al., and these principal components were carried forward for two main purposes: to perform cell clustering and to enhance visualization (Macosko et al., 2015). Cells were grouped into an optimal number of clusters for de novo cell type discovery using Seurat’s FindNeighbors() and FindClusters() functions (resolution = 0.4), graph-based clustering approach with visualization of cells being achieved through the use of UMAP, which reduces the information captured in the selected significant principal components to two dimensions (Stuart et al., 2019). Differential expression of relevant cell marker genes was visualized on UMAP plots to reveal specific individual cell types.

RNA velocity

RNA velocity of Pdgfra-expressing cells was computed using velocyto.py and velocyto.R (La Manno et al., 2018). The aligned BAM file was processed using velocyto.py v0.17 to obtain the counts of unspliced and spliced reads in loom format. The loom file was processed using velocyto.R v0.6 in combination with an R package SeuratWrappers v0.3.0 (https://github.com/satijalab/seurat-wrappers). Twenty nearest neighbors in slope calculation smoothing were used for RunVelocity().

QUANTIFICATION AND STATISTICAL ANALYSIS

When comparing two groups, statistical analysis was performed using a two-tailed, unpaired Student’s t test unless otherwise noted. A p value of <0.05 was considered significant. When dealing with batched animals in cytokine treated sorted cells, paired Student’s t test were used to account for the observed batch effects. This difference is noted in the corresponding figure legend. When comparing multiple groups, a two-way ANOVA with post-hoc Tukey’s test was used, as noted in the figure legend.

KEY RESOURCES TABLE

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
anti-CD31-FITC	Biolegend	Cat#102405; RRID: AB_312900
anti-CD45-FITC	Biolegend	Cat#103107; RRIS: AB_314005
anti-F4/80-FITC	Invitrogen	Cat#11-4801-82; RRID: AB_2538124
anti-Ly6A/E-APC	Biolegend	Cat#108111
anti-CD140a-PE/Cy7	Biolegend	Cat#135911
goat anti-rabbit IgG (H + L), APC	ThermoFisher	Cat#A-10931
APC anti-mouse CD124	Biolegend	Cat#144807; RRID: AB_2750451
anti-HSD11B1	Abcam	Cat#ab39364; RRID: AB_940037
anti-ADAM12	ProteinTech	Cat#14139-1-AP; RRID:AB_2880689
anti-GAP43	Invitrogen	Cat#PA1-16729; RRID: AB_568546
anti-CLU	Abcam	Cat#ab69644; RRID:AB_1267705
anti-CD26	Invitrogen	Cat#MA5-32643; RRID:AB_1955200
DAPI	ThermoFisher	Cat#D1306; RRID:AB_2629482
anti-αSMA	Abcam	Cat#ab5694; RRID:AB_2305167
anti-RFP	Abcam	Cat#ab62341; RRID: AB_945213
anti-TNMD	Abcam	Cat#ab203676; RRID: AB_2722782
anti-CD10	Invitrogen	Cat#MA5-14050; RRID:AB_10983979
αBungarotoxin, Alexa Fluor 488	ThermoFisher	Cat#B13422
αBungarotoxin, Alexa Fluor 647	ThermoFisher	Cat#B35450
Alexa Fluor 647 Phalloidin	ThermoFisher	Cat#A22287
Alexa Fluor 488 goat anti-rabbit IgG	ThermoFisher	Cat#A32731
Alexa Fluor 594 goat anti-rabbit IgG	ThermoFisher	Cat#A32740
Alexa Fluor 488 goat anti-mouse IgG	ThermoFisher	Cat#A32723
Alexa Fluor 594 goat anti-mouse IgG	ThermoFisher	Cat#A32742
ProLong Gold Antifade Mountant with DAPI	Invitrogen	Cat#P36935
HCS LipidTOX Neutral Lipid Stain	ThermoFisher	Cat#H34475
Chemicals, peptides, and recombinant proteins
Recombinant Mouse IL4	R&D Systems	Cat#404-ML
Recombinant Mouse IL13	R&D Systems	Cat#413-ML
Recombinant Mouse BMP2	R&D Systems	Cat#335-BM
Recombinant Mouse TGFβ	R&D Systems	Cat#7666-MB-005
DMEM	ThermoFisher (Gibco)
HBSS	ThermoFisher (Gibco)
Dexamethasone	Sigma-Aldrich	Cat#D1756
Ascorbic Acid	Sigma-Aldrich	Cat#A4544
B-Galactosidase	Sigma-Aldrich	Cat#G9422
Collagenase B	Millipore Sigma (Roche)	Cat#110888
Dispase II	Millipore Sigma (Roche)	Cat#4942078001
Barium Chloride	Millipore Sigma (Sigma Aldrich)	Cat#B0750
Glycerin	Millipore Sigma (Sigma Aldrich)	Cat#G2289
Critical commercial assays
RNeasy Mini Kit	Qiagen	Cat#74104
MesenCult Adipogenic Differentiation Kit	STEMCELL Technologies	Cat#05507
Deposited data
De Micheli AJ, Laurilliard EJ, Heinke CL, Ravichandran H et al. Single-Cell Analysis of the Muscle Stem Cell Hierarchy Identifies Heterotypic Communication Signals Involved in Skeletal Muscle Regeneration. Cell Rep2020 Mar 10; 30(10):3583-3595.e5. PMID: 32160558	GEO	GSE143437
Oprescu SN, Yue F, Qiu J, Brito LF et al. Temporal Dynamics and Heterogeneity of Cell Populations during Skeletal Muscle Regeneration. iScience 2020 Apr 24;23(4):100993. PMID: 32248062	GEO	GSE138826
This paper	GEO	GSE200234
Experimental models: Organisms/strains
Mouse: Prrx1Cre	Jackson Laboratories	005584
Mouse: Rosa26-tdTomatof/+	Jackson Laboratories	Ai9 007909
Mouse: Pax7-eGFP	Dr. David Kirsch (Duke University)
Mouse: D2.B10-Dmdmdx/J (D2-Mdx)	Jackson Laboratories	013141
Oligonucleotides
Fn1-F	Invitrogen	GCGACTCTGACTGGCCTTAC
M-GAPDH-F	IDT	GCACAGTCAAGGCCGAGAAT
M-GAPDH-R	IDT	GCCTTCTCCATGGTGGTGAA
Fn1-R	Invitrogen	CCGTGTAAGGGTCAAAGCAT
M-Col1a1-F	IDT	GCATGGCCAAGAAGACATCC
M-Col1a1-R	IDT	CCTCGGGTTTCCACGTCTC
Software and algorithms
Seurat	Hao et al., “Integrated analysis of multimodal single-cell data”, Cell (2021), Vol. 184, Issue 13. Pages 3573-3587 e.29 ISSN 0092-8674
RNA Velocity	La Manno et al., “RNA Velocity of single cells”, Nature (2018), Vol. 560, Issue 7719, Pages 494-498