Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

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

A reference single-cell transcriptomic atlas of human skeletal muscle tissue reveals bifurcated... Single-cell RNA-sequencing (scRNA-seq) facilitates the unbiased reconstruction of multicellular tissue systems in health and disease. Here, we present a curated scRNA-seq dataset of human muscle samples from 10 adult donors with diverse anatomical locations. We integrated ~ 22,000 single-cell transcriptomes using Scanorama to account for technical and biological variation and resolved 16 distinct populations of muscle-resident cells using unsupervised clustering of the data compendium. These cell populations included muscle stem/progenitor cells (MuSCs), which bifurcated into discrete “quiescent” and “early-activated” MuSC subpopulations. Differential expression analysis identified transcriptional profiles altered in the activated MuSCs including genes associated with aging, obesity, diabetes, and impaired muscle regeneration, as well as long non-coding RNAs previously undescribed in human myogenic cells. Further, we modeled ligand-receptor cell-communication interactions and observed enrichment of the TWEAK-FN14 pathway in activated MuSCs, a characteristic signature of muscle wasting diseases. In contrast, the quiescent MuSCs have enhanced expression of the EGFR receptor, a recognized human MuSC marker. This work provides a new benchmark reference resource to examine human muscle tissue heterogeneity and identify potential targets in MuSC diversity and dysregulation in disease contexts. Introduction β1-integrin (CD29), NCAM (CD56), EGFR, and CD82 to Skeletal muscles are essential to daily functions such as varying purities [6–10]. With aging, human MuSCs locomotion, respiration, and metabolism. Upon damage, exhibit a heterogeneous expression of the senescence Ink4a resident muscle stem cells (MuSCs) repair the tissue in marker p16 and accumulate other cell-intrinsic coordination with supporting non-myogenic cell types alterations in myogenic gene expression programs, cell such as immune cells, fibroblasts, and endothelial cells [1]. cycle control, and metabolic regulation [2, 11]. However, However, with age and disease, the repair capacity of given their varied molecular and functional states, our MuSCs declines, leading to complications such as fibrotic understanding of MuSCs in adult human muscle tissue scarring, reduced muscle mass and strength [2, 3], fat remains incompletely defined. In addition, cellular coordin- accumulation, and decreased insulin sensitivity [4], all of ation in the regulation of human muscle homeostasis and which severely affect mobility and quality of life [5]. regeneration remains poorly understood due to the lack of Human MuSCs are defined by the expression of the experimentally tractable models with multiple human paired box family transcription factor PAX7 and can be muscle cell types. Given these challenges, we posited that isolated using various surface marker proteins including an unbiased single-cell reference atlas of skeletal muscle could provide a useful framework to explore MuSC vari- * Correspondence: bdc68@cornell.edu ability and communication in adult humans. Meinig School of Biomedical Engineering, Cornell University, Ithaca, NY Here, we deeply profiled the transcriptome of thou- 14853, USA Full list of author information is available at the end of the article sands of individual MuSCs and muscle-resident cells © The Author(s). 2020 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Micheli et al. Skeletal Muscle (2020) 10:19 Page 2 of 13 from diverse adult human muscle samples using single- expression analysis between the 16 distinct subpopulations cell RNA-sequencing (scRNA-seq). After integrating identified an extensive set of unique markers that we these donor datasets to conserve biological information grouped into 4 categories (Fig. 1d). and overcome technical variation, we resolved two subpopulations of MuSCs with distinct gene expression scRNA-seq resolves the cellular diversity of human muscle signatures. Using differential gene expression analysis and novel markers and ligand-receptor interaction modeling, we extend the We annotated and interpreted the consensus cell atlas known repertoire of human MuSC gene expression (Fig. 1b, d) into cell type subpopulations as follows. We programs, suggesting new regulatory programs that may identify four types of stromal cells starting with adipocytes be associated with human MuSC activation, as well as found to be expressing apolipoprotein D (APOD)[15]), the features of human muscle aging and disease. brown fat tissue adipokine CXCL14 [16], GPX3,and GLUL. Among the 3 other subpopulations of fibroblast-like cells, Results Fibroblast 1 expresses high levels of collagen 1 (COL1A1), Collection and integration of a diverse human scRNA-seq SFRP4, SERPINE1,and CCL2; Fibroblast 2 expresses fibro- dataset nectin (FBN1), the microfibril-associated glycoprotein We used scRNA-seq to collect and annotate a single-cell MFAP5,and CD55 known to be expressed by synoviocytes transcriptomic dataset of diverse adult human muscle [17]; and Fibroblast 3 is mainly characterized by SMOC2 samples under homeostatic conditions. The muscle identified in tendon fibroblasts [18]. The Fibroblast 3 clus- samples were from surgically discarded tissue from n =10 ter is similar to the adipocytes cluster though exhibits lower donors (range 41 to 81 years old) undergoing reconstruct- expression levels and frequencies of the marker genes ive procedures and originating from a wide variety of APOD, CXCL12,and GLUL, and contain pre-adipocytes. anatomical sites in otherwise healthy patients (Fig. 1a). We also identify 5 types of vascular cells, including 3 Each sample was ~ 50 mg after removal of extraneous fat endothelial subpopulations, and a subpopulation of peri- and connective tissue. Muscle samples were enzymatically cytes and smooth muscle cells (SMCs). Pericytes and digested into single-cell suspensions and independently SMCs express the canonical markers RGS5 and MYH11. loaded into the 10X Chromium system. All together, we Endothelial 1 express E-selectin (SELE), IL6, ICAM1,and collected over 22,000 human muscle single-cell transcrip- VCAM1. These genes are upregulated at sites of inflamma- tomes (2206 ± 1961 cells per dataset) into a single data tion to facilitate immune cell recruitment, suggesting this compendium. Using unsupervised clustering, we resolved Endothelial 1 cell population may be involved in homeo- 16 types of cells of immune, vascular, and stromal origin, static muscle tissue remodeling [19, 20]. Endothelial 2 cells as well as two distinct subpopulations of MuSCs and some are distinguished by expressing high levels of claudin-5 myofiber myonuclei (Fig. 1b). (CLDN5), ICAM2, and the chemokine CXCL2.Endothelial Given important differences in anatomical site, donor 3 expresses high levels of the platelet-recruiting Von Will- health history, age, sex, and surgical procedures, the ebrand Factor (VWF) and caveolin-1 (CAV1), a protein muscle samples were highly heterogeneous in terms of known to regulate cholesterol metabolism, atherosclerosis cell-type diversity and underlying gene expression progression, and MuSC activation [21, 22]. Endothelial 3 profiles. Comparing the resulting scRNA-seq datasets is cells are enriched for expression of BTLN9, suggesting they therefore a challenge that we addressed using recently might represent a lymphatic endothelial phenotype [23]. developed bioinformatic integration methods [12–14]. We also noted two types of myeloid immune cells: first, Our goal was to assemble a unified dataset of human tissue-resident and anti-inflammatory macrophages which muscle tissue that faithfully conserved sources of bio- express CD74 and histocompatibility complex HLA proteins; logical variability such as donor, anatomical location, second, activated macrophages and monocytes that express and cell composition heterogeneity, while accounting for inflammatory markers such as S100A9 (calgranulin) and technical biases. We tested four different scRNA-seq LYZ (lysozyme). Moreover, S100A9 transcript abundance data integration methods (Fig. S1 and S3) and found that levels have been shown to be a feature in aging and chronic Scanorama [13] followed by scaling the output by regres- inflammation [24]. We also identified a pool of T/B lympho- sing against the library chemistry technical variable (“10X cytes and NK cells characterized by IL7R and NKG7,re- chemistry”) and the number of genes detected per single- spectively, as well as a small subset of HBA1 erythroblasts. cell best satisfied this goal. Detailed information on our Finally, we identified two subpopulations of MuSCs methodology is provided in Fig. S1.After integratingthe (henceforth called “MuSC1” and “MuSC2”). MuSC1 highly 10 datasets, we noted remarkable consistency amid cell expressed the canonical myogenic transcription factor PAX7 types across donors (Fig. 1c, e), owing to the robustness of [25], as well as chordin-like protein 2 (CHRDL2)and Delta- scRNA-seq technology, the bioinformatic method chosen, like non-canonical Notch ligand 1 (DLK1). CHRDL2 has and our sample preparation protocol. Differential gene been shown previously to be expressed in freshly isolated Micheli et al. Skeletal Muscle (2020) 10:19 Page 3 of 13 Fig. 1 (See legend on next page.) Micheli et al. Skeletal Muscle (2020) 10:19 Page 4 of 13 (See figure on previous page.) Fig. 1 Single-cell transcriptomic map of human muscle tissue biopsies. a Metadata (sex, age, anatomical site, and the number of single-cell transcriptomes after quality control (QC) filtering) from n = 10 donors. Colors indicate sample anatomical sites. b Scanorama-integrated and batch-effect corrected transcriptomic atlas revealing a consensus description of 16 distinct muscle-tissue cell populations. c Transcriptomic atlas colored by donor and anatomical location. d Dot-plot showing differentially expressed genes that distinguish the cell populations. Grouped in four compartments: muscle, endothelial/vascular, stromal, and immune. e Cell type proportions as annotated in (b) across the 10 donors and grouped by body sections. L, leg (donors 02, 07, 08); T, trunk (donors 01, 05, 06, 09, 10); F, face (donors 03, 04) quiescent human MuSCs [7], though its function is still to In addition, IL32 has been shown to have inflammatory be understood. DLK1 is an inhibitor of adipogenesis whose properties in human obesity [37] and have a negative impact role in muscle has mainly been recognized in the on insulin sensitivityand myogenesis [38], while embryo but remains controversial in adult muscle re- TNFRSF12/FN14 has been implicated in various muscle generation [26–28]. In contrast to MuSC1, MuSC2 wasting diseases [39, 40] and metabolic dysfunction [41]. expressed lower levels of PAX7 but maintain expression Furthermore, the MuSC2 population is enriched for riboso- of MYF5 (a marker of activated MuSCs) and APOC1 mal gene expression (e.g. RPLP1 and RPS6; data not shown), (Fig. 2b). Interestingly, the MuSC2 population also had indicating that these cells may have elevated translational elevated expression of two long non-coding RNAs mechanisms. Lastly, the MuSC1 population has enriched ex- (lncRNAs), LINC00152,and MIR4435-2HG.LncRNAs pression of the myogenic gene PAX7 and, to a lesser extent, are involved in regulating myogenesis [29]. Surprisingly, MYF5, compared the MuSC2 population. These observa- we detected low expression of the myogenic commitment tions suggest that MuSC1 is comprised of quiescent MuSCs, factors MYOD1 and MYOG (Fig. 2b), in contrast to scRNA- and MuSC2 is comprised of an early-activated MuSCs. seq analyses of adult mouse muscle [30, 31]. These observa- We performed Ingenuity Pathway Analysis (IPA) to com- tions suggest that the MuSC1 and MuSC2 populations are pare biological processes differentially activated between both comprised largely of muscle stem cells, not commit- the MuSC1 and MuSC2 populations. The IPA gene group ted myogenic progenitors. In addition, we noted that “Oxidative Phosphorylation” is enriched in MuSC1 [42], “Myonuclei” population (Fig. 1b) was enriched for myosin while “EIF2 Signaling,” associated with protein translation light chain (MYLFP), skeletal alpha-actin (ATCA1), and processes, is enriched in MuSC2 (Fig. 2c). Furthermore, troponin C (TNNC2), proteins involved in muscle con- Gene Set Enrichment Analysis (GSEA) also found MuSC1 traction. This multiple-donor scRNA-seq atlas highlights to be enriched for “myogenesis,”“muscle cell differenti- the cellular diversity of human muscle tissue and revealed ation,”“hypoxia,” and “response to mechanical stimulus” two distinct MuSC subpopulations along with specific gene sets, supporting the observation that these cells are myogenic expression programs. both less differentiated and may exhibit enhanced tran- scriptional responses to mechanical disruption due to tissue Homeostatic human muscle contains two distinct MuSC dissociation [32–34](Fig. 2d). MuSC2 cells are enriched for subpopulations “ribosome and translational initiation,”“MYC targets,” and We examined genes that were differentially expressed be- “E2F (cell proliferation),”“G2M checkpoint (cell division),” tween the MuSC1 and MuSC2 subpopulations and the bio- and “inflammation” gene sets, further supporting the inter- logical processes that characterize them (Fig. 2a, b). The pretation that these cells may be in an early activated or MuSC1 subpopulation was enriched for PAX7, DLK1,and partially differentiated state within an inflammatory envir- CHRDL2, as well as for the cyclin-dependent kinase inhibi- onment (Fig. 2d). Taken together, these observations KIP2 tor CKDN1C (encoding P57 ), suggesting that these cells suggest that the MuSC1 population is comprised of quies- are quiescent and not cycling. In addition, this subpopula- cent MuSCs, while the MuSC2 population is comprised of tion expresses the transcription factor BTG2,which was active, proliferating, and/or dysregulated MuSCs, with identified in mouse to be enriched in quiescent MuSCs [30]. expression alterations associated with inflammation, aging, We also note that the MuSC1 subpopulation expressed and muscle wasting. Differentially expressed genes such as elevated levels of mitochondrial genes as well as FOS, JUN, IL32, CXCL1, CCL2,and TNFRSF12/FN14 may constitute and ERG1. Upregulation of these genes has been shown to a marker set for MuSC variation in chronic muscle inflam- be potential artefacts of the enzymatic digestion during the mation in various pathologies. sample preparation [32–34]. The MuSC2 subpopulation was enriched for multiple Ligand-receptor interaction model identifies potential markers of inflammation including CCL2, CXCL1, IL32,and surface markers and cell-communication channels in surface receptor TNFRSF12/FN14. In particular, CCL2 and human skeletal muscle homeostasis CXCL1 are inflammatory cytokines known to be upregu- We used a ligand-receptor (LR) interaction model and a lated in muscle repair, exercise, and fat metabolism [35, 36]. database of LR pairs [43] to map cell signaling Micheli et al. Skeletal Muscle (2020) 10:19 Page 5 of 13 Fig. 2 Gene expression and pathway analysis comparison between two MuSC subpopulations. a Volcano plot from comparing transcript levels between all cells within the “MuSC1” and “MuSC2” subpopulations. Log fold-change in normalized gene expression versus −log adjusted p value 2 10 plotted. Differentially expressed genes (adjusted p value <0.05) are colored dark or light blue (based on their enrichment in MuSC1 or MuSC2, respectively). Genes with log fold-change > 0.75 are labeled. b Normalized expression values of select differentially expressed genes. q values reported in inset. c Top activated canonical pathways by Ingenuity Pathway Analysis (IPA) based on differentially expressed genes and ranked by p value. Pathways significantly enriched in either population with |z score| > 1 are indicated in blue. d Select gene ontology (GO) terms and hallmark pathways enriched between the MuSC subpopulations as identified by gene set enrichment analysis (GSEA) and ranked by enrichment score (ES) communication channels in human muscle and uncover expressed receptors on a given cell population (e.g., differences between MuSC1 and MuSC2 subpopula- “MuSC1”) relative to all other population and ligands tions (Fig. 3). The model also identifies interacting li- expressed by other cell types. The MuSC1 and MuSC2 gand(s) and is restricted to receptor genes differentially subpopulations are involved in numerous LR interac- expressed by a specific cell type within the consensus tions, as both ligand- and receptor-expressing cells (Fig. human muscle cell atlas (Fig. 1b). For each LR pair, the 3a), though a majority of all LR interaction pairs instead model calculates an interaction score from differentially involve other cell types. This suggests that only a small Micheli et al. Skeletal Muscle (2020) 10:19 Page 6 of 13 Fig. 3 (See legend on next page.) Micheli et al. Skeletal Muscle (2020) 10:19 Page 7 of 13 (See figure on previous page.) Fig. 3 Differentially expressed receptors and ligand-receptor interaction between cell populations. a Chord plot of all ligand-receptor (LR) interactions across cell populations/types within the consensus atlas based on co-expression. Each cell type is color-coded with its receptor genes more displaced from the perimeter than its ligand genes. All interactions not involving MuSC1 or MuSC2 are presented in gray. b List of differentially expressed genes between the MuSC1 and MuSC2 subpopulation ranked by log fold-change in expression. Positive average values correspond to genes that are upregulated in MuSC1, whereas negative values are upregulated in MuSC2. Receptors that are statistically significant (FDR-corrected q value < 0.05) are colored in blue. Receptors that are not statistically significant are in gray. c Heatmap representing row- normalized (Z-score) LR interaction scores. Rows represent ligand-receptor interaction pairs in the format LIGAND_RECEPTOR, where the receptor is either differentially expressed in the MuSC1 or MuSC2 populations compared to all the other cell types. Columns identify cell types expressing the ligand. Asterisks after the pair name also indicate that the ligand is differentially expressed by the other cell type and that interaction is likely cell-type specific. Red pairs involve the EGFR receptor, purple pairs the NOTCH3 receptor. A positive value indicates that the interaction has a high score for a particular ligand and cell type compared to other cell types subset of potential paracrine interactions in human (EREG). Other EGFR ligands include brevican (BCAN), and muscle may include MuSCs. betacellulin (BTC) produced by endothelial cells; ECM pro- Given the distinct expression profiles between the teins fibulin 3 (EFEMP1), decorin (DCN), and tenascin C MuSC1 and MuSC2 populations, we sought to identify (TNC) expressed by fibroblasts; and FGF13, AHM, NRG4, genes that could facilitate surface antigen-based separation and EGF, expressed by mature skeletal myofibers. We also of these two human MuSC populations for prospective iso- detect seven interactions involving NOTCH3 with a variety lation strategies. We identified surface receptor genes that of ligands. Notch3 signaling is involved in maintaining were differentially expressed between the MuSC1 and MuSC quiescence, in particular through interaction with MuSC2 populations, using a database of 542 human sur- DLL4 [50], whichwefound differentlyexpressedby endo- face “receptor” genes [43](Fig. 3). MuSC1 exhibit elevated thelial cells along with JAG2. In addition, NOTCH3 also in- expression of EGFR, ITGB1, FGFR4, SDC2,aswellasthe teracts with the ECM protein thrombospondin-2 (THBS2). three tetraspanins CD81, CD82,and CD151(Fig. 3b). EGFR Only two receptors, TNFRSF12/FN14 and RPSA,were is a recently established human MuSC marker and is re- found differentially expressed in MuSC2 compared to quired for basal-apical asymmetric cell division [7, 10]. The other cell types. The first, TNFRSF12/FN14, interacts with tetraspanin CD82 is also a recently recognized human the TWEAK cytokine ligand. While typically recognized to MuSC maker [6], while CD9 and CD81 have been identi- be expressed by macrophages and other immune cells fied to control muscle myoblast fusion [44]. Furthermore, [51], our model suggests that TWEAK is also expressed by Syndecans (SDCs) have been identified in mouse to be het- the Fibroblast 2 and pericyte cell populations, though not erogeneously expressed on MuSCs and myoblasts during in a statistically significant manner. The second, RPSA,is muscle repair [30] and have been shown to form co- surface ribosomal protein that interacts with laminins receptor complexes with integrin β1 (ITGB1) and FGFR4 (LAM), a dual-specificity phosphatase 18 (DUSP18), and upstream of signaling pathways regulating myogenesis prion protein 2 (PRND), which taken together may sug- [45]. Only SDC4 and SDC3 have yet been identified to gest various pathological processes such as prion diseases mark adult mouse MuSCs [46]. In comparison, the MuSC2 and cancer [52, 53]. Together, this ligand-receptor analysis subpopulation has elevated expression of CD44 and identifies a broad set of surface markers that could refine TNFRSF12/FN14 as previously noted (Fig. 3b). The CD44 the molecular definition of human MuSCs and their receptor has been shown to regulate myoblast migration subpopulations, as well as candidate cell-communication andfusioninmouse,but also mark MuSCsinosteoarthritis channels differentially involved in healthy and diseased patients [47, 48]. muscle tissues. Next, we focused the LR analysis on the MuSC1 and Lastly, we performed a comparative analysis of receptor MuSC2 populations. We identified 73 and 6 significant LR gene expression between mouse and human MuSCs. We interactions for the MuSC1 and MuSC2 populations, integrated the human scRNA-seq datasets described in respectively (Fig. 3c). Over one third of all interactions in Fig. 1 and an adult mouse muscle injury-response scRNA- the MuSC1 subpopulation involve the EGFR receptor, seq time-course previously reported [30]byconverting which has recently been shown to play a critical role in mouse genes to their corresponding human ortholog. The directing MuSC asymmetric division in regenerating multi-species scRNA-seq atlas was integrated with Sca- muscle [10]. A limited number of EGFR ligands have been norama and corrected with Harmony (Fig. S2A-B) [54]. identified in muscle repair, for example, amphiregulin. From this integrated atlas, we annotated all clusters as in (AREG) secreted by T cells [49]. According to our model Fig. 1. We identified two MuSC clusters which both con- reg findings, EGFR may also interact with ligands expressed by tained cells from both mouse and human samples. We immune cells, such as with TGF-α (TGFA), heparin-biding then performed differential expression analysis between EGF (HBEGF), amphiregulin (AREG), and epiregulin species comparing aggregated human MuSC1 and MuSC2 Micheli et al. Skeletal Muscle (2020) 10:19 Page 8 of 13 cells to mouse MuSCs from the uninjured timepoint (Fig. due to its differential expression of BTLN9, a marker of S2C). We found that EGFR and CD99 were most differen- lymphatic endothelial cells [23]. tially expressed by human MuSCs and, conversely, CRLF1 Most notably, this analysis suggests that human and SDC4 were most enriched in mouse MuSCs. This muscle may contain two distinct MuSC subpopulations findings suggest that mouse and human MuSC exhibit (Fig. 2). This finding contrasts with Rubenstein et al. species-specific receptor expression signatures. which observed a single MuSC (“satellite cell”) popula- tion from dissociated whole muscle samples and Barruet et al. which observed ~ 12 clusters from human MuSCs Discussion prospectively enriched by CXCR4+/CD29+/CD56+ Here we present an annotated multi-donor single-cell FACS. Since cluster distinction depends on both the cel- RNA sequencing dataset consisting of 22,000 single-cell lular diversity and sample complexity, it is expected that transcriptomes from 10 different donors and unique variation in study design and methods will yield differing anatomical sites, some of which difficult to access out- conclusions regarding sub-population resolution. In this side of reconstructive surgeries. Our study complements work, we found a “MuSC1” subpopulation to be largely other recent reports by Rubenstein et al. and Barruet comprised of “quiescent” MuSCs, owning to high levels et al., which collected dissociated whole vastus lateralis of PAX7, the mitotic inhibitor CDKN1C, and DLK1. muscles and FACS-sorted MuSC samples mostly from Interestingly, DLK1 may be an important regulator for vastus lateralis muscles, respectively, by providing more human MuSC maintenance and a marker of healthy tis- diversity in anatomical sites and donor demographics sue given its role in inhibiting adipogenesis [26]. Con- [55, 56]. As such, these scRNA-seq data exhibited versely, we identified in the “MuSC2” population notable biological and technical variation, and therefore, signatures of inflammation and increased fat metabolism we applied the bioinformatic method Scanorama to assem- (CCL2 and CXCL1), reduced insulin sensitivity (IL32), ble an integrated cellular atlas with minimal technical cell cycle (EIF2 Signaling terms), and muscle wasting biases so that we could examine the cellular heterogeneity (TNFRSF12/FN14), thereby suggesting that these cells across diverse adult human muscle tissue samples. We ob- may constitute an “early-activated” and possibly dysfunc- served that Scanorama performed more successfully than tional MuSC pool. These markers are consistent with other data integration approaches, especially when includ- prior observations that excessive fat accumulation in ing a scaling regression for sequencing chemistry (Fig. S1 muscle can be attributed to obesity, diabetes, and aging and S3). Notably, even after performing Scanorama with [4]. In addition, we identify two upregulated lncRNAs scaling, we still observed that integrated atlas exhibited bio- that warrant further investigation as candidate non- logical (donor) and technical (sequencing chemistry) coding regulators of myogenesis [29]. Moreover, the biases, but retained some degree of donor-specific cell-type finding of two human MuSC subpopulations mirrors subpopulations. similar observations made from mouse muscle scRNA-seq We describe the muscle tissue cellular heterogeneity analyses [30, 31] and agrees with the general conceptual and provide a comprehensive analysis of differentially framework that MuSCs transition between quiescent, acti- expressed genes for 16 resolved cell subpopulations (Fig. vated, and cycling states [1]. Future studies comparative 1), adding to a growing literature documenting human analysis of these MuSC subpopulations across species may muscle cell transcriptional diversity [55–57]. Compared reveal human-specific aspects of myogenesis. to other studies, the broader variety of muscle tissue Ligand-receptor interaction models from scRNA-seq samples combined with the lack of FACS selection data can help formulate new hypotheses about cell- allowed us to identify candidate subpopulations of communication channels that regulate muscle function muscle fibroblasts and vascular endothelial cells that [30]. Identifying new MuSCs surface receptors will also may provide unique perspective to human muscle physi- help us refine MuSC purification protocols for prospect- ology. In particular, we remark that Endothelial 1 ive isolation studies used for in vitro and transplantation expressed DARC/ACKR1, a gene identified in mouse and models. Our LR model revealed a set of 40 surface human [56, 58] to mark cells of post-capillary venular receptor genes that are distinctly expressed between origin (Fig. 1d). Rubenstein et al. also found a DARC/ MuSC1 and MuSC2, confirming some prior reports but ACKR1+ post-capillary venular endothelial cluster and a also providing new candidate surface antigens for human second VWF+ FABP+ cluster, which overlaps with the MuSC subpopulation fractionation (Fig. 3). For example, Endothelial 2 and 3 clusters reported here. We suggest we identify that SDC2 may mark “quiescent” MuSCs that the Endothelial 2 cluster may contain both arterial while CD44, TNFRSF12, and RPSA “early-activated” and capillary endothelial cells, but could not further par- MuSCs in aging and disease contexts. In addition, our tition and classify this cluster. We suggest that the Endo- model proposed 79 cell-communication signals that may thelial 3 cluster may represent lymphatic endothelium act between MuSCs and other cell types, in particular Micheli et al. Skeletal Muscle (2020) 10:19 Page 9 of 13 with fibroblasts, myofibers and immune cells through would not be possible for some muscles presented in this the EGFR receptor, and with vascular cells through the study. These biopsies would allow for aging and disease NOTCH3 receptor. These interactions may be critical comparative analyses. Indeed, a recent report by Ruben- regulators of muscle homeostasis and should be further stein et al. [56] performed scRNA-seq on four human vas- investigated. tus lateralis muscle biopsies found that myofiber type This study presents a new set of candidate receptor composition and gene expression alterations based on expression signatures that may define human MuSC sub- donor age. populations (Fig. 3b) and provide human-specific receptor Nevertheless, our dataset offers a new transcriptomic patterns (Fig. S2C). This approach is complimentary to re- cell reference atlas and computational data integration ap- ceptor screening approaches, which have previously been proaches as a benchmark resource to examine human useful to identify EGFR and CD82 as human MuSC muscle cell diversity in health, aging, and disease. receptor markers for flow cytometry [6, 7, 9]. The subpopulation-specific receptor genes identified here may Methods allow for further comparison of molecular and functional Human participation for muscle sample collection human MuSC diversity across muscle groups [59, 60]. All procedures were approved by the Institutional Review Our study has some limitations. First, the sample size is Board at Weill Cornell Medical College (WCMC IRB small, and donors are very diverse, thus limiting our ability Protocol # 1510016712) and were performed in accordance to control for age and sex. Therefore, we could not exam- with relevant guidelines and regulations. All specimens ine cell composition or gene expression trends based on were obtained at the New York-Presbyterian/Weill Cornell muscle group, donor sex, or donor age. Even for samples campus. All subjects provided written informed consent from the same muscle (e.g., flexor hallucis longus [donors prior to participation. Samples were de-identified in 2 and 7] or external oblique [donors 6 and 9]), we were un- accordance to IRB guidelines, and only details concerning able to perform these comparions with statistical power. age, sex, and anatomic origin were included. Sample ana- Further, we performed differential expression and gene set tomic locations and donor details are provided in Fig. 1a. enrichment analyses within the MuSC1 and MuSC2 popu- lations between the four middle-age (43–69 years old) and Muscle digestion and single-cell sequencing library six aged (70–81 years old) donors, but found few age- preparation cohort specific differences (data not shown). Second, future After collection from donors during surgery, the muscle studies should aim at collecting muscle specimens in a samples were cleared from excessive fat and connective more controlled manner, for example using a Bergström tissue and weighted. About 50–65 mg of tissue was then needle [61, 62] from a unique anatomical site; though this digested into a single-cell suspension following a previously Table 1 List of reagents and other resources used in this study Reagents Dispase II (neutral protease, grade II) Sigma-Aldrich 04942078001 Collagenase D, from Clostridium histolyticum Sigma-Aldrich 11088866001 Commercial kits Chromium Single Cell 3' Library & Gel Bead Kit v2 10X Genomics CG00052 (protocol) Chromium Single Cell 3' Library & Gel Bead Kit v3 10X Genomics CG000183 (protocol) Deposited data Human ligand-receptor database [43] https://www.ncbi.nlm.nih.gov/pubmed/26198319 Human scRNA-seq dataset This paper GSE143704 Mouse scRNA-seq dataset [29] GSE143437 Software packages and algorithms Cell Ranger 3.1.0 (July 24, 2019) 10X Genomics https://support.10xgenomics.com/single-cell-gene-expression/software/ downloads/latest Seurat 3.1.0 [14] https://github.com/satijalab/seurat Scanorama (online version as of 2019-11-19) [13] https://github.com/brianhie/scanorama Harmony (online version as of 2019-11-19) [54] https://github.com/immunogenomics/harmony biomaRt 2.43.1 (online version as of 2020-01-08) [64] https://bioconductor.org/packages/release/bioc/html/biomaRt.html Gene Set Enrichment Analysis (4.0.3) [65] http://software.broadinstitute.org/gsea/index.jsp Ingenuity Pathway Analysis (IPA, 2019-08-30) QIAGEN https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis/ Micheli et al. Skeletal Muscle (2020) 10:19 Page 10 of 13 reported protocol [63]. Briefly, the specimen was digested andthe uninjuredmouse MuSCsto focus on cellsfrom in 8 mg/mL Collagenase D (Roche) and 4.8 U/mL Dispase the homeostatic conditions. II (Roche) for 1 h followed by manual dissociation, filtra- tion, and red blood cell lysis (Table 1). All single-cell sus- Pathway and gene set enrichment analysis pensions were then frozen at -80 °C in 90% FBS, 10% The list of differentially expressed genes between DMSO and were re-filtered after thawing and prior to gen- MuSC1 and MuSC2 (Fig. 2a) was used in Ingenuity erating scRNA-seq libraries. The sequencing libraries were Pathway Analysis (IPA) (QIAGEN, 2019-08-30). Acti- prepared using the Chromium Single Cell 3' reagent V2 or vated (canonical) pathways were calculated by “Core V3 kit (10X Genomics) in accordance with the manufac- Analysis” setting a q value cutoff of 0.05, which yielded turer’s protocol and diluted as to yield a recovery of ~ 6000 964 genes (366 down, 598 up). Top canonical pathways single-cell transcriptomes with < 5% doublet rate (Table 1). were chosen based of − log(p value) and z score values. The libraries were sequenced in multiplex (n =2perse- Gene setenrichmentanalysis (GSEA, v.4.0.3)[65]was quencing run) on the NextSeq 500 sequencer (Illumina) to ran on the same gene list as IPA ranked by log fold- produce between 200 and 250 million reads per library. change and with default program settings (Table 1). Gene sets database used the following: h.all.v7.0.sym- bols.gmt, c2.all.v7.0.symbols.gmt, c5.all.v7.0.sym- Single-cell data analysis bols.gmt (Broad Institute). Gene sets enriched in Sequencing reads were processed with the Cell Ranger phenotype were selected based on q value and enrich- version 3.1 (10X Genomics) using the human reference ment score (ES). transcriptome GRCh38. The downstream analysis was carried out with R 3.6.1 (2019-07-05). Quality control Ligand-receptor cell communication model filtering, data clustering, visualization, and differential The model aims at scoring potential ligand-receptor gene expression analysis was carried out using Seurat interactions between MuSCs (receptor) and other cell 3.1.0 R package [14]. Each of the 10 datasets was first types (ligand). We used the ligand-receptor interaction analyzed and annotated independently before integration database from Ramilowski et al. [43] (Table 1). From the with Scanorama [13] (Table 1). Filtering retained cells database, we considered 1915 ligand-receptor pairs with > 1000 unique molecular identifiers (UMIs), < 20% (from 542 receptors and 518 ligands) to test for differen- UMIs mapped to mitochondrial genes, and genes tial expression in our scRNA-seq dataset. To calculate expressed in at least 3 cells (Fig. S4). Unsupervised the score for a given ligand-receptor pair, we multiply shared nearest neighbor (SNN) clustering was performed the average receptor expression in MuSCs by the aver- with a resolution of 0.4 following which clusters were age ligand expression per other cell type. We only con- annotated with a common nomenclature of 12 cell type sidered receptors that are differentially expressed in terms (Fig. S1). Differential expression analysis was either the MuSC1 or MuSC2 subpopulation when com- achieved using either Seurat’s “FindAllMarkers” (Fig. 1d) pared individually to all other cell types. or “FindMarkers” (Fig. 2a) function using a Wilcoxon Rank Sum test and only considering genes with > Supplementary information log (0.25) fold-change and expressed in at least 25% of Supplementary information accompanies this paper at https://doi.org/10. cells in the cluster. P values were corrected for false- 1186/s13395-020-00236-3. discovery (FDR) and then reported as q values. Integration of raw counts was achieved using the “scanorama.correct” Additional file 1: Figure S1. Comparison of scRNA-seq integration and batch correction methods. We compared four scRNA-seq data in- function from Scanorama. The integrated values were fi- tegration methods to evaluate which most faithfully conserves donor, nally scaled in Seurat regressing out the 10X chemistry anatomical, and biological information while minimizes technical biases. type and the number of genes per cell. Visualization was (A) The n = 10 donor datasets were first annotated independently using a nomenclature of 12 common cell type terms following unsupervised done using uniform manifold approximation and projec- SNN clustering. Then we evaluated the integration method by UMAP and tion (UMAP) [66]. In Fig. S2, we integrated these human by coloring the data either by cell type, donor ID, or 10X library chemistry scRNA-seq datasets with a cohort of adult mouse muscle used. First, we integrated the data by merging the individually normalized gene expression matrices without any further correction. We saw strong scRNA-seq datasets collected 0–7 days post-notexin injury technical biases that overwhelmed biological information as the different [30]. For multi-species integration, scRNA-seq datasets cell populations segregate by sample/donor and chemistry type. For in- were integrated using first Scanorama and then Harmony stance, the two MuSC and progenitor subpopulations are grouped with fibroblasts and endothelial cells. Second, we tested the Seurat SCT inte- [54] to align related cell populations across species. Mouse gration method [14] . This method first calculates a cross-correlation sub- genes were converted to human orthologs using biomaRt space from genes that are shared between datasets. We noticed that this Bioconductor R package [64](Table 1). For differential ex- method better “aligns” donor and chemistry type but at the expense of masking biological variability. For instance, we observed that the two pression analysis between human and mouse samples, we MuSC and four stromal subpopulations (Fibroblast 1,2,3 and Adipocytes) compared human MuSCs (combining MuSC1 + 2 clusters) Micheli et al. Skeletal Muscle (2020) 10:19 Page 11 of 13 Availability of data and materials were grouped together, hiding important biological heterogeneity. Al- The human muscle scRNA-seq datasets supporting the conclusions of this though certainly useful to validate reproducibility in scRNA-seq experi- article are archived at the NIH GEO repository under accession number ments, the Seurat SCT integration approach overcorrected biological GSE143704. heterogeneity for heterogeneous samples. Third, we tested the Scanor- ama method [13], which relies on a computer vision algorithm that Ethics approval and consent to participate “stitches” datasets together even when the cell type composition be- All procedures were approved by the Institutional Review Board at Weill tween dataset is considerably different. We see that this method groups Cornell Medical College (WCMC IRB Protocol # 1510016712) and were similar cell populations together while acknowledging donor differences. performed in accordance with relevant guidelines and regulations. All Yet, surprisingly, this method is also very sensitive at picking up differ- specimens were obtained at the New York-Presbyterian/Weill Cornell ences in chemistry. To correct this chemistry effect, we scaled the Scanor- campus. All subjects provided written informed consent prior to ama output by regressing out the chemistry and the number of genes participation. Samples were de-identified in accordance to IRB guidelines detected per cell (significantly different between chemistry type) (B). and only details concerning age, sex, and anatomic origin were included. Using this integration method, we observed clear separation of the inde- pendently annotated cell populations. We present the resulting Consent for publication Scanorama-integrated dataset as a “consensus atlas” (see Fig. 1b-c) of hu- man muscle that describes donor-to-donor differences while grouping Not applicable. cells that are similar together and removing technical biases. Figure S2. Integration of human and mouse scRNA-seq data sets allows com- Competing interests parison of MuSC receptor gene expression across species. We gen- The authors declare no conflicts of interest. erated an integrated scRNA-seq atlas including human sample datasets from Fig. 1 and an adult mouse muscle regeneration time-course from Author details De Micheli et al. [29]. These datasets were integrated using first Scanor- Meinig School of Biomedical Engineering, Cornell University, Ithaca, NY ama and then Harmony for alignment across species. (A) Multi-species in- 14853, USA. Englander Institute for Precision Medicine, Weill Cornell tegrated atlas presented by UMAP plot a colored by sample type. (B) Medicine, New York, NY 10021, USA. Division of Plastic Surgery, Weill Cornell Multi-species integrated atlas presented by UMAP plot and annotated by Medical College, New York, NY 10021, USA. cell-type clusters. (C) The human MuSC1 and MuSC2 clusters were grouped into a cumulative human MuSC cell population, which was Received: 30 January 2020 Accepted: 10 June 2020 compared to mouse MuSCs from the uninjured samples only. Receptor genes were analyzed between the mouse and human MuSC cells for dif- ferential expression. Differentially expressed genes with an FDR-corrected References q-value < 0.05 are shown in (C). Figure S3. Composition of single-cell 1. Bentzinger CF, Wang YX, Dumont NA, Rudnicki MA. Cellular dynamics in the reference atlas as a whole and in cell-type clusters by donor. (A) muscle satellite cell niche. EMBO Rep. 2013;14:1062–72. Visualization of donor (n = 10) contributions to the whole single-cell ref- 2. Blau HM, Cosgrove BD, Ho ATV. The central role of muscle stem cells in erence atlas. In each panel, the full atlas is presented as a UMAP plot, regenerative failure with aging. Nat Med. 2015;21:854. with the cells for an individual donor are colored and overlaid on cells 3. Järvinen TA, Järvinen M, Kalimo H. Regeneration of injured skeletal muscle from all other donors (in gray). Note the total number of cells assayed dif- after the injury. Muscles Ligaments Tendons J. 2014;3:337–45. fers for each donor (see Fig. 1a). (B) Bar plot representing the relative 4. Addison O, Marcus RL, LaStayo PC, Ryan AS. Intermuscular fat: a review of contribution of cells with each cell type cluster from each donor. Note the consequences and causes. Int J Endocrinol. 2014;2014:1–11. that the MuSC1 and MuSC2 clusters are also plotted as a combined clus- 5. Larsson L, Degens H, Li M, Salviati L, Lee YI, Thompson W, Kirkland JL, Sandri ter on the left side of the bar plot for reference. Figure S4. Transcrip- M. Sarcopenia: aging-related loss of muscle mass and function. Physiol Rev. tomic detection variation within human muscle reference atlas. 2018;99:427–511. UMAP plots featuring (left) the number of unique molecular identifiers 6. Alexander MS, Rozkalne A, Colletta A, Spinazzola JM, Johnson S, Rahimov F, (UMIs) and (right) number of genes detected per single cell. Note that QC Meng H, Lawlor MW, Estrella E, Kunkel LM, et al. CD82 is a marker for filtering removed all cells with less than 1000 UMIs (see Methods). prospective isolation of human muscle satellite cells and is linked to muscular dystrophies. Cell Stem Cell. 2016;19:800–7. 7. Charville GW, Cheung TH, Yoo B, Santos PJ, Lee GK, Shrager JB, Rando TA. Ex vivo expansion and in vivo self-renewal of human muscle stem cells. Stem Cell Reports. 2015;5:621–32. Acknowledgements 8. Pisani DF, Clement N, Loubat A, Plaisant M, Sacconi S, Kurzenne J-Y, The authors acknowledge helpful advice from colleagues in the Cosgrove Desnuelle C, Dani C, Dechesne CA. Hierarchization of myogenic and and Elemento groups, as well as Christopher Mendias at the Hospital for adipogenic progenitors within human skeletal muscle. Stem Cells. 2010;28: Special Surgery and Peter Schweitzer of Genomics Facility at the Cornell 2182–94. University Biotechnology Resource Center. Lastly, the authors are grateful for 9. Uezumi A, Nakatani M, Ikemoto-Uezumi M, Yamamoto N, Morita M, the human tissue donors. Yamaguchi A, Yamada H, Kasai T, Masuda S, Narita A, et al. Cell-surface protein profiling identifies distinctive markers of progenitor cells in human skeletal muscle. Stem Cell Reports. 2016;7:263–78. Authors’ contributions 10. Wang YX, Feige P, Brun CE, Hekmatnejad B, Dumont NA, Renaud J-M, A.J.D. and B.D.C. designed the study and wrote the manuscript. J.A.S. Faulkes S, Guindon DE, Rudnicki MA. EGFR-Aurka signaling rescues polarity obtained the human tissue samples. A.J.D. performed the tissue dissociations, and regeneration defects in dystrophin-deficient muscle stem cells by scRNA-seq, and data analysis, with supervision and assistance from B.D.C. and increasing asymmetric divisions. Cell Stem Cell. 2019;24:419–432.e6. O.E. All authors reviewed, read, and approved the final manuscript. 11. Sousa-Victor P, Gutarra S, García-Prat L, Rodriguez-Ubreva J, Ortet L, Ruiz-Bonilla V, Jardí M, Ballestar E, González S, Serrano AL, et al. Geriatric muscle stem cells switch reversible quiescence into senescence. Nature. 2014;506:316–21. Funding 12. Stuart T, Satija R. Integrative single-cell analysis. Nat Rev Genet. 2019;20:257–72. This work was financially supported by the National Institutes of Health 13. Hie B, Bryson B, Berger B. Efficient integration of heterogeneous single-cell under award R01AG058630 (to B.D.C.), a Glenn Medical Research Foundation transcriptomes using Scanorama. Nat Biotechnol. 2019;37:685–91. and American Federation for Aging Research Grant for Junior Faculty (to 14. Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM, Hao Y, B.D.C.), and a US Department of Education Graduate Assistantship in Areas of Stoeckius M, Smibert P, Satija R. Comprehensive integration of single-cell National Need under Award P200A150273 (to A.J.D.). The content is solely data. Cell. 2019;177:1888–1902.e21. the responsibility of the authors and does not necessarily represent the 15. Muffat J, Walker DW. Apolipoprotein D: an overview of its role in aging and official views of any of these funding sources. age-related diseases. Cell Cycle. 2010;9:269–73. Micheli et al. Skeletal Muscle (2020) 10:19 Page 12 of 13 16. Cereijo R, Gavaldà-Navarro A, Cairó M, Quesada-López T, Villarroya J, Morón- 37. Catalán V, Gómez-Ambrosi J, Rodríguez A, Ramírez B, Ortega VA, Ros S, Sánchez-Infantes D, Peyrou M, Iglesias R, Mampel T, et al. CXCL14, a Hernández-Lizoain JL, Baixauli J, Becerril S, Rotellar F, Valentí V, et al. IL-32α- brown adipokine that mediates brown-fat-to-macrophage communication induced inflammation constitutes a link between obesity and colon cancer. in thermogenic adaptation. Cell Metab. 2018, 28:750–763.e6. Oncoimmunology. 2017;6:e1328338. 17. Karpus ON, Kiener HP, Niederreiter B, Yilmaz-Elis AS, van der Kaa J, Ramaglia 38. Davegårdh C, Broholm C, Perfilyev A, Henriksen T, García-Calzón S, Peijs L, V, Arens R, Smolen JS, Botto M, Tak PP, et al. CD55 deposited on synovial Hansen NS, Volkov P, Kjøbsted R, Wojtaszewski JFP, et al. Abnormal collagen fibers protects from immune complex-mediated arthritis. Arthritis epigenetic changes during differentiation of human skeletal muscle stem Research & Therapy. 2015;17:6. cells from obese subjects. BMC Med. 2017;15:39. 39. Enwere EK, Lacasse EC, Adam NJ, Korneluk RG. Role of the TWEAK-Fn14- 18. De Micheli AJ, Swanson JB, Disser NP, Martinez LM, Walker NR, Oliver DJ, Cosgrove BD, Mendias CL. Single-cell transcriptomics identify extensive cIAP1-NF-κB signaling axis in the regulation of myogenesis and muscle heterogeneity in the cellular composition of mouse Achilles tendons. homeostasis. Front Immunol. 2014;5:34. BioRxiv. 2020b;801266. 40. Mittal A, Kumar A, Lach-Trifilieff E, Wauters S, Li H, Makonchuk D, Glass D, 19. Goncharov NV, Nadeev AD, Jenkins RO, Avdonin PV. Markers and Kumar A. The TWEAK-Fn14 system is a critical regulator of denervation- biomarkers of endothelium: when something is rotten in the state. induced skeletal muscle atrophy in mice. J Cell Biol. 2010;188:833–49. Oxidative Med Cell Longev. 2017;2017:9759735. 41. Sato S, Ogura Y, Kumar A. TWEAK/Fn14 signaling axis mediates skeletal 20. Watson C, Whittaker S, Smith N, Vora AJ, Dumonde DC, Brown KA. IL-6 acts muscle atrophy and metabolic dysfunction. Front Immunol. 2014;5:18. on endothelial cells to preferentially increase their adherence for 42. Ryall JG, Dell’Orso S, Derfoul A, Juan A, Zare H, Feng X, Clermont D, Koulnis lymphocytes. Clin Exp Immunol. 1996;105(1):112–9. M, Gutierrez-Cruz G, Fulco M, et al. The NAD + -dependent SIRT1 21. Fernández-Hernando C, Yu J, Dávalos A, Prendergast J, Sessa WC. deacetylase translates a metabolic switch into regulatory epigenetics in Endothelial-specific overexpression of caveolin-1 accelerates atherosclerosis skeletal muscle stem cells. Cell Stem Cell. 2015;16:171–83. in apolipoprotein E-deficient mice. Am J Pathol. 2010;177:998–1003. 43. Ramilowski JA, Goldberg T, Harshbarger J, Kloppmann E, Lizio M, 22. Volonte D, Liu Y, Galbiati F. The modulation of caveolin-1 expression controls Satagopam VP, Itoh M, Kawaji H, Carninci P, Rost B, et al. A draft network of satellite cell activation during muscle repair. FASEB J. 2004;19:237–9. ligand–receptor-mediated multicellular signalling in human. Nat Commun. 23. Fujimoto N, He Y, D’Addio M, Tacconi C, Detmar M, Dieterich LC. Single-cell 2015;6:7866. mapping reveals new markers and functions of lymphatic endothelial cells 44. Charrin S, Latil M, Soave S, Polesskaya A, Chrétien F, Boucheix C, Rubinstein in lymph nodes. PLoS Biol. 2020;18:e3000704. E. Normal muscle regeneration requires tight control of muscle cell fusion 24. Swindell WR, Johnston A, Xing X, Little A, Robichaud P, Voorhees JJ, Fisher by tetraspanins CD9 and CD81. Nat Commun. 2013;4:1674. G, Gudjonsson JE. Robust shifts in S100a9 expression with aging: a novel 45. Pawlikowski B, Vogler TO, Gadek K, Olwin BB. Regulation of skeletal muscle mechanism for chronic inflammation. Sci Rep. 2013;3:1215. stem cells by fibroblast growth factors. Dev Dyn. 2017;246:359–67. 46. Pisconti A, Bernet JD, Olwin BB. Syndecans in skeletal muscle development, 25. Kuang S, Chargé SB, Seale P, Huh M, Rudnicki MA. Distinct roles for Pax7 regeneration and homeostasis. Muscles Ligaments Tendons J. 2012;2:1–9. and Pax3 in adult regenerative myogenesis. J Cell Biol. 2006;172:103. 26. Andersen DC, Laborda J, Baladron V, Kassem M, Sheikh SP, Jensen CH. Dual 47. Mylona E, Jones KA, Mills ST, Pavlath GK. CD44 regulates myoblast migration role of delta-like 1 homolog (DLK1) in skeletal muscle development and and differentiation. J Cell Physiol. 2006;209:314–21. adult muscle regeneration. Development. 2013;140:3743. 48. Scimeca M, Bonanno E, Piccirilli E, Baldi J, Mauriello A, Orlandi A, Tancredi V, 27. Waddell JN, Zhang P, Wen Y, Gupta SK, Yevtodiyenko A, Schmidt JV, Bidwell Gasbarra E, Tarantino U. Satellite cells CD44 positive drive muscle CA, Kumar A, Kuang S. Dlk1 is necessary for proper skeletal muscle regeneration in osteoarthritis patients. Stem Cells Int. 2015;2015:469459. development and regeneration. PLoS One. 2010;5:e15055. 49. Burzyn D, Kuswanto W, Kolodin D, Shadrach JL, Cerletti M, Jang Y, Sefik E, Tan TG, Wagers AJ, Benoist C, et al. A special population of regulatory T 28. Zhang L, Uezumi A, Kaji T, Tsujikawa K, Andersen DC, Jensen CH, Fukada S. cells potentiates muscle repair. Cell. 2013;155:1282–95. Expression and Functional Analyses of Dlk1 in Muscle stem cells and 50. Low S, Barnes JL, Zammit PS, Beauchamp JR. Delta-like 4 activates Notch 3 mesenchymal progenitors during muscle regeneration. Int J Mol Sci. 2019; 20:3269. to regulate self-renewal in skeletal muscle stem cells. Stem Cells. 2018;36: 29. Hagan M, Zhou M, Ashraf M, Kim I-M, Su H, Weintraub NL, Tang Y. Long 458–66. noncoding RNAs and their roles in skeletal muscle fate determination. 51. Tajrishi MM, Zheng TS, Burkly LC, Kumar A. The TWEAK-Fn14 pathway: a Noncoding RNA Investig. 2017;1:24. potent regulator of skeletal muscle biology in health and disease. Cytokine 30. De Micheli AJ, Laurilliard EJ, Heinke CL, Ravichandran H, Fraczek P, Soueid- Growth Factor Rev. 2014;25:215–25. Baumgarten S, De Vlaminck I, Elemento O, Cosgrove BD. Single-cell analysis of 52. Pampeno C, Derkatch IL, Meruelo D. Interaction of human laminin receptor the muscle stem cell hierarchy identifies heterotypic communication signals with Sup35, the [PSI+] prion-forming protein from S. cerevisiae: a yeast involved in skeletal muscle regeneration. Cell Rep. 2020;30:3583–3595.e5. model for studies of LamR interactions with amyloidogenic proteins. PLoS 31. Dell’Orso, S., Juan, A.H., Ko, K.-D., Naz, F., Gutierrez-Cruz, G., Feng, X., and One. 2014;9:e86013. Sartorelli, V. (2019). Single-cell analysis of adult skeletal muscle stem cells in 53. Wu Y, Tan X, Liu P, Yang Y, Huang Y, Liu X, Meng X, Yu B, Wu M, Jin H. ITGA6 homeostatic and regenerative conditions. Development dev.174177. and RPSA synergistically promote pancreatic cancer invasion and metastasis via PI3K and MAPK signaling pathways. Exp Cell Res. 2019;379:30–47. 32. Machado L, Esteves de Lima J, Fabre O, Proux C, Legendre R, Szegedi A, Varet H, Ingerslev LR, Barrès R, Relaix F, et al. In situ fixation redefines 54. Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, Baglaenko Y, quiescence and early activation of skeletal muscle stem cells. Cell Rep. 2017; Brenner M, Loh P, Raychaudhuri S. Fast, sensitive and accurate integration 21:1982–93. of single-cell data with harmony. Nat Methods. 2019;16:1289–96. 33. van den Brink SC, Sage F, Vértesy Á, Spanjaard B, Peterson-Maduro J, Baron 55. Barruet E, Garcia SM, Striedinger K, Wu J, Lee S, Byrnes L, Wong A, Xuefeng CS, Robin C, van Oudenaarden A. Single-cell sequencing reveals S, Tamaki S, Brack AS, Pomerantz JH. Functionally heterogeneous human dissociation-induced gene expression in tissue subpopulations. Nat satellite cells identified by single cell RNA sequencing. eLife. 2020;9:e51576. Methods. 2017;14:935–6. 56. Rubenstein AB, Smith GR, Raue U, Begue G, Minchev K, Ruf-Zamojski F, Nair 34. van Velthoven CTJ, de Morree A, Egner IM, Brett JO, Rando TA. VD, Wang X, Zhou L, Zaslavsky E, Trappe TA, Sealfon SC. Single-cell Transcriptional profiling of quiescent muscle stem cells in vivo. Cell Rep. transcriptional profiles of human skeletal muscle. Sci Rep. 2020;10:229. 2017;21:P1994–2004. 57. Riddle ES, Bender EL, Thalacker-Mercer AE. Transcript profile distinguishes variability in human myogenic progenitor cell expansion capacity. Physiol 35. Harmon, B.T., Orkunoglu-Suer, E.F., Adham, K., Larkin, J.S., Gordish-Dressman, Genomics. 2018;50:817–27. H., Clarkson, P.M., Thompson, P.D., Angelopoulos, T.J., Gordon, P.M., Moyna, N.M., et al. (2010). CCL2 and CCR2 variants are associated with skeletal 58. Thiriot A, Perdomo C, Cheng G, Novitzky-Basso I, McArdle S, Kishimoto JK, muscle strength and change in strength with resistance training. J Appl Barreiro O, Mazo I, Triboulet R, Ley K, et al. Differential DARC/ACKR1 Physiol (1985) 109, 1779–1785. expression distinguishes venular from non-venular endothelial cells in 36. Pedersen L, Olsen CH, Pedersen BK, Hojman P. Muscle-derived expression of murine tissues. BMC Biol. 2017;15:45. the chemokine CXCL1 attenuates diet-induced obesity and improves fatty 59. Garcia SM, Tamaki S, Lee S, Wong A, Jose A, Dreux J, Kouklis G, Sbitany H, Seth acid oxidation in the muscle. American Journal of Physiology-Endocrinology R, Knott PD, et al. High-yield purification, preservation, and serial and Metabolism. 2012;302:E831–40. Transplantation of Human Satellite Cells. Stem Cell Reports. 2018;10:1160–74. Micheli et al. Skeletal Muscle (2020) 10:19 Page 13 of 13 60. Xu X, Wilschut KJ, Kouklis G, Tian H, Hesse R, Garland C, Sbitany H, Hansen S, Seth R, Knott PD, Hoffman WY, Pomerantz JH. Human satellite cell transplantation and regeneration from diverse skeletal muscles. Stem Cell Reports. 2015;5:419–34. 61. Sarver DC, Sugg KB, Disser NP, Enselman ERS, Awan TM, Mendias CL. Local cryotherapy minimally impacts the metabolome and transcriptome of human skeletal muscle. Sci Rep. 2017;7. 62. Tarnopolsky MA, Pearce E, Smith K, Lach B. Suction-modified Bergström muscle biopsy technique: Experience with 13,500 procedures. Muscle Nerve. 2011;43:716–25. 63. Spinazzola JM, Gussoni E. Isolation of primary human skeletal muscle cells. Bio-Protocol. 2017;7:e2591. 64. Durinck S, Spellman P, Birney E, Huber W. Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt. Nat Protoc. 2009;4:1184–91. 65. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102:15545. 66. Becht E, McInnes L, Healy J, Dutertre CA, Kwok IWH, Ng LG, Ginhoux F, Newell EW. Dimensionality reduction for visualizing single-cell data using UMAP. Nat Biotechnol. 2018. Publisher’sNote Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Skeletal Muscle Springer Journals

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

Download PDF
13 pages

Loading next page...
 
/lp/springer-journals/a-reference-single-cell-transcriptomic-atlas-of-human-skeletal-muscle-CNeJkrl0Io

References (60)

Publisher
Springer Journals
Copyright
Copyright © The Author(s) 2020
eISSN
2044-5040
DOI
10.1186/s13395-020-00236-3
Publisher site
See Article on Publisher Site

Abstract

Single-cell RNA-sequencing (scRNA-seq) facilitates the unbiased reconstruction of multicellular tissue systems in health and disease. Here, we present a curated scRNA-seq dataset of human muscle samples from 10 adult donors with diverse anatomical locations. We integrated ~ 22,000 single-cell transcriptomes using Scanorama to account for technical and biological variation and resolved 16 distinct populations of muscle-resident cells using unsupervised clustering of the data compendium. These cell populations included muscle stem/progenitor cells (MuSCs), which bifurcated into discrete “quiescent” and “early-activated” MuSC subpopulations. Differential expression analysis identified transcriptional profiles altered in the activated MuSCs including genes associated with aging, obesity, diabetes, and impaired muscle regeneration, as well as long non-coding RNAs previously undescribed in human myogenic cells. Further, we modeled ligand-receptor cell-communication interactions and observed enrichment of the TWEAK-FN14 pathway in activated MuSCs, a characteristic signature of muscle wasting diseases. In contrast, the quiescent MuSCs have enhanced expression of the EGFR receptor, a recognized human MuSC marker. This work provides a new benchmark reference resource to examine human muscle tissue heterogeneity and identify potential targets in MuSC diversity and dysregulation in disease contexts. Introduction β1-integrin (CD29), NCAM (CD56), EGFR, and CD82 to Skeletal muscles are essential to daily functions such as varying purities [6–10]. With aging, human MuSCs locomotion, respiration, and metabolism. Upon damage, exhibit a heterogeneous expression of the senescence Ink4a resident muscle stem cells (MuSCs) repair the tissue in marker p16 and accumulate other cell-intrinsic coordination with supporting non-myogenic cell types alterations in myogenic gene expression programs, cell such as immune cells, fibroblasts, and endothelial cells [1]. cycle control, and metabolic regulation [2, 11]. However, However, with age and disease, the repair capacity of given their varied molecular and functional states, our MuSCs declines, leading to complications such as fibrotic understanding of MuSCs in adult human muscle tissue scarring, reduced muscle mass and strength [2, 3], fat remains incompletely defined. In addition, cellular coordin- accumulation, and decreased insulin sensitivity [4], all of ation in the regulation of human muscle homeostasis and which severely affect mobility and quality of life [5]. regeneration remains poorly understood due to the lack of Human MuSCs are defined by the expression of the experimentally tractable models with multiple human paired box family transcription factor PAX7 and can be muscle cell types. Given these challenges, we posited that isolated using various surface marker proteins including an unbiased single-cell reference atlas of skeletal muscle could provide a useful framework to explore MuSC vari- * Correspondence: bdc68@cornell.edu ability and communication in adult humans. Meinig School of Biomedical Engineering, Cornell University, Ithaca, NY Here, we deeply profiled the transcriptome of thou- 14853, USA Full list of author information is available at the end of the article sands of individual MuSCs and muscle-resident cells © The Author(s). 2020 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Micheli et al. Skeletal Muscle (2020) 10:19 Page 2 of 13 from diverse adult human muscle samples using single- expression analysis between the 16 distinct subpopulations cell RNA-sequencing (scRNA-seq). After integrating identified an extensive set of unique markers that we these donor datasets to conserve biological information grouped into 4 categories (Fig. 1d). and overcome technical variation, we resolved two subpopulations of MuSCs with distinct gene expression scRNA-seq resolves the cellular diversity of human muscle signatures. Using differential gene expression analysis and novel markers and ligand-receptor interaction modeling, we extend the We annotated and interpreted the consensus cell atlas known repertoire of human MuSC gene expression (Fig. 1b, d) into cell type subpopulations as follows. We programs, suggesting new regulatory programs that may identify four types of stromal cells starting with adipocytes be associated with human MuSC activation, as well as found to be expressing apolipoprotein D (APOD)[15]), the features of human muscle aging and disease. brown fat tissue adipokine CXCL14 [16], GPX3,and GLUL. Among the 3 other subpopulations of fibroblast-like cells, Results Fibroblast 1 expresses high levels of collagen 1 (COL1A1), Collection and integration of a diverse human scRNA-seq SFRP4, SERPINE1,and CCL2; Fibroblast 2 expresses fibro- dataset nectin (FBN1), the microfibril-associated glycoprotein We used scRNA-seq to collect and annotate a single-cell MFAP5,and CD55 known to be expressed by synoviocytes transcriptomic dataset of diverse adult human muscle [17]; and Fibroblast 3 is mainly characterized by SMOC2 samples under homeostatic conditions. The muscle identified in tendon fibroblasts [18]. The Fibroblast 3 clus- samples were from surgically discarded tissue from n =10 ter is similar to the adipocytes cluster though exhibits lower donors (range 41 to 81 years old) undergoing reconstruct- expression levels and frequencies of the marker genes ive procedures and originating from a wide variety of APOD, CXCL12,and GLUL, and contain pre-adipocytes. anatomical sites in otherwise healthy patients (Fig. 1a). We also identify 5 types of vascular cells, including 3 Each sample was ~ 50 mg after removal of extraneous fat endothelial subpopulations, and a subpopulation of peri- and connective tissue. Muscle samples were enzymatically cytes and smooth muscle cells (SMCs). Pericytes and digested into single-cell suspensions and independently SMCs express the canonical markers RGS5 and MYH11. loaded into the 10X Chromium system. All together, we Endothelial 1 express E-selectin (SELE), IL6, ICAM1,and collected over 22,000 human muscle single-cell transcrip- VCAM1. These genes are upregulated at sites of inflamma- tomes (2206 ± 1961 cells per dataset) into a single data tion to facilitate immune cell recruitment, suggesting this compendium. Using unsupervised clustering, we resolved Endothelial 1 cell population may be involved in homeo- 16 types of cells of immune, vascular, and stromal origin, static muscle tissue remodeling [19, 20]. Endothelial 2 cells as well as two distinct subpopulations of MuSCs and some are distinguished by expressing high levels of claudin-5 myofiber myonuclei (Fig. 1b). (CLDN5), ICAM2, and the chemokine CXCL2.Endothelial Given important differences in anatomical site, donor 3 expresses high levels of the platelet-recruiting Von Will- health history, age, sex, and surgical procedures, the ebrand Factor (VWF) and caveolin-1 (CAV1), a protein muscle samples were highly heterogeneous in terms of known to regulate cholesterol metabolism, atherosclerosis cell-type diversity and underlying gene expression progression, and MuSC activation [21, 22]. Endothelial 3 profiles. Comparing the resulting scRNA-seq datasets is cells are enriched for expression of BTLN9, suggesting they therefore a challenge that we addressed using recently might represent a lymphatic endothelial phenotype [23]. developed bioinformatic integration methods [12–14]. We also noted two types of myeloid immune cells: first, Our goal was to assemble a unified dataset of human tissue-resident and anti-inflammatory macrophages which muscle tissue that faithfully conserved sources of bio- express CD74 and histocompatibility complex HLA proteins; logical variability such as donor, anatomical location, second, activated macrophages and monocytes that express and cell composition heterogeneity, while accounting for inflammatory markers such as S100A9 (calgranulin) and technical biases. We tested four different scRNA-seq LYZ (lysozyme). Moreover, S100A9 transcript abundance data integration methods (Fig. S1 and S3) and found that levels have been shown to be a feature in aging and chronic Scanorama [13] followed by scaling the output by regres- inflammation [24]. We also identified a pool of T/B lympho- sing against the library chemistry technical variable (“10X cytes and NK cells characterized by IL7R and NKG7,re- chemistry”) and the number of genes detected per single- spectively, as well as a small subset of HBA1 erythroblasts. cell best satisfied this goal. Detailed information on our Finally, we identified two subpopulations of MuSCs methodology is provided in Fig. S1.After integratingthe (henceforth called “MuSC1” and “MuSC2”). MuSC1 highly 10 datasets, we noted remarkable consistency amid cell expressed the canonical myogenic transcription factor PAX7 types across donors (Fig. 1c, e), owing to the robustness of [25], as well as chordin-like protein 2 (CHRDL2)and Delta- scRNA-seq technology, the bioinformatic method chosen, like non-canonical Notch ligand 1 (DLK1). CHRDL2 has and our sample preparation protocol. Differential gene been shown previously to be expressed in freshly isolated Micheli et al. Skeletal Muscle (2020) 10:19 Page 3 of 13 Fig. 1 (See legend on next page.) Micheli et al. Skeletal Muscle (2020) 10:19 Page 4 of 13 (See figure on previous page.) Fig. 1 Single-cell transcriptomic map of human muscle tissue biopsies. a Metadata (sex, age, anatomical site, and the number of single-cell transcriptomes after quality control (QC) filtering) from n = 10 donors. Colors indicate sample anatomical sites. b Scanorama-integrated and batch-effect corrected transcriptomic atlas revealing a consensus description of 16 distinct muscle-tissue cell populations. c Transcriptomic atlas colored by donor and anatomical location. d Dot-plot showing differentially expressed genes that distinguish the cell populations. Grouped in four compartments: muscle, endothelial/vascular, stromal, and immune. e Cell type proportions as annotated in (b) across the 10 donors and grouped by body sections. L, leg (donors 02, 07, 08); T, trunk (donors 01, 05, 06, 09, 10); F, face (donors 03, 04) quiescent human MuSCs [7], though its function is still to In addition, IL32 has been shown to have inflammatory be understood. DLK1 is an inhibitor of adipogenesis whose properties in human obesity [37] and have a negative impact role in muscle has mainly been recognized in the on insulin sensitivityand myogenesis [38], while embryo but remains controversial in adult muscle re- TNFRSF12/FN14 has been implicated in various muscle generation [26–28]. In contrast to MuSC1, MuSC2 wasting diseases [39, 40] and metabolic dysfunction [41]. expressed lower levels of PAX7 but maintain expression Furthermore, the MuSC2 population is enriched for riboso- of MYF5 (a marker of activated MuSCs) and APOC1 mal gene expression (e.g. RPLP1 and RPS6; data not shown), (Fig. 2b). Interestingly, the MuSC2 population also had indicating that these cells may have elevated translational elevated expression of two long non-coding RNAs mechanisms. Lastly, the MuSC1 population has enriched ex- (lncRNAs), LINC00152,and MIR4435-2HG.LncRNAs pression of the myogenic gene PAX7 and, to a lesser extent, are involved in regulating myogenesis [29]. Surprisingly, MYF5, compared the MuSC2 population. These observa- we detected low expression of the myogenic commitment tions suggest that MuSC1 is comprised of quiescent MuSCs, factors MYOD1 and MYOG (Fig. 2b), in contrast to scRNA- and MuSC2 is comprised of an early-activated MuSCs. seq analyses of adult mouse muscle [30, 31]. These observa- We performed Ingenuity Pathway Analysis (IPA) to com- tions suggest that the MuSC1 and MuSC2 populations are pare biological processes differentially activated between both comprised largely of muscle stem cells, not commit- the MuSC1 and MuSC2 populations. The IPA gene group ted myogenic progenitors. In addition, we noted that “Oxidative Phosphorylation” is enriched in MuSC1 [42], “Myonuclei” population (Fig. 1b) was enriched for myosin while “EIF2 Signaling,” associated with protein translation light chain (MYLFP), skeletal alpha-actin (ATCA1), and processes, is enriched in MuSC2 (Fig. 2c). Furthermore, troponin C (TNNC2), proteins involved in muscle con- Gene Set Enrichment Analysis (GSEA) also found MuSC1 traction. This multiple-donor scRNA-seq atlas highlights to be enriched for “myogenesis,”“muscle cell differenti- the cellular diversity of human muscle tissue and revealed ation,”“hypoxia,” and “response to mechanical stimulus” two distinct MuSC subpopulations along with specific gene sets, supporting the observation that these cells are myogenic expression programs. both less differentiated and may exhibit enhanced tran- scriptional responses to mechanical disruption due to tissue Homeostatic human muscle contains two distinct MuSC dissociation [32–34](Fig. 2d). MuSC2 cells are enriched for subpopulations “ribosome and translational initiation,”“MYC targets,” and We examined genes that were differentially expressed be- “E2F (cell proliferation),”“G2M checkpoint (cell division),” tween the MuSC1 and MuSC2 subpopulations and the bio- and “inflammation” gene sets, further supporting the inter- logical processes that characterize them (Fig. 2a, b). The pretation that these cells may be in an early activated or MuSC1 subpopulation was enriched for PAX7, DLK1,and partially differentiated state within an inflammatory envir- CHRDL2, as well as for the cyclin-dependent kinase inhibi- onment (Fig. 2d). Taken together, these observations KIP2 tor CKDN1C (encoding P57 ), suggesting that these cells suggest that the MuSC1 population is comprised of quies- are quiescent and not cycling. In addition, this subpopula- cent MuSCs, while the MuSC2 population is comprised of tion expresses the transcription factor BTG2,which was active, proliferating, and/or dysregulated MuSCs, with identified in mouse to be enriched in quiescent MuSCs [30]. expression alterations associated with inflammation, aging, We also note that the MuSC1 subpopulation expressed and muscle wasting. Differentially expressed genes such as elevated levels of mitochondrial genes as well as FOS, JUN, IL32, CXCL1, CCL2,and TNFRSF12/FN14 may constitute and ERG1. Upregulation of these genes has been shown to a marker set for MuSC variation in chronic muscle inflam- be potential artefacts of the enzymatic digestion during the mation in various pathologies. sample preparation [32–34]. The MuSC2 subpopulation was enriched for multiple Ligand-receptor interaction model identifies potential markers of inflammation including CCL2, CXCL1, IL32,and surface markers and cell-communication channels in surface receptor TNFRSF12/FN14. In particular, CCL2 and human skeletal muscle homeostasis CXCL1 are inflammatory cytokines known to be upregu- We used a ligand-receptor (LR) interaction model and a lated in muscle repair, exercise, and fat metabolism [35, 36]. database of LR pairs [43] to map cell signaling Micheli et al. Skeletal Muscle (2020) 10:19 Page 5 of 13 Fig. 2 Gene expression and pathway analysis comparison between two MuSC subpopulations. a Volcano plot from comparing transcript levels between all cells within the “MuSC1” and “MuSC2” subpopulations. Log fold-change in normalized gene expression versus −log adjusted p value 2 10 plotted. Differentially expressed genes (adjusted p value <0.05) are colored dark or light blue (based on their enrichment in MuSC1 or MuSC2, respectively). Genes with log fold-change > 0.75 are labeled. b Normalized expression values of select differentially expressed genes. q values reported in inset. c Top activated canonical pathways by Ingenuity Pathway Analysis (IPA) based on differentially expressed genes and ranked by p value. Pathways significantly enriched in either population with |z score| > 1 are indicated in blue. d Select gene ontology (GO) terms and hallmark pathways enriched between the MuSC subpopulations as identified by gene set enrichment analysis (GSEA) and ranked by enrichment score (ES) communication channels in human muscle and uncover expressed receptors on a given cell population (e.g., differences between MuSC1 and MuSC2 subpopula- “MuSC1”) relative to all other population and ligands tions (Fig. 3). The model also identifies interacting li- expressed by other cell types. The MuSC1 and MuSC2 gand(s) and is restricted to receptor genes differentially subpopulations are involved in numerous LR interac- expressed by a specific cell type within the consensus tions, as both ligand- and receptor-expressing cells (Fig. human muscle cell atlas (Fig. 1b). For each LR pair, the 3a), though a majority of all LR interaction pairs instead model calculates an interaction score from differentially involve other cell types. This suggests that only a small Micheli et al. Skeletal Muscle (2020) 10:19 Page 6 of 13 Fig. 3 (See legend on next page.) Micheli et al. Skeletal Muscle (2020) 10:19 Page 7 of 13 (See figure on previous page.) Fig. 3 Differentially expressed receptors and ligand-receptor interaction between cell populations. a Chord plot of all ligand-receptor (LR) interactions across cell populations/types within the consensus atlas based on co-expression. Each cell type is color-coded with its receptor genes more displaced from the perimeter than its ligand genes. All interactions not involving MuSC1 or MuSC2 are presented in gray. b List of differentially expressed genes between the MuSC1 and MuSC2 subpopulation ranked by log fold-change in expression. Positive average values correspond to genes that are upregulated in MuSC1, whereas negative values are upregulated in MuSC2. Receptors that are statistically significant (FDR-corrected q value < 0.05) are colored in blue. Receptors that are not statistically significant are in gray. c Heatmap representing row- normalized (Z-score) LR interaction scores. Rows represent ligand-receptor interaction pairs in the format LIGAND_RECEPTOR, where the receptor is either differentially expressed in the MuSC1 or MuSC2 populations compared to all the other cell types. Columns identify cell types expressing the ligand. Asterisks after the pair name also indicate that the ligand is differentially expressed by the other cell type and that interaction is likely cell-type specific. Red pairs involve the EGFR receptor, purple pairs the NOTCH3 receptor. A positive value indicates that the interaction has a high score for a particular ligand and cell type compared to other cell types subset of potential paracrine interactions in human (EREG). Other EGFR ligands include brevican (BCAN), and muscle may include MuSCs. betacellulin (BTC) produced by endothelial cells; ECM pro- Given the distinct expression profiles between the teins fibulin 3 (EFEMP1), decorin (DCN), and tenascin C MuSC1 and MuSC2 populations, we sought to identify (TNC) expressed by fibroblasts; and FGF13, AHM, NRG4, genes that could facilitate surface antigen-based separation and EGF, expressed by mature skeletal myofibers. We also of these two human MuSC populations for prospective iso- detect seven interactions involving NOTCH3 with a variety lation strategies. We identified surface receptor genes that of ligands. Notch3 signaling is involved in maintaining were differentially expressed between the MuSC1 and MuSC quiescence, in particular through interaction with MuSC2 populations, using a database of 542 human sur- DLL4 [50], whichwefound differentlyexpressedby endo- face “receptor” genes [43](Fig. 3). MuSC1 exhibit elevated thelial cells along with JAG2. In addition, NOTCH3 also in- expression of EGFR, ITGB1, FGFR4, SDC2,aswellasthe teracts with the ECM protein thrombospondin-2 (THBS2). three tetraspanins CD81, CD82,and CD151(Fig. 3b). EGFR Only two receptors, TNFRSF12/FN14 and RPSA,were is a recently established human MuSC marker and is re- found differentially expressed in MuSC2 compared to quired for basal-apical asymmetric cell division [7, 10]. The other cell types. The first, TNFRSF12/FN14, interacts with tetraspanin CD82 is also a recently recognized human the TWEAK cytokine ligand. While typically recognized to MuSC maker [6], while CD9 and CD81 have been identi- be expressed by macrophages and other immune cells fied to control muscle myoblast fusion [44]. Furthermore, [51], our model suggests that TWEAK is also expressed by Syndecans (SDCs) have been identified in mouse to be het- the Fibroblast 2 and pericyte cell populations, though not erogeneously expressed on MuSCs and myoblasts during in a statistically significant manner. The second, RPSA,is muscle repair [30] and have been shown to form co- surface ribosomal protein that interacts with laminins receptor complexes with integrin β1 (ITGB1) and FGFR4 (LAM), a dual-specificity phosphatase 18 (DUSP18), and upstream of signaling pathways regulating myogenesis prion protein 2 (PRND), which taken together may sug- [45]. Only SDC4 and SDC3 have yet been identified to gest various pathological processes such as prion diseases mark adult mouse MuSCs [46]. In comparison, the MuSC2 and cancer [52, 53]. Together, this ligand-receptor analysis subpopulation has elevated expression of CD44 and identifies a broad set of surface markers that could refine TNFRSF12/FN14 as previously noted (Fig. 3b). The CD44 the molecular definition of human MuSCs and their receptor has been shown to regulate myoblast migration subpopulations, as well as candidate cell-communication andfusioninmouse,but also mark MuSCsinosteoarthritis channels differentially involved in healthy and diseased patients [47, 48]. muscle tissues. Next, we focused the LR analysis on the MuSC1 and Lastly, we performed a comparative analysis of receptor MuSC2 populations. We identified 73 and 6 significant LR gene expression between mouse and human MuSCs. We interactions for the MuSC1 and MuSC2 populations, integrated the human scRNA-seq datasets described in respectively (Fig. 3c). Over one third of all interactions in Fig. 1 and an adult mouse muscle injury-response scRNA- the MuSC1 subpopulation involve the EGFR receptor, seq time-course previously reported [30]byconverting which has recently been shown to play a critical role in mouse genes to their corresponding human ortholog. The directing MuSC asymmetric division in regenerating multi-species scRNA-seq atlas was integrated with Sca- muscle [10]. A limited number of EGFR ligands have been norama and corrected with Harmony (Fig. S2A-B) [54]. identified in muscle repair, for example, amphiregulin. From this integrated atlas, we annotated all clusters as in (AREG) secreted by T cells [49]. According to our model Fig. 1. We identified two MuSC clusters which both con- reg findings, EGFR may also interact with ligands expressed by tained cells from both mouse and human samples. We immune cells, such as with TGF-α (TGFA), heparin-biding then performed differential expression analysis between EGF (HBEGF), amphiregulin (AREG), and epiregulin species comparing aggregated human MuSC1 and MuSC2 Micheli et al. Skeletal Muscle (2020) 10:19 Page 8 of 13 cells to mouse MuSCs from the uninjured timepoint (Fig. due to its differential expression of BTLN9, a marker of S2C). We found that EGFR and CD99 were most differen- lymphatic endothelial cells [23]. tially expressed by human MuSCs and, conversely, CRLF1 Most notably, this analysis suggests that human and SDC4 were most enriched in mouse MuSCs. This muscle may contain two distinct MuSC subpopulations findings suggest that mouse and human MuSC exhibit (Fig. 2). This finding contrasts with Rubenstein et al. species-specific receptor expression signatures. which observed a single MuSC (“satellite cell”) popula- tion from dissociated whole muscle samples and Barruet et al. which observed ~ 12 clusters from human MuSCs Discussion prospectively enriched by CXCR4+/CD29+/CD56+ Here we present an annotated multi-donor single-cell FACS. Since cluster distinction depends on both the cel- RNA sequencing dataset consisting of 22,000 single-cell lular diversity and sample complexity, it is expected that transcriptomes from 10 different donors and unique variation in study design and methods will yield differing anatomical sites, some of which difficult to access out- conclusions regarding sub-population resolution. In this side of reconstructive surgeries. Our study complements work, we found a “MuSC1” subpopulation to be largely other recent reports by Rubenstein et al. and Barruet comprised of “quiescent” MuSCs, owning to high levels et al., which collected dissociated whole vastus lateralis of PAX7, the mitotic inhibitor CDKN1C, and DLK1. muscles and FACS-sorted MuSC samples mostly from Interestingly, DLK1 may be an important regulator for vastus lateralis muscles, respectively, by providing more human MuSC maintenance and a marker of healthy tis- diversity in anatomical sites and donor demographics sue given its role in inhibiting adipogenesis [26]. Con- [55, 56]. As such, these scRNA-seq data exhibited versely, we identified in the “MuSC2” population notable biological and technical variation, and therefore, signatures of inflammation and increased fat metabolism we applied the bioinformatic method Scanorama to assem- (CCL2 and CXCL1), reduced insulin sensitivity (IL32), ble an integrated cellular atlas with minimal technical cell cycle (EIF2 Signaling terms), and muscle wasting biases so that we could examine the cellular heterogeneity (TNFRSF12/FN14), thereby suggesting that these cells across diverse adult human muscle tissue samples. We ob- may constitute an “early-activated” and possibly dysfunc- served that Scanorama performed more successfully than tional MuSC pool. These markers are consistent with other data integration approaches, especially when includ- prior observations that excessive fat accumulation in ing a scaling regression for sequencing chemistry (Fig. S1 muscle can be attributed to obesity, diabetes, and aging and S3). Notably, even after performing Scanorama with [4]. In addition, we identify two upregulated lncRNAs scaling, we still observed that integrated atlas exhibited bio- that warrant further investigation as candidate non- logical (donor) and technical (sequencing chemistry) coding regulators of myogenesis [29]. Moreover, the biases, but retained some degree of donor-specific cell-type finding of two human MuSC subpopulations mirrors subpopulations. similar observations made from mouse muscle scRNA-seq We describe the muscle tissue cellular heterogeneity analyses [30, 31] and agrees with the general conceptual and provide a comprehensive analysis of differentially framework that MuSCs transition between quiescent, acti- expressed genes for 16 resolved cell subpopulations (Fig. vated, and cycling states [1]. Future studies comparative 1), adding to a growing literature documenting human analysis of these MuSC subpopulations across species may muscle cell transcriptional diversity [55–57]. Compared reveal human-specific aspects of myogenesis. to other studies, the broader variety of muscle tissue Ligand-receptor interaction models from scRNA-seq samples combined with the lack of FACS selection data can help formulate new hypotheses about cell- allowed us to identify candidate subpopulations of communication channels that regulate muscle function muscle fibroblasts and vascular endothelial cells that [30]. Identifying new MuSCs surface receptors will also may provide unique perspective to human muscle physi- help us refine MuSC purification protocols for prospect- ology. In particular, we remark that Endothelial 1 ive isolation studies used for in vitro and transplantation expressed DARC/ACKR1, a gene identified in mouse and models. Our LR model revealed a set of 40 surface human [56, 58] to mark cells of post-capillary venular receptor genes that are distinctly expressed between origin (Fig. 1d). Rubenstein et al. also found a DARC/ MuSC1 and MuSC2, confirming some prior reports but ACKR1+ post-capillary venular endothelial cluster and a also providing new candidate surface antigens for human second VWF+ FABP+ cluster, which overlaps with the MuSC subpopulation fractionation (Fig. 3). For example, Endothelial 2 and 3 clusters reported here. We suggest we identify that SDC2 may mark “quiescent” MuSCs that the Endothelial 2 cluster may contain both arterial while CD44, TNFRSF12, and RPSA “early-activated” and capillary endothelial cells, but could not further par- MuSCs in aging and disease contexts. In addition, our tition and classify this cluster. We suggest that the Endo- model proposed 79 cell-communication signals that may thelial 3 cluster may represent lymphatic endothelium act between MuSCs and other cell types, in particular Micheli et al. Skeletal Muscle (2020) 10:19 Page 9 of 13 with fibroblasts, myofibers and immune cells through would not be possible for some muscles presented in this the EGFR receptor, and with vascular cells through the study. These biopsies would allow for aging and disease NOTCH3 receptor. These interactions may be critical comparative analyses. Indeed, a recent report by Ruben- regulators of muscle homeostasis and should be further stein et al. [56] performed scRNA-seq on four human vas- investigated. tus lateralis muscle biopsies found that myofiber type This study presents a new set of candidate receptor composition and gene expression alterations based on expression signatures that may define human MuSC sub- donor age. populations (Fig. 3b) and provide human-specific receptor Nevertheless, our dataset offers a new transcriptomic patterns (Fig. S2C). This approach is complimentary to re- cell reference atlas and computational data integration ap- ceptor screening approaches, which have previously been proaches as a benchmark resource to examine human useful to identify EGFR and CD82 as human MuSC muscle cell diversity in health, aging, and disease. receptor markers for flow cytometry [6, 7, 9]. The subpopulation-specific receptor genes identified here may Methods allow for further comparison of molecular and functional Human participation for muscle sample collection human MuSC diversity across muscle groups [59, 60]. All procedures were approved by the Institutional Review Our study has some limitations. First, the sample size is Board at Weill Cornell Medical College (WCMC IRB small, and donors are very diverse, thus limiting our ability Protocol # 1510016712) and were performed in accordance to control for age and sex. Therefore, we could not exam- with relevant guidelines and regulations. All specimens ine cell composition or gene expression trends based on were obtained at the New York-Presbyterian/Weill Cornell muscle group, donor sex, or donor age. Even for samples campus. All subjects provided written informed consent from the same muscle (e.g., flexor hallucis longus [donors prior to participation. Samples were de-identified in 2 and 7] or external oblique [donors 6 and 9]), we were un- accordance to IRB guidelines, and only details concerning able to perform these comparions with statistical power. age, sex, and anatomic origin were included. Sample ana- Further, we performed differential expression and gene set tomic locations and donor details are provided in Fig. 1a. enrichment analyses within the MuSC1 and MuSC2 popu- lations between the four middle-age (43–69 years old) and Muscle digestion and single-cell sequencing library six aged (70–81 years old) donors, but found few age- preparation cohort specific differences (data not shown). Second, future After collection from donors during surgery, the muscle studies should aim at collecting muscle specimens in a samples were cleared from excessive fat and connective more controlled manner, for example using a Bergström tissue and weighted. About 50–65 mg of tissue was then needle [61, 62] from a unique anatomical site; though this digested into a single-cell suspension following a previously Table 1 List of reagents and other resources used in this study Reagents Dispase II (neutral protease, grade II) Sigma-Aldrich 04942078001 Collagenase D, from Clostridium histolyticum Sigma-Aldrich 11088866001 Commercial kits Chromium Single Cell 3' Library & Gel Bead Kit v2 10X Genomics CG00052 (protocol) Chromium Single Cell 3' Library & Gel Bead Kit v3 10X Genomics CG000183 (protocol) Deposited data Human ligand-receptor database [43] https://www.ncbi.nlm.nih.gov/pubmed/26198319 Human scRNA-seq dataset This paper GSE143704 Mouse scRNA-seq dataset [29] GSE143437 Software packages and algorithms Cell Ranger 3.1.0 (July 24, 2019) 10X Genomics https://support.10xgenomics.com/single-cell-gene-expression/software/ downloads/latest Seurat 3.1.0 [14] https://github.com/satijalab/seurat Scanorama (online version as of 2019-11-19) [13] https://github.com/brianhie/scanorama Harmony (online version as of 2019-11-19) [54] https://github.com/immunogenomics/harmony biomaRt 2.43.1 (online version as of 2020-01-08) [64] https://bioconductor.org/packages/release/bioc/html/biomaRt.html Gene Set Enrichment Analysis (4.0.3) [65] http://software.broadinstitute.org/gsea/index.jsp Ingenuity Pathway Analysis (IPA, 2019-08-30) QIAGEN https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis/ Micheli et al. Skeletal Muscle (2020) 10:19 Page 10 of 13 reported protocol [63]. Briefly, the specimen was digested andthe uninjuredmouse MuSCsto focus on cellsfrom in 8 mg/mL Collagenase D (Roche) and 4.8 U/mL Dispase the homeostatic conditions. II (Roche) for 1 h followed by manual dissociation, filtra- tion, and red blood cell lysis (Table 1). All single-cell sus- Pathway and gene set enrichment analysis pensions were then frozen at -80 °C in 90% FBS, 10% The list of differentially expressed genes between DMSO and were re-filtered after thawing and prior to gen- MuSC1 and MuSC2 (Fig. 2a) was used in Ingenuity erating scRNA-seq libraries. The sequencing libraries were Pathway Analysis (IPA) (QIAGEN, 2019-08-30). Acti- prepared using the Chromium Single Cell 3' reagent V2 or vated (canonical) pathways were calculated by “Core V3 kit (10X Genomics) in accordance with the manufac- Analysis” setting a q value cutoff of 0.05, which yielded turer’s protocol and diluted as to yield a recovery of ~ 6000 964 genes (366 down, 598 up). Top canonical pathways single-cell transcriptomes with < 5% doublet rate (Table 1). were chosen based of − log(p value) and z score values. The libraries were sequenced in multiplex (n =2perse- Gene setenrichmentanalysis (GSEA, v.4.0.3)[65]was quencing run) on the NextSeq 500 sequencer (Illumina) to ran on the same gene list as IPA ranked by log fold- produce between 200 and 250 million reads per library. change and with default program settings (Table 1). Gene sets database used the following: h.all.v7.0.sym- bols.gmt, c2.all.v7.0.symbols.gmt, c5.all.v7.0.sym- Single-cell data analysis bols.gmt (Broad Institute). Gene sets enriched in Sequencing reads were processed with the Cell Ranger phenotype were selected based on q value and enrich- version 3.1 (10X Genomics) using the human reference ment score (ES). transcriptome GRCh38. The downstream analysis was carried out with R 3.6.1 (2019-07-05). Quality control Ligand-receptor cell communication model filtering, data clustering, visualization, and differential The model aims at scoring potential ligand-receptor gene expression analysis was carried out using Seurat interactions between MuSCs (receptor) and other cell 3.1.0 R package [14]. Each of the 10 datasets was first types (ligand). We used the ligand-receptor interaction analyzed and annotated independently before integration database from Ramilowski et al. [43] (Table 1). From the with Scanorama [13] (Table 1). Filtering retained cells database, we considered 1915 ligand-receptor pairs with > 1000 unique molecular identifiers (UMIs), < 20% (from 542 receptors and 518 ligands) to test for differen- UMIs mapped to mitochondrial genes, and genes tial expression in our scRNA-seq dataset. To calculate expressed in at least 3 cells (Fig. S4). Unsupervised the score for a given ligand-receptor pair, we multiply shared nearest neighbor (SNN) clustering was performed the average receptor expression in MuSCs by the aver- with a resolution of 0.4 following which clusters were age ligand expression per other cell type. We only con- annotated with a common nomenclature of 12 cell type sidered receptors that are differentially expressed in terms (Fig. S1). Differential expression analysis was either the MuSC1 or MuSC2 subpopulation when com- achieved using either Seurat’s “FindAllMarkers” (Fig. 1d) pared individually to all other cell types. or “FindMarkers” (Fig. 2a) function using a Wilcoxon Rank Sum test and only considering genes with > Supplementary information log (0.25) fold-change and expressed in at least 25% of Supplementary information accompanies this paper at https://doi.org/10. cells in the cluster. P values were corrected for false- 1186/s13395-020-00236-3. discovery (FDR) and then reported as q values. Integration of raw counts was achieved using the “scanorama.correct” Additional file 1: Figure S1. Comparison of scRNA-seq integration and batch correction methods. We compared four scRNA-seq data in- function from Scanorama. The integrated values were fi- tegration methods to evaluate which most faithfully conserves donor, nally scaled in Seurat regressing out the 10X chemistry anatomical, and biological information while minimizes technical biases. type and the number of genes per cell. Visualization was (A) The n = 10 donor datasets were first annotated independently using a nomenclature of 12 common cell type terms following unsupervised done using uniform manifold approximation and projec- SNN clustering. Then we evaluated the integration method by UMAP and tion (UMAP) [66]. In Fig. S2, we integrated these human by coloring the data either by cell type, donor ID, or 10X library chemistry scRNA-seq datasets with a cohort of adult mouse muscle used. First, we integrated the data by merging the individually normalized gene expression matrices without any further correction. We saw strong scRNA-seq datasets collected 0–7 days post-notexin injury technical biases that overwhelmed biological information as the different [30]. For multi-species integration, scRNA-seq datasets cell populations segregate by sample/donor and chemistry type. For in- were integrated using first Scanorama and then Harmony stance, the two MuSC and progenitor subpopulations are grouped with fibroblasts and endothelial cells. Second, we tested the Seurat SCT inte- [54] to align related cell populations across species. Mouse gration method [14] . This method first calculates a cross-correlation sub- genes were converted to human orthologs using biomaRt space from genes that are shared between datasets. We noticed that this Bioconductor R package [64](Table 1). For differential ex- method better “aligns” donor and chemistry type but at the expense of masking biological variability. For instance, we observed that the two pression analysis between human and mouse samples, we MuSC and four stromal subpopulations (Fibroblast 1,2,3 and Adipocytes) compared human MuSCs (combining MuSC1 + 2 clusters) Micheli et al. Skeletal Muscle (2020) 10:19 Page 11 of 13 Availability of data and materials were grouped together, hiding important biological heterogeneity. Al- The human muscle scRNA-seq datasets supporting the conclusions of this though certainly useful to validate reproducibility in scRNA-seq experi- article are archived at the NIH GEO repository under accession number ments, the Seurat SCT integration approach overcorrected biological GSE143704. heterogeneity for heterogeneous samples. Third, we tested the Scanor- ama method [13], which relies on a computer vision algorithm that Ethics approval and consent to participate “stitches” datasets together even when the cell type composition be- All procedures were approved by the Institutional Review Board at Weill tween dataset is considerably different. We see that this method groups Cornell Medical College (WCMC IRB Protocol # 1510016712) and were similar cell populations together while acknowledging donor differences. performed in accordance with relevant guidelines and regulations. All Yet, surprisingly, this method is also very sensitive at picking up differ- specimens were obtained at the New York-Presbyterian/Weill Cornell ences in chemistry. To correct this chemistry effect, we scaled the Scanor- campus. All subjects provided written informed consent prior to ama output by regressing out the chemistry and the number of genes participation. Samples were de-identified in accordance to IRB guidelines detected per cell (significantly different between chemistry type) (B). and only details concerning age, sex, and anatomic origin were included. Using this integration method, we observed clear separation of the inde- pendently annotated cell populations. We present the resulting Consent for publication Scanorama-integrated dataset as a “consensus atlas” (see Fig. 1b-c) of hu- man muscle that describes donor-to-donor differences while grouping Not applicable. cells that are similar together and removing technical biases. Figure S2. Integration of human and mouse scRNA-seq data sets allows com- Competing interests parison of MuSC receptor gene expression across species. We gen- The authors declare no conflicts of interest. erated an integrated scRNA-seq atlas including human sample datasets from Fig. 1 and an adult mouse muscle regeneration time-course from Author details De Micheli et al. [29]. These datasets were integrated using first Scanor- Meinig School of Biomedical Engineering, Cornell University, Ithaca, NY ama and then Harmony for alignment across species. (A) Multi-species in- 14853, USA. Englander Institute for Precision Medicine, Weill Cornell tegrated atlas presented by UMAP plot a colored by sample type. (B) Medicine, New York, NY 10021, USA. Division of Plastic Surgery, Weill Cornell Multi-species integrated atlas presented by UMAP plot and annotated by Medical College, New York, NY 10021, USA. cell-type clusters. (C) The human MuSC1 and MuSC2 clusters were grouped into a cumulative human MuSC cell population, which was Received: 30 January 2020 Accepted: 10 June 2020 compared to mouse MuSCs from the uninjured samples only. Receptor genes were analyzed between the mouse and human MuSC cells for dif- ferential expression. Differentially expressed genes with an FDR-corrected References q-value < 0.05 are shown in (C). Figure S3. Composition of single-cell 1. Bentzinger CF, Wang YX, Dumont NA, Rudnicki MA. Cellular dynamics in the reference atlas as a whole and in cell-type clusters by donor. (A) muscle satellite cell niche. EMBO Rep. 2013;14:1062–72. Visualization of donor (n = 10) contributions to the whole single-cell ref- 2. Blau HM, Cosgrove BD, Ho ATV. The central role of muscle stem cells in erence atlas. In each panel, the full atlas is presented as a UMAP plot, regenerative failure with aging. Nat Med. 2015;21:854. with the cells for an individual donor are colored and overlaid on cells 3. Järvinen TA, Järvinen M, Kalimo H. Regeneration of injured skeletal muscle from all other donors (in gray). Note the total number of cells assayed dif- after the injury. Muscles Ligaments Tendons J. 2014;3:337–45. fers for each donor (see Fig. 1a). (B) Bar plot representing the relative 4. Addison O, Marcus RL, LaStayo PC, Ryan AS. Intermuscular fat: a review of contribution of cells with each cell type cluster from each donor. Note the consequences and causes. Int J Endocrinol. 2014;2014:1–11. that the MuSC1 and MuSC2 clusters are also plotted as a combined clus- 5. Larsson L, Degens H, Li M, Salviati L, Lee YI, Thompson W, Kirkland JL, Sandri ter on the left side of the bar plot for reference. Figure S4. Transcrip- M. Sarcopenia: aging-related loss of muscle mass and function. Physiol Rev. tomic detection variation within human muscle reference atlas. 2018;99:427–511. UMAP plots featuring (left) the number of unique molecular identifiers 6. Alexander MS, Rozkalne A, Colletta A, Spinazzola JM, Johnson S, Rahimov F, (UMIs) and (right) number of genes detected per single cell. Note that QC Meng H, Lawlor MW, Estrella E, Kunkel LM, et al. CD82 is a marker for filtering removed all cells with less than 1000 UMIs (see Methods). prospective isolation of human muscle satellite cells and is linked to muscular dystrophies. Cell Stem Cell. 2016;19:800–7. 7. Charville GW, Cheung TH, Yoo B, Santos PJ, Lee GK, Shrager JB, Rando TA. Ex vivo expansion and in vivo self-renewal of human muscle stem cells. Stem Cell Reports. 2015;5:621–32. Acknowledgements 8. Pisani DF, Clement N, Loubat A, Plaisant M, Sacconi S, Kurzenne J-Y, The authors acknowledge helpful advice from colleagues in the Cosgrove Desnuelle C, Dani C, Dechesne CA. Hierarchization of myogenic and and Elemento groups, as well as Christopher Mendias at the Hospital for adipogenic progenitors within human skeletal muscle. Stem Cells. 2010;28: Special Surgery and Peter Schweitzer of Genomics Facility at the Cornell 2182–94. University Biotechnology Resource Center. Lastly, the authors are grateful for 9. Uezumi A, Nakatani M, Ikemoto-Uezumi M, Yamamoto N, Morita M, the human tissue donors. Yamaguchi A, Yamada H, Kasai T, Masuda S, Narita A, et al. Cell-surface protein profiling identifies distinctive markers of progenitor cells in human skeletal muscle. Stem Cell Reports. 2016;7:263–78. Authors’ contributions 10. Wang YX, Feige P, Brun CE, Hekmatnejad B, Dumont NA, Renaud J-M, A.J.D. and B.D.C. designed the study and wrote the manuscript. J.A.S. Faulkes S, Guindon DE, Rudnicki MA. EGFR-Aurka signaling rescues polarity obtained the human tissue samples. A.J.D. performed the tissue dissociations, and regeneration defects in dystrophin-deficient muscle stem cells by scRNA-seq, and data analysis, with supervision and assistance from B.D.C. and increasing asymmetric divisions. Cell Stem Cell. 2019;24:419–432.e6. O.E. All authors reviewed, read, and approved the final manuscript. 11. Sousa-Victor P, Gutarra S, García-Prat L, Rodriguez-Ubreva J, Ortet L, Ruiz-Bonilla V, Jardí M, Ballestar E, González S, Serrano AL, et al. Geriatric muscle stem cells switch reversible quiescence into senescence. Nature. 2014;506:316–21. Funding 12. Stuart T, Satija R. Integrative single-cell analysis. Nat Rev Genet. 2019;20:257–72. This work was financially supported by the National Institutes of Health 13. Hie B, Bryson B, Berger B. Efficient integration of heterogeneous single-cell under award R01AG058630 (to B.D.C.), a Glenn Medical Research Foundation transcriptomes using Scanorama. Nat Biotechnol. 2019;37:685–91. and American Federation for Aging Research Grant for Junior Faculty (to 14. Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM, Hao Y, B.D.C.), and a US Department of Education Graduate Assistantship in Areas of Stoeckius M, Smibert P, Satija R. Comprehensive integration of single-cell National Need under Award P200A150273 (to A.J.D.). The content is solely data. Cell. 2019;177:1888–1902.e21. the responsibility of the authors and does not necessarily represent the 15. Muffat J, Walker DW. Apolipoprotein D: an overview of its role in aging and official views of any of these funding sources. age-related diseases. Cell Cycle. 2010;9:269–73. Micheli et al. Skeletal Muscle (2020) 10:19 Page 12 of 13 16. Cereijo R, Gavaldà-Navarro A, Cairó M, Quesada-López T, Villarroya J, Morón- 37. Catalán V, Gómez-Ambrosi J, Rodríguez A, Ramírez B, Ortega VA, Ros S, Sánchez-Infantes D, Peyrou M, Iglesias R, Mampel T, et al. CXCL14, a Hernández-Lizoain JL, Baixauli J, Becerril S, Rotellar F, Valentí V, et al. IL-32α- brown adipokine that mediates brown-fat-to-macrophage communication induced inflammation constitutes a link between obesity and colon cancer. in thermogenic adaptation. Cell Metab. 2018, 28:750–763.e6. Oncoimmunology. 2017;6:e1328338. 17. Karpus ON, Kiener HP, Niederreiter B, Yilmaz-Elis AS, van der Kaa J, Ramaglia 38. Davegårdh C, Broholm C, Perfilyev A, Henriksen T, García-Calzón S, Peijs L, V, Arens R, Smolen JS, Botto M, Tak PP, et al. CD55 deposited on synovial Hansen NS, Volkov P, Kjøbsted R, Wojtaszewski JFP, et al. Abnormal collagen fibers protects from immune complex-mediated arthritis. Arthritis epigenetic changes during differentiation of human skeletal muscle stem Research & Therapy. 2015;17:6. cells from obese subjects. BMC Med. 2017;15:39. 39. Enwere EK, Lacasse EC, Adam NJ, Korneluk RG. Role of the TWEAK-Fn14- 18. De Micheli AJ, Swanson JB, Disser NP, Martinez LM, Walker NR, Oliver DJ, Cosgrove BD, Mendias CL. Single-cell transcriptomics identify extensive cIAP1-NF-κB signaling axis in the regulation of myogenesis and muscle heterogeneity in the cellular composition of mouse Achilles tendons. homeostasis. Front Immunol. 2014;5:34. BioRxiv. 2020b;801266. 40. Mittal A, Kumar A, Lach-Trifilieff E, Wauters S, Li H, Makonchuk D, Glass D, 19. Goncharov NV, Nadeev AD, Jenkins RO, Avdonin PV. Markers and Kumar A. The TWEAK-Fn14 system is a critical regulator of denervation- biomarkers of endothelium: when something is rotten in the state. induced skeletal muscle atrophy in mice. J Cell Biol. 2010;188:833–49. Oxidative Med Cell Longev. 2017;2017:9759735. 41. Sato S, Ogura Y, Kumar A. TWEAK/Fn14 signaling axis mediates skeletal 20. Watson C, Whittaker S, Smith N, Vora AJ, Dumonde DC, Brown KA. IL-6 acts muscle atrophy and metabolic dysfunction. Front Immunol. 2014;5:18. on endothelial cells to preferentially increase their adherence for 42. Ryall JG, Dell’Orso S, Derfoul A, Juan A, Zare H, Feng X, Clermont D, Koulnis lymphocytes. Clin Exp Immunol. 1996;105(1):112–9. M, Gutierrez-Cruz G, Fulco M, et al. The NAD + -dependent SIRT1 21. Fernández-Hernando C, Yu J, Dávalos A, Prendergast J, Sessa WC. deacetylase translates a metabolic switch into regulatory epigenetics in Endothelial-specific overexpression of caveolin-1 accelerates atherosclerosis skeletal muscle stem cells. Cell Stem Cell. 2015;16:171–83. in apolipoprotein E-deficient mice. Am J Pathol. 2010;177:998–1003. 43. Ramilowski JA, Goldberg T, Harshbarger J, Kloppmann E, Lizio M, 22. Volonte D, Liu Y, Galbiati F. The modulation of caveolin-1 expression controls Satagopam VP, Itoh M, Kawaji H, Carninci P, Rost B, et al. A draft network of satellite cell activation during muscle repair. FASEB J. 2004;19:237–9. ligand–receptor-mediated multicellular signalling in human. Nat Commun. 23. Fujimoto N, He Y, D’Addio M, Tacconi C, Detmar M, Dieterich LC. Single-cell 2015;6:7866. mapping reveals new markers and functions of lymphatic endothelial cells 44. Charrin S, Latil M, Soave S, Polesskaya A, Chrétien F, Boucheix C, Rubinstein in lymph nodes. PLoS Biol. 2020;18:e3000704. E. Normal muscle regeneration requires tight control of muscle cell fusion 24. Swindell WR, Johnston A, Xing X, Little A, Robichaud P, Voorhees JJ, Fisher by tetraspanins CD9 and CD81. Nat Commun. 2013;4:1674. G, Gudjonsson JE. Robust shifts in S100a9 expression with aging: a novel 45. Pawlikowski B, Vogler TO, Gadek K, Olwin BB. Regulation of skeletal muscle mechanism for chronic inflammation. Sci Rep. 2013;3:1215. stem cells by fibroblast growth factors. Dev Dyn. 2017;246:359–67. 46. Pisconti A, Bernet JD, Olwin BB. Syndecans in skeletal muscle development, 25. Kuang S, Chargé SB, Seale P, Huh M, Rudnicki MA. Distinct roles for Pax7 regeneration and homeostasis. Muscles Ligaments Tendons J. 2012;2:1–9. and Pax3 in adult regenerative myogenesis. J Cell Biol. 2006;172:103. 26. Andersen DC, Laborda J, Baladron V, Kassem M, Sheikh SP, Jensen CH. Dual 47. Mylona E, Jones KA, Mills ST, Pavlath GK. CD44 regulates myoblast migration role of delta-like 1 homolog (DLK1) in skeletal muscle development and and differentiation. J Cell Physiol. 2006;209:314–21. adult muscle regeneration. Development. 2013;140:3743. 48. Scimeca M, Bonanno E, Piccirilli E, Baldi J, Mauriello A, Orlandi A, Tancredi V, 27. Waddell JN, Zhang P, Wen Y, Gupta SK, Yevtodiyenko A, Schmidt JV, Bidwell Gasbarra E, Tarantino U. Satellite cells CD44 positive drive muscle CA, Kumar A, Kuang S. Dlk1 is necessary for proper skeletal muscle regeneration in osteoarthritis patients. Stem Cells Int. 2015;2015:469459. development and regeneration. PLoS One. 2010;5:e15055. 49. Burzyn D, Kuswanto W, Kolodin D, Shadrach JL, Cerletti M, Jang Y, Sefik E, Tan TG, Wagers AJ, Benoist C, et al. A special population of regulatory T 28. Zhang L, Uezumi A, Kaji T, Tsujikawa K, Andersen DC, Jensen CH, Fukada S. cells potentiates muscle repair. Cell. 2013;155:1282–95. Expression and Functional Analyses of Dlk1 in Muscle stem cells and 50. Low S, Barnes JL, Zammit PS, Beauchamp JR. Delta-like 4 activates Notch 3 mesenchymal progenitors during muscle regeneration. Int J Mol Sci. 2019; 20:3269. to regulate self-renewal in skeletal muscle stem cells. Stem Cells. 2018;36: 29. Hagan M, Zhou M, Ashraf M, Kim I-M, Su H, Weintraub NL, Tang Y. Long 458–66. noncoding RNAs and their roles in skeletal muscle fate determination. 51. Tajrishi MM, Zheng TS, Burkly LC, Kumar A. The TWEAK-Fn14 pathway: a Noncoding RNA Investig. 2017;1:24. potent regulator of skeletal muscle biology in health and disease. Cytokine 30. De Micheli AJ, Laurilliard EJ, Heinke CL, Ravichandran H, Fraczek P, Soueid- Growth Factor Rev. 2014;25:215–25. Baumgarten S, De Vlaminck I, Elemento O, Cosgrove BD. Single-cell analysis of 52. Pampeno C, Derkatch IL, Meruelo D. Interaction of human laminin receptor the muscle stem cell hierarchy identifies heterotypic communication signals with Sup35, the [PSI+] prion-forming protein from S. cerevisiae: a yeast involved in skeletal muscle regeneration. Cell Rep. 2020;30:3583–3595.e5. model for studies of LamR interactions with amyloidogenic proteins. PLoS 31. Dell’Orso, S., Juan, A.H., Ko, K.-D., Naz, F., Gutierrez-Cruz, G., Feng, X., and One. 2014;9:e86013. Sartorelli, V. (2019). Single-cell analysis of adult skeletal muscle stem cells in 53. Wu Y, Tan X, Liu P, Yang Y, Huang Y, Liu X, Meng X, Yu B, Wu M, Jin H. ITGA6 homeostatic and regenerative conditions. Development dev.174177. and RPSA synergistically promote pancreatic cancer invasion and metastasis via PI3K and MAPK signaling pathways. Exp Cell Res. 2019;379:30–47. 32. Machado L, Esteves de Lima J, Fabre O, Proux C, Legendre R, Szegedi A, Varet H, Ingerslev LR, Barrès R, Relaix F, et al. In situ fixation redefines 54. Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, Baglaenko Y, quiescence and early activation of skeletal muscle stem cells. Cell Rep. 2017; Brenner M, Loh P, Raychaudhuri S. Fast, sensitive and accurate integration 21:1982–93. of single-cell data with harmony. Nat Methods. 2019;16:1289–96. 33. van den Brink SC, Sage F, Vértesy Á, Spanjaard B, Peterson-Maduro J, Baron 55. Barruet E, Garcia SM, Striedinger K, Wu J, Lee S, Byrnes L, Wong A, Xuefeng CS, Robin C, van Oudenaarden A. Single-cell sequencing reveals S, Tamaki S, Brack AS, Pomerantz JH. Functionally heterogeneous human dissociation-induced gene expression in tissue subpopulations. Nat satellite cells identified by single cell RNA sequencing. eLife. 2020;9:e51576. Methods. 2017;14:935–6. 56. Rubenstein AB, Smith GR, Raue U, Begue G, Minchev K, Ruf-Zamojski F, Nair 34. van Velthoven CTJ, de Morree A, Egner IM, Brett JO, Rando TA. VD, Wang X, Zhou L, Zaslavsky E, Trappe TA, Sealfon SC. Single-cell Transcriptional profiling of quiescent muscle stem cells in vivo. Cell Rep. transcriptional profiles of human skeletal muscle. Sci Rep. 2020;10:229. 2017;21:P1994–2004. 57. Riddle ES, Bender EL, Thalacker-Mercer AE. Transcript profile distinguishes variability in human myogenic progenitor cell expansion capacity. Physiol 35. Harmon, B.T., Orkunoglu-Suer, E.F., Adham, K., Larkin, J.S., Gordish-Dressman, Genomics. 2018;50:817–27. H., Clarkson, P.M., Thompson, P.D., Angelopoulos, T.J., Gordon, P.M., Moyna, N.M., et al. (2010). CCL2 and CCR2 variants are associated with skeletal 58. Thiriot A, Perdomo C, Cheng G, Novitzky-Basso I, McArdle S, Kishimoto JK, muscle strength and change in strength with resistance training. J Appl Barreiro O, Mazo I, Triboulet R, Ley K, et al. Differential DARC/ACKR1 Physiol (1985) 109, 1779–1785. expression distinguishes venular from non-venular endothelial cells in 36. Pedersen L, Olsen CH, Pedersen BK, Hojman P. Muscle-derived expression of murine tissues. BMC Biol. 2017;15:45. the chemokine CXCL1 attenuates diet-induced obesity and improves fatty 59. Garcia SM, Tamaki S, Lee S, Wong A, Jose A, Dreux J, Kouklis G, Sbitany H, Seth acid oxidation in the muscle. American Journal of Physiology-Endocrinology R, Knott PD, et al. High-yield purification, preservation, and serial and Metabolism. 2012;302:E831–40. Transplantation of Human Satellite Cells. Stem Cell Reports. 2018;10:1160–74. Micheli et al. Skeletal Muscle (2020) 10:19 Page 13 of 13 60. Xu X, Wilschut KJ, Kouklis G, Tian H, Hesse R, Garland C, Sbitany H, Hansen S, Seth R, Knott PD, Hoffman WY, Pomerantz JH. Human satellite cell transplantation and regeneration from diverse skeletal muscles. Stem Cell Reports. 2015;5:419–34. 61. Sarver DC, Sugg KB, Disser NP, Enselman ERS, Awan TM, Mendias CL. Local cryotherapy minimally impacts the metabolome and transcriptome of human skeletal muscle. Sci Rep. 2017;7. 62. Tarnopolsky MA, Pearce E, Smith K, Lach B. Suction-modified Bergström muscle biopsy technique: Experience with 13,500 procedures. Muscle Nerve. 2011;43:716–25. 63. Spinazzola JM, Gussoni E. Isolation of primary human skeletal muscle cells. Bio-Protocol. 2017;7:e2591. 64. Durinck S, Spellman P, Birney E, Huber W. Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt. Nat Protoc. 2009;4:1184–91. 65. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102:15545. 66. Becht E, McInnes L, Healy J, Dutertre CA, Kwok IWH, Ng LG, Ginhoux F, Newell EW. Dimensionality reduction for visualizing single-cell data using UMAP. Nat Biotechnol. 2018. Publisher’sNote Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Journal

Skeletal Muscle – Springer Journals

Published: Jul 6, 2020

There are no references for this article.