A reference single-cell transcriptomic atlas of human skeletal muscle tissue reveals bifurcated muscle stem cell populations

1Skeletal muscles are essential to daily functions such as locomotion, respiration, and metabolism. Upon damage, resident muscle stem cells (MuSCs) repair the tissue in coordination with supporting non-myogenic cell types such as immune cells, fibroblasts, and endothelial cells. However, with age and disease, the repair capacity of MuSCs declines, leading to complications such as fibrotic scarring, reduced muscle mass and strength, fat accumulation, and decreased insulin sensitivity, all of which severely affect mobility and quality of life.

2Human MuSCs are defined by the expression of the paired box family transcription factor PAX7 and can be isolated using various surface marker proteins including β1-integrin (CD29), NCAM (CD56), EGFR, and CD82 to varying purities. With aging, human MuSCs exhibit a heterogeneous expression of the senescence marker p16Ink4a and accumulate other cell-intrinsic alterations in myogenic gene expression programs, cell cycle control, and metabolic regulation. However, given their varied molecular and functional states, our understanding of MuSCs in adult human muscle tissue remains incompletely defined. In addition, cellular coordination in the regulation of human muscle homeostasis and regeneration remains poorly understood due to the lack of experimentally tractable models with multiple human muscle cell types. Given these challenges, we posited that an unbiased single-cell reference atlas of skeletal muscle could provide a useful framework to explore MuSC variability and communication in adult humans.

3Here, we deeply profiled the transcriptome of thousands of individual MuSCs and muscle-resident cells from diverse adult human muscle samples using single-cell RNA-sequencing (scRNA-seq). After integrating these donor datasets to conserve biological information and overcome technical variation, we resolved two subpopulations of MuSCs with distinct gene expression signatures. Using differential gene expression analysis and ligand-receptor interaction modeling, we extend the known repertoire of human MuSC gene expression programs, suggesting new regulatory programs that may be associated with human MuSC activation, as well as features of human muscle aging and disease.

4We used scRNA-seq to collect and annotate a single-cell transcriptomic dataset of diverse adult human muscle samples under homeostatic conditions. The muscle samples were from surgically discarded tissue from n = 10 donors (range 41 to 81 years old) undergoing reconstructive procedures and originating from a wide variety of anatomical sites in otherwise healthy patients (Fig. 1a). Each sample was ~ 50 mg after removal of extraneous fat and connective tissue. Muscle samples were enzymatically digested into single-cell suspensions and independently loaded into the 10X Chromium system. All together, we collected over 22,000 human muscle single-cell transcriptomes (2206 ± 1961 cells per dataset) into a single data compendium. Using unsupervised clustering, we resolved 16 types of cells of immune, vascular, and stromal origin, as well as two distinct subpopulations of MuSCs and some myofiber myonuclei (Fig. 1b).

5Given important differences in anatomical site, donor health history, age, sex, and surgical procedures, the muscle samples were highly heterogeneous in terms of cell-type diversity and underlying gene expression profiles. Comparing the resulting scRNA-seq datasets is therefore a challenge that we addressed using recently developed bioinformatic integration methods. Our goal was to assemble a unified dataset of human muscle tissue that faithfully conserved sources of biological variability such as donor, anatomical location, and cell composition heterogeneity, while accounting for technical biases. We tested four different scRNA-seq data integration methods (Fig. S1 and S3) and found that Scanorama followed by scaling the output by regressing against the library chemistry technical variable (“10X chemistry”) and the number of genes detected per single-cell best satisfied this goal. Detailed information on our methodology is provided in Fig. S1. After integrating the 10 datasets, we noted remarkable consistency amid cell types across donors (Fig. 1c, e), owing to the robustness of scRNA-seq technology, the bioinformatic method chosen, and our sample preparation protocol. Differential gene expression analysis between the 16 distinct subpopulations identified an extensive set of unique markers that we grouped into 4 categories (Fig. 1d).

6We annotated and interpreted the consensus cell atlas (Fig. 1b, d) into cell type subpopulations as follows. We identify four types of stromal cells starting with adipocytes found to be expressing apolipoprotein D (APOD)), the brown fat tissue adipokine CXCL14, GPX3, and GLUL. Among the 3 other subpopulations of fibroblast-like cells, Fibroblast 1 expresses high levels of collagen 1 (COL1A1), SFRP4, SERPINE1, and CCL2; Fibroblast 2 expresses fibronectin (FBN1), the microfibril-associated glycoprotein MFAP5, and CD55 known to be expressed by synoviocytes; and Fibroblast 3 is mainly characterized by SMOC2 identified in tendon fibroblasts. The Fibroblast 3 cluster is similar to the adipocytes cluster though exhibits lower expression levels and frequencies of the marker genes APOD, CXCL12, and GLUL, and contain pre-adipocytes.

7We also identify 5 types of vascular cells, including 3 endothelial subpopulations, and a subpopulation of pericytes and smooth muscle cells (SMCs). Pericytes and SMCs express the canonical markers RGS5 and MYH11. Endothelial 1 express E-selectin (SELE), IL6, ICAM1, and VCAM1. These genes are upregulated at sites of inflammation to facilitate immune cell recruitment, suggesting this Endothelial 1 cell population may be involved in homeostatic muscle tissue remodeling. Endothelial 2 cells are distinguished by expressing high levels of claudin-5 (CLDN5), ICAM2, and the chemokine CXCL2. Endothelial 3 expresses high levels of the platelet-recruiting Von Willebrand Factor (VWF) and caveolin-1 (CAV1), a protein known to regulate cholesterol metabolism, atherosclerosis progression, and MuSC activation. Endothelial 3 cells are enriched for expression of BTLN9, suggesting they might represent a lymphatic endothelial phenotype.

8We also noted two types of myeloid immune cells: first, tissue-resident and anti-inflammatory macrophages which express CD74 and histocompatibility complex HLA proteins; second, activated macrophages and monocytes that express inflammatory markers such as S100A9 (calgranulin) and LYZ (lysozyme). Moreover, S100A9 transcript abundance levels have been shown to be a feature in aging and chronic inflammation. We also identified a pool of T/B lymphocytes and NK cells characterized by IL7R and NKG7, respectively, as well as a small subset of HBA1+ erythroblasts.

9Finally, we identified two subpopulations of MuSCs (henceforth called “MuSC1” and “MuSC2”). MuSC1 highly expressed the canonical myogenic transcription factor PAX7, as well as chordin-like protein 2 (CHRDL2) and Delta-like non-canonical Notch ligand 1 (DLK1). CHRDL2 has been shown previously to be expressed in freshly isolated quiescent human MuSCs, though its function is still to be understood. DLK1 is an inhibitor of adipogenesis whose role in muscle has mainly been recognized in the embryo but remains controversial in adult muscle regeneration. In contrast to MuSC1, MuSC2 expressed lower levels of PAX7 but maintain expression of MYF5 (a marker of activated MuSCs) and APOC1 (Fig. 2b). Interestingly, the MuSC2 population also had elevated expression of two long non-coding RNAs (lncRNAs), LINC00152, and MIR4435-2HG. LncRNAs are involved in regulating myogenesis. Surprisingly, we detected low expression of the myogenic commitment factors MYOD1 and MYOG (Fig. 2b), in contrast to scRNA-seq analyses of adult mouse muscle. These observations suggest that the MuSC1 and MuSC2 populations are both comprised largely of muscle stem cells, not committed myogenic progenitors. In addition, we noted that “Myonuclei” population (Fig. 1b) was enriched for myosin light chain (MYLFP), skeletal alpha-actin (ATCA1), and troponin C (TNNC2), proteins involved in muscle contraction. This multiple-donor scRNA-seq atlas highlights the cellular diversity of human muscle tissue and revealed two distinct MuSC subpopulations along with specific myogenic expression programs.

10We examined genes that were differentially expressed between the MuSC1 and MuSC2 subpopulations and the biological processes that characterize them (Fig. 2a, b). The MuSC1 subpopulation was enriched for PAX7, DLK1, and CHRDL2, as well as for the cyclin-dependent kinase inhibitor CKDN1C (encoding P57KIP2), suggesting that these cells are quiescent and not cycling. In addition, this subpopulation expresses the transcription factor BTG2, which was identified in mouse to be enriched in quiescent MuSCs. We also note that the MuSC1 subpopulation expressed elevated levels of mitochondrial genes as well as FOS, JUN, and ERG1. Upregulation of these genes has been shown to be potential artefacts of the enzymatic digestion during the sample preparation.

11The MuSC2 subpopulation was enriched for multiple markers of inflammation including CCL2, CXCL1, IL32, and surface receptor TNFRSF12/FN14. In particular, CCL2 and CXCL1 are inflammatory cytokines known to be upregulated in muscle repair, exercise, and fat metabolism. In addition, IL32 has been shown to have inflammatory properties in human obesity and have a negative impact on insulin sensitivity and myogenesis, while TNFRSF12/FN14 has been implicated in various muscle wasting diseases and metabolic dysfunction. Furthermore, the MuSC2 population is enriched for ribosomal gene expression (e.g. RPLP1 and RPS6; data not shown), indicating that these cells may have elevated translational mechanisms. Lastly, the MuSC1 population has enriched expression of the myogenic gene PAX7 and, to a lesser extent, MYF5, compared the MuSC2 population. These observations suggest that MuSC1 is comprised of quiescent MuSCs, and MuSC2 is comprised of an early-activated MuSCs.

12We performed Ingenuity Pathway Analysis (IPA) to compare biological processes differentially activated between the MuSC1 and MuSC2 populations. The IPA gene group “Oxidative Phosphorylation” is enriched in MuSC1, while “EIF2 Signaling,” associated with protein translation processes, is enriched in MuSC2 (Fig. 2c). Furthermore, Gene Set Enrichment Analysis (GSEA) also found MuSC1 to be enriched for “myogenesis,” “muscle cell differentiation,” “hypoxia,” and “response to mechanical stimulus” gene sets, supporting the observation that these cells are both less differentiated and may exhibit enhanced transcriptional responses to mechanical disruption due to tissue dissociation (Fig. 2d). MuSC2 cells are enriched for “ribosome and translational initiation,” “MYC targets,” and “E2F (cell proliferation),” “G2M checkpoint (cell division),” and “inflammation” gene sets, further supporting the interpretation that these cells may be in an early activated or partially differentiated state within an inflammatory environment (Fig. 2d). Taken together, these observations suggest that the MuSC1 population is comprised of quiescent MuSCs, while the MuSC2 population is comprised of active, proliferating, and/or dysregulated MuSCs, with expression alterations associated with inflammation, aging, and muscle wasting. Differentially expressed genes such as IL32, CXCL1, CCL2, and TNFRSF12/FN14 may constitute a marker set for MuSC variation in chronic muscle inflammation in various pathologies.

13We used a ligand-receptor (LR) interaction model and a database of LR pairs to map cell signaling communication channels in human muscle and uncover differences between MuSC1 and MuSC2 subpopulations (Fig. 3). The model also identifies interacting ligand(s) and is restricted to receptor genes differentially expressed by a specific cell type within the consensus human muscle cell atlas (Fig. 1b). For each LR pair, the model calculates an interaction score from differentially expressed receptors on a given cell population (e.g., “MuSC1”) relative to all other population and ligands expressed by other cell types. The MuSC1 and MuSC2 subpopulations are involved in numerous LR interactions, as both ligand- and receptor-expressing cells (Fig. 3a), though a majority of all LR interaction pairs instead involve other cell types. This suggests that only a small subset of potential paracrine interactions in human muscle may include MuSCs.

14Given the distinct expression profiles between the MuSC1 and MuSC2 populations, we sought to identify genes that could facilitate surface antigen-based separation of these two human MuSC populations for prospective isolation strategies. We identified surface receptor genes that were differentially expressed between the MuSC1 and MuSC2 populations, using a database of 542 human surface “receptor” genes (Fig. 3). MuSC1 exhibit elevated expression of EGFR, ITGB1, FGFR4, SDC2, as well as the three tetraspanins CD81, CD82, and CD151(Fig. 3b). EGFR is a recently established human MuSC marker and is required for basal-apical asymmetric cell division. The tetraspanin CD82 is also a recently recognized human MuSC maker, while CD9 and CD81 have been identified to control muscle myoblast fusion. Furthermore, Syndecans (SDCs) have been identified in mouse to be heterogeneously expressed on MuSCs and myoblasts during muscle repair and have been shown to form co-receptor complexes with integrin β1 (ITGB1) and FGFR4 upstream of signaling pathways regulating myogenesis. Only SDC4 and SDC3 have yet been identified to mark adult mouse MuSCs. In comparison, the MuSC2 subpopulation has elevated expression of CD44 and TNFRSF12/FN14 as previously noted (Fig. 3b). The CD44 receptor has been shown to regulate myoblast migration and fusion in mouse, but also mark MuSCs inosteoarthritis patients.

15Next, we focused the LR analysis on the MuSC1 and MuSC2 populations. We identified 73 and 6 significant LR interactions for the MuSC1 and MuSC2 populations, respectively (Fig. 3c). Over one third of all interactions in the MuSC1 subpopulation involve the EGFR receptor, which has recently been shown to play a critical role in directing MuSC asymmetric division in regenerating muscle. A limited number of EGFR ligands have been identified in muscle repair, for example, amphiregulin. (AREG) secreted by Treg cells. According to our model findings, EGFR may also interact with ligands expressed by immune cells, such as with TGF-α (TGFA), heparin-biding EGF (HBEGF), amphiregulin (AREG), and epiregulin (EREG). Other EGFR ligands include brevican (BCAN), and betacellulin (BTC) produced by endothelial cells; ECM proteins fibulin 3 (EFEMP1), decorin (DCN), and tenascin C (TNC) expressed by fibroblasts; and FGF13, AHM, NRG4, and EGF, expressed by mature skeletal myofibers. We also detect seven interactions involving NOTCH3 with a variety of ligands. Notch3 signaling is involved in maintaining MuSC quiescence, in particular through interaction with DLL4, which we found differently expressed by endothelial cells along with JAG2. In addition, NOTCH3 also interacts with the ECM protein thrombospondin-2 (THBS2).

16Only two receptors, TNFRSF12/FN14 and RPSA, were found differentially expressed in MuSC2 compared to other cell types. The first, TNFRSF12/FN14, interacts with the TWEAK cytokine ligand. While typically recognized to be expressed by macrophages and other immune cells, our model suggests that TWEAK is also expressed by the Fibroblast 2 and pericyte cell populations, though not in a statistically significant manner. The second, RPSA, is surface ribosomal protein that interacts with laminins (LAM), a dual-specificity phosphatase 18 (DUSP18), and prion protein 2 (PRND), which taken together may suggest various pathological processes such as prion diseases and cancer. Together, this ligand-receptor analysis identifies a broad set of surface markers that could refine the molecular definition of human MuSCs and their subpopulations, as well as candidate cell-communication channels differentially involved in healthy and diseased muscle tissues.

17Lastly, we performed a comparative analysis of receptor gene expression between mouse and human MuSCs. We integrated the human scRNA-seq datasets described in Fig. 1 and an adult mouse muscle injury-response scRNA-seq time-course previously reported by converting mouse genes to their corresponding human ortholog. The multi-species scRNA-seq atlas was integrated with Scanorama and corrected with Harmony (Fig. S2A-B). From this integrated atlas, we annotated all clusters as in Fig. 1. We identified two MuSC clusters which both contained cells from both mouse and human samples. We then performed differential expression analysis between species comparing aggregated human MuSC1 and MuSC2 cells to mouse MuSCs from the uninjured timepoint (Fig. S2C). We found that EGFR and CD99 were most differentially expressed by human MuSCs and, conversely, CRLF1 and SDC4 were most enriched in mouse MuSCs. This findings suggest that mouse and human MuSC exhibit species-specific receptor expression signatures.

18Here we present an annotated multi-donor single-cell RNA sequencing dataset consisting of 22,000 single-cell transcriptomes from 10 different donors and unique anatomical sites, some of which difficult to access outside of reconstructive surgeries. Our study complements other recent reports by Rubenstein et al. and Barruet et al., which collected dissociated whole vastus lateralis muscles and FACS-sorted MuSC samples mostly from vastus lateralis muscles, respectively, by providing more diversity in anatomical sites and donor demographics. As such, these scRNA-seq data exhibited notable biological and technical variation, and therefore, we applied the bioinformatic method Scanorama to assemble an integrated cellular atlas with minimal technical biases so that we could examine the cellular heterogeneity across diverse adult human muscle tissue samples. We observed that Scanorama performed more successfully than other data integration approaches, especially when including a scaling regression for sequencing chemistry (Fig. S1 and S3). Notably, even after performing Scanorama with scaling, we still observed that integrated atlas exhibited biological (donor) and technical (sequencing chemistry) biases, but retained some degree of donor-specific cell-type subpopulations.

19We describe the muscle tissue cellular heterogeneity and provide a comprehensive analysis of differentially expressed genes for 16 resolved cell subpopulations (Fig. 1), adding to a growing literature documenting human muscle cell transcriptional diversity. Compared to other studies, the broader variety of muscle tissue samples combined with the lack of FACS selection allowed us to identify candidate subpopulations of muscle fibroblasts and vascular endothelial cells that may provide unique perspective to human muscle physiology. In particular, we remark that Endothelial 1 expressed DARC/ACKR1, a gene identified in mouse and human to mark cells of post-capillary venular origin (Fig. 1d). Rubenstein et al. also found a DARC/ACKR1+ post-capillary venular endothelial cluster and a second VWF+ FABP+ cluster, which overlaps with the Endothelial 2 and 3 clusters reported here. We suggest that the Endothelial 2 cluster may contain both arterial and capillary endothelial cells, but could not further partition and classify this cluster. We suggest that the Endothelial 3 cluster may represent lymphatic endothelium due to its differential expression of BTLN9, a marker of lymphatic endothelial cells.

20Most notably, this analysis suggests that human muscle may contain two distinct MuSC subpopulations (Fig. 2). This finding contrasts with Rubenstein et al. which observed a single MuSC (“satellite cell”) population from dissociated whole muscle samples and Barruet et al. which observed ~ 12 clusters from human MuSCs prospectively enriched by CXCR4+/CD29+/CD56+ FACS. Since cluster distinction depends on both the cellular diversity and sample complexity, it is expected that variation in study design and methods will yield differing conclusions regarding sub-population resolution. In this work, we found a “MuSC1” subpopulation to be largely comprised of “quiescent” MuSCs, owning to high levels of PAX7, the mitotic inhibitor CDKN1C, and DLK1. Interestingly, DLK1 may be an important regulator for human MuSC maintenance and a marker of healthy tissue given its role in inhibiting adipogenesis. Conversely, we identified in the “MuSC2” population signatures of inflammation and increased fat metabolism (CCL2 and CXCL1), reduced insulin sensitivity (IL32), cell cycle (EIF2 Signaling terms), and muscle wasting (TNFRSF12/FN14), thereby suggesting that these cells may constitute an “early-activated” and possibly dysfunctional MuSC pool. These markers are consistent with prior observations that excessive fat accumulation in muscle can be attributed to obesity, diabetes, and aging. In addition, we identify two upregulated lncRNAs that warrant further investigation as candidate non-coding regulators of myogenesis. Moreover, the finding of two human MuSC subpopulations mirrors similar observations made from mouse muscle scRNA-seq analyses and agrees with the general conceptual framework that MuSCs transition between quiescent, activated, and cycling states. Future studies comparative analysis of these MuSC subpopulations across species may reveal human-specific aspects of myogenesis.

21Ligand-receptor interaction models from scRNA-seq data can help formulate new hypotheses about cell-communication channels that regulate muscle function. Identifying new MuSCs surface receptors will also help us refine MuSC purification protocols for prospective isolation studies used for in vitro and transplantation models. Our LR model revealed a set of 40 surface receptor genes that are distinctly expressed between MuSC1 and MuSC2, confirming some prior reports but also providing new candidate surface antigens for human MuSC subpopulation fractionation (Fig. 3). For example, we identify that SDC2 may mark “quiescent” MuSCs while CD44, TNFRSF12, and RPSA “early-activated” MuSCs in aging and disease contexts. In addition, our model proposed 79 cell-communication signals that may act between MuSCs and other cell types, in particular with fibroblasts, myofibers and immune cells through the EGFR receptor, and with vascular cells through the NOTCH3 receptor. These interactions may be critical regulators of muscle homeostasis and should be further investigated.

22This study presents a new set of candidate receptor expression signatures that may define human MuSC subpopulations (Fig. 3b) and provide human-specific receptor patterns (Fig. S2C). This approach is complimentary to receptor screening approaches, which have previously been useful to identify EGFR and CD82 as human MuSC receptor markers for flow cytometry. The subpopulation-specific receptor genes identified here may allow for further comparison of molecular and functional human MuSC diversity across muscle groups.

23Our study has some limitations. First, the sample size is small, and donors are very diverse, thus limiting our ability to control for age and sex. Therefore, we could not examine cell composition or gene expression trends based on muscle group, donor sex, or donor age. Even for samples from the same muscle (e.g., flexor hallucis longus [donors 2 and 7] or external oblique [donors 6 and 9]), we were unable to perform these comparions with statistical power. Further, we performed differential expression and gene set enrichment analyses within the MuSC1 and MuSC2 populations between the four middle-age (43–69 years old) and six aged (70–81 years old) donors, but found few age-cohort specific differences (data not shown). Second, future studies should aim at collecting muscle specimens in a more controlled manner, for example using a Bergström needle from a unique anatomical site; though this would not be possible for some muscles presented in this study. These biopsies would allow for aging and disease comparative analyses. Indeed, a recent report by Rubenstein et al. performed scRNA-seq on four human vastus lateralis muscle biopsies found that myofiber type composition and gene expression alterations based on donor age.

24Nevertheless, our dataset offers a new transcriptomic cell reference atlas and computational data integration approaches as a benchmark resource to examine human muscle cell diversity in health, aging, and disease.

25All procedures were approved by the Institutional Review Board at Weill Cornell Medical College (WCMC IRB Protocol # 1510016712) and were performed in accordance with relevant guidelines and regulations. All specimens were obtained at the New York-Presbyterian/Weill Cornell campus. All subjects provided written informed consent prior to participation. Samples were de-identified in accordance to IRB guidelines, and only details concerning age, sex, and anatomic origin were included. Sample anatomic locations and donor details are provided in Fig. 1a.

26After collection from donors during surgery, the muscle samples were cleared from excessive fat and connective tissue and weighted. About 50–65 mg of tissue was then digested into a single-cell suspension following a previously reported protocol. Briefly, the specimen was digested in 8 mg/mL Collagenase D (Roche) and 4.8 U/mL Dispase II (Roche) for 1 h followed by manual dissociation, filtration, and red blood cell lysis (Table 1). All single-cell suspensions were then frozen at -80 °C in 90% FBS, 10% DMSO and were re-filtered after thawing and prior to generating scRNA-seq libraries. The sequencing libraries were prepared using the Chromium Single Cell 3' reagent V2 or V3 kit (10X Genomics) in accordance with the manufacturer’s protocol and diluted as to yield a recovery of ~ 6000 single-cell transcriptomes with < 5% doublet rate (Table 1). The libraries were sequenced in multiplex (n = 2 per sequencing run) on the NextSeq 500 sequencer (Illumina) to produce between 200 and 250 million reads per library.

27Sequencing reads were processed with the Cell Ranger version 3.1 (10X Genomics) using the human reference transcriptome GRCh38. The downstream analysis was carried out with R 3.6.1 (2019-07-05). Quality control filtering, data clustering, visualization, and differential gene expression analysis was carried out using Seurat 3.1.0 R package. Each of the 10 datasets was first analyzed and annotated independently before integration with Scanorama (Table 1). Filtering retained cells with > 1000 unique molecular identifiers (UMIs), < 20% UMIs mapped to mitochondrial genes, and genes expressed in at least 3 cells (Fig. S4). Unsupervised shared nearest neighbor (SNN) clustering was performed with a resolution of 0.4 following which clusters were annotated with a common nomenclature of 12 cell type terms (Fig. S1). Differential expression analysis was achieved using either Seurat’s “FindAllMarkers” (Fig. 1d) or “FindMarkers” (Fig. 2a) function using a Wilcoxon Rank Sum test and only considering genes with > log2(0.25) fold-change and expressed in at least 25% of cells in the cluster. P values were corrected for false-discovery (FDR) and then reported as q values. Integration of raw counts was achieved using the “scanorama.correct” function from Scanorama. The integrated values were finally scaled in Seurat regressing out the 10X chemistry type and the number of genes per cell. Visualization was done using uniform manifold approximation and projection (UMAP). In Fig. S2, we integrated these human scRNA-seq datasets with a cohort of adult mouse muscle scRNA-seq datasets collected 0–7 days post-notexin injury. For multi-species integration, scRNA-seq datasets were integrated using first Scanorama and then Harmony to align related cell populations across species. Mouse genes were converted to human orthologs using biomaRt Bioconductor R package (Table 1). For differential expression analysis between human and mouse samples, we compared human MuSCs (combining MuSC1 + 2 clusters) and the uninjured mouse MuSCs to focus on cells from the homeostatic conditions.

28The list of differentially expressed genes between MuSC1 and MuSC2 (Fig. 2a) was used in Ingenuity Pathway Analysis (IPA) (QIAGEN, 2019-08-30). Activated (canonical) pathways were calculated by “Core Analysis” setting a q value cutoff of 0.05, which yielded 964 genes (366 down, 598 up). Top canonical pathways were chosen based of − log(p value) and z score values. Gene set enrichment analysis (GSEA, v.4.0.3) was ran on the same gene list as IPA ranked by log2 fold-change and with default program settings (Table 1). Gene sets database used the following: h.all.v7.0.symbols.gmt, c2.all.v7.0.symbols.gmt, c5.all.v7.0.symbols.gmt (Broad Institute). Gene sets enriched in phenotype were selected based on q value and enrichment score (ES).

29The model aims at scoring potential ligand-receptor interactions between MuSCs (receptor) and other cell types (ligand). We used the ligand-receptor interaction database from Ramilowski et al. (Table 1). From the database, we considered 1915 ligand-receptor pairs (from 542 receptors and 518 ligands) to test for differential expression in our scRNA-seq dataset. To calculate the score for a given ligand-receptor pair, we multiply the average receptor expression in MuSCs by the average ligand expression per other cell type. We only considered receptors that are differentially expressed in either the MuSC1 or MuSC2 subpopulation when compared individually to all other cell types.