Single cell analysis of adult mouse skeletal muscle stem cells in homeostatic and regenerative conditions

Growth and homeostasis of the skeletal muscle system, a tissue/organ constituting approximately 35% of the human total body mass (Janssen et al., 2000), rely on tissue-resident skeletal muscle stem cells (MuSCs). Following a period of postnatal proliferation leading to a muscle mass increase of several-fold (White et al., 2010), adult MuSCs enter quiescence (Chargé and Rudnicki, 2004). In response to acute or chronic muscle injury, MuSCs are activated, proliferate and differentiate into either newly formed or pre-existing muscle fibers to ensure compensatory muscle regeneration and repair (Brack and Rando, 2012; Comai and Tajbakhsh, 2014). Whereas quiescent and activated MuSCs express the transcription factor Pax7 and accumulate its cognate protein (Seale et al., 2000; Oustanina et al., 2004), transcripts of the myogenic regulators Myf5 and MyoD are present in quiescent MuSCs (Beauchamp et al., 2000; Crist et al., 2012; Liu et al., 2013; Ryall et al., 2015b; de Morree et al., 2017) but their corresponding proteins are detected only after MuSCs become activated (Zammit et al., 2002). MuSCs retain regenerative capabilities when isolated from donors and transplanted into recipient hosts (Sherwood et al., 2004; Collins et al., 2005; Montarras et al., 2005; Kuang et al., 2007). In one study, a single MuSC transplanted into the muscles of mice extensively proliferated giving rise to tens of thousands of cells contributing to muscle fibers. In contrast, transplanted primary myoblasts were ineffective at engrafting and proliferating (Sacco et al., 2008). Failure of MuSCs to efficiently undergo cell division and regenerate is observed during aging (Brack and Rando, 2007; Brack et al., 2007; Kuang et al., 2007; Chakkalakal et al., 2012; Sousa-Victor et al., 2014; Bernet et al., 2014; Cosgrove et al., 2014; Price et al., 2014; Tierney et al., 2014; Sacco and Puri, 2015; García-Prat et al., 2016; Lukjanenko et al., 2016). Moreover, dystrophin deficiency compromises the ability of MuSCs to undergo asymmetric division resulting in inefficient generation of myogenic precursors (Dumont et al., 2015; Feige et al., 2018). The transitions experienced by MuSCs during development and regeneration are regulated by partially overlapping and dynamic transcriptional programs (Buckingham and Relaix, 2007; Sincennes et al., 2016; Chang et al., 2018). Bulk transcriptome analyses of MuSCs have contributed important mechanistic insights by providing static snapshots of averaged gene expression (Pallafacchina et al., 2010; Liu et al., 2013; Sousa-Victor et al., 2014; Yennek et al., 2014; Ryall et al., 2015b; Alonso-Martin et al., 2016; Aguilar et al., 2016; Machado et al., 2017; van Velthoven et al., 2017; Pala et al., 2018). The granularity required to identify distinct cell types, states, and their transcriptional dynamics can, however, be provided by single cell analysis.

Here, we use droplet-based single cell RNA-seq (scRNA-seq) to capture the transcriptional state of thousands of individual MuSCs and primary myoblasts. Using computational approaches, we infer developmental dynamic trajectories of homeostatic and regenerative adult MuSCs and of primary myoblasts, and describe the relative transcriptional changes, including those relative to metabolic pathways, associated with discrete cell state transitions.

We performed scRNA-seq to identify and enumerate cell populations of mouse adult hindlimb skeletal muscles. Muscles of 3-month-old C56BL/6J mice were dissociated into single cells by enzymatic digestion and live cells were isolated by fluorescence-activated cell sorting (FACS) (Fig. 1A). Cells obtained from two mice were sequenced separately. The two samples were tested for similarity and merged for further analysis (Fig. S1A). After quality control, we retained 4414 cells for downstream analysis. On average, we detected 1221 gene transcripts in each individual cell (156,277 mean reads per cell) (Fig. S1B). These cell and transcript numbers are comparable to those recently reported by the Tabula Muris Consortium for droplet scRNA-seq in limb muscle (Tabula Muris Consortium et al., 2018). The sequencing data were analyzed by t-distributed stochastic neighbor embedding (t-SNE) (Satija et al., 2015). Unsupervised clustering of the quality-controlled cells of the two combined mouse samples returned nine cell clusters, transcriptomes of which were related to fibroadipogenic precursors (FAPs), tenocyte-like cells, smooth muscle cells, endothelial cells, adaptive and innate immunity cells (B- and T-lymphocytes, macrophages, monocytes), and a small cluster expressing genes present in mature skeletal muscle (e.g. creatine kinase muscle, Ckm) (Fig. 1B). These findings are similar to those reported in other studies (Verma et al., 2018; Giordani et al., 2018 preprint).

Visualization of the top 20 most variably expressed genes between cell clusters documented distinct transcriptional programs of the nine clusters (Fig. 1C) and expression of known cell lineage-enriched or -specific transcripts was observed in each of the cell clusters (Fig. 1D). Of the total 4414 cells, 217 (5%) revealed a gene expression pattern that could be assigned to MuSCs. Within this cluster, cells were observed to express variable levels of the MuSC markers vascular cell adhesion molecule 1 (Vcam1), Myod (Myod1), Myf5 and desmin (Des) (Chargé and Rudnicki, 2004) (Fig. S1C). Skeletal muscle dissociation and FACS isolation affect a subset of transcripts in a subpopulation of MuSCs (van den Brink et al., 2017). These manipulations, necessary to obtain MuSCs that can be employed in transplantation protocols (Montarras et al., 2005; Sacco et al., 2008), decrease expression of transcripts represented in niche-associated quiescent MuSCs, such as that of the paired box/homeodomain Pax7 (Seale et al., 2000), and promote accumulation of transcripts enriched in activated MuSCs, such as that encoding MyoD (Machado et al., 2017; van Velthoven et al., 2017). Consistent with these observations, t-SNE analysis of the MuSC cluster revealed the presence of two juxtaposed clusters (clusters 1 and 2; Fig. 1E). Using non-parametric kernel density estimation (KDE), we determined the probability density function of cells expressing Pax7 and Myod in MuSC cluster 1 and cluster 2. In cluster 1, KDE identified two cell populations: one with low, the other with high Pax7. Cluster 2 contained only low-expressing Pax7 cells. KDE for MyoD transcripts identified two cell populations in cluster 2: one with lower, the other with higher Myod. Cluster 1 was enriched for low-expressing Myod cells (Fig. 1F). A heatmap illustrating expression of the top 50 most variable genes in the two MuSC subpopulations is shown in Fig. 1G and the corresponding data are reported in Table S1. Overall, scRNA-seq of mononucleated cells obtained from dissociated hindlimb skeletal muscles permitted the identification of distinct cell lineages, including MuSCs.

MuSCs derived from hindlimb muscles of two 3-month-old C56BL/6J mice were prospectively FACS-purified as described (Liu et al., 2015) [VCAM1⁺/CD31 (PECAM1)⁻/CD45 (PTPRC)⁻/Sca1 (Ly6a)⁻] and immediately sequenced (Fig. 2A, Table S1). The two samples were tested for similarity and merged for further analysis (Table S1, MuSCs1 versus MuSCs2 expression sheet; Fig. S2A,B). After quality control, we retained 3081 MuSCs for downstream scRNA-seq analysis. On average, we detected 994 expressed genes in each individual MuSC (Fig. S2C). Using the Chromium platform (10x Genomics), 50,000-70,000 mean reads per cell are generally sufficient to approach saturation, and primary cells with low RNA content and complexity, such as MuSCs, may require less sequencing to achieve saturation reads of 80-90% (https://kb.10xgenomics.com/hc/en-us/articles/115002474263-How-much-sequencing-saturation-should-I-aim-for-; Zhang et al. 2019).

To approximate sequencing saturation, we obtained 236,841 mean reads per cell (Fig. S2C). t-SNE analysis identified two MuSC clusters (Fig. 2B). Transcripts of the cell cycle inhibitor Cdkn1a were present in both clusters whereas those of cyclin E were not detected suggesting that all FACS-isolated MuSCs remained G0-arrested (Fig. 2C). Differentially expressed genes in the two MuSC clusters included growth arrest-specific 1 (Gas1), the Notch target Hes1, Pax7 and the calcitonin receptor (Calcr) in one cluster (MuSCs close-to-quiescence, MuSC cQ) (Fig. 2C,D, Table S1). Genes associated with early activation, such as Myod, Myf5 and ribosomal genes, were instead enriched in the other MuSC cluster (MuSCs early-activation, MuSC eA) (Fig. 2D, Table S1). The MuSC cQ cluster comprised 975 cells (975/3081; 32% of total MuSCs) and the MuSC eA cluster 2108 cells (2108/3081; 68% of total MuSCs). Gene ontology (GO) analysis confirmed that the two MuSC clusters are transcriptionally distinct (Fig. S2D, Table S1). GO terms related to resistance to stress and response to unfolded protein, cell cycle arrest and circadian rhythm were enriched in the MuSC cQ cluster whereas terms indicating activation of ribosome biogenesis, mRNA processing, translation, and protein stabilization were enriched in the MuSC eA cluster (Fig. 2D, Fig. S2D, Table S1). Transcriptome comparison between FACS-isolated MuSCs (ex vivo MuSCs) and in vivo quiescent MuSCs indicated that the MuSC ex vivo transcriptome remains largely reflective of the in vivo transcriptome (van Velthoven et al., 2017). However, another study has reported marked transcriptional differences between MuSCs FACS-isolated from unfixed or paraformaldehyde (PFA)-fixed muscles (Machado et al., 2017). To evaluate the transcriptional state of MuSC cQ and MuSC eA clusters, we compared their respective transcriptomes with that of MuSCs isolated from PFA-fixed muscles to remove genes for which transcription is affected by muscle dissection and FACS isolation (Machado et al., 2017). This analysis revealed that although MuSC eA and MuSCs derived from PFA-fixed muscles shared only 1.8% of transcripts, the percentage increased to 23% in MuSC cQ (Fig. 2E,F).

Myod transcripts are present in quiescent MuSCs but their corresponding protein can be detected only in activated MuSCs. To determine whether FACS-isolated MuSCs cells translated Myod mRNA into the corresponding protein, we captured 217 MuSCs for single cell western blot analysis. Every MuSC was positive for the histone H3 protein. However, MyoD could be detected in only seven out of 217 MuSCs (3.2%) (Fig. 2G). Thus, although tissue and cell manipulations induce transcriptional modifications, these are not immediately followed by MuSC activation as indicated by the paucity of freshly isolated MuSCs expressing detectable levels of the MyoD protein.

We wished to compare scRNA-seq profiles of freshly isolated MuSCs, which have not fully committed to the muscle cell lineage (Chargé and Rudnicki, 2004), with those of committed primary myoblasts (PMs). To this end, PMs isolated from the hindlimb muscles of 3-month-old C56BL/6J mice were cultured in growth medium before scRNA-seq. After quality control, we retained 4429 PMs for downstream analysis. On average, we detected 4933 expressed genes in each individual PM (Fig. S3A). To compare MuSCs and PMs directly, we performed t-SNE analysis and computationally merged their corresponding scRNA-seq datasets (Fig. 3A, Fig. S2E,F). t-SNE analysis segregated the transcriptomes corresponding to MuSCs and PMs into five clusters: two belonging to MuSCs and three to PMs (Fig. 3A, Fig. S2G,H). Analysis of the top 20 genes differentially expressed in each of the five clusters is shown in Fig. 3B. Cells present in both PM clusters (PM1 and PM2) expressed transcripts of cyclins Ccnd1 and Ccnd2 (Fig. 3C). Cells of the PM1 cluster were enriched for cell cycle regulators of the G1/S phase (Table S2). A cluster of PMs containing cells that had activated transcription reflective of differentiation (differentiating myocytes, DM) was characterized by myogenin expression (Fig. 3C). In addition to being negative for Ccnd1 and Ccnd2, the cluster corresponding to MuSC cQ expressed growth arrest-specific 1 (Gas1) and the cell cycle inhibitors Cdkn1a and Cdkn2d, whereas Myf5 was preferentially detected in the cluster corresponding to MuSC eA (Table S2). In addition to terms related to the cell cycle, GO analysis of the PM clusters revealed an enrichment for terms related to metabolic pathways (oxidative and glycolytic enzymes, serine biosynthesis enzymes, enzymes of the pentose phosphate pathway), proteostasis (proteome quality and homeostasis), and activation of translation (Fig. 3D, Table S2).

To obtain temporal resolution of MuSCs and PMs, we employed the pseudotemporal ordering algorithm Monocle2 (Qiu et al., 2017). This algorithm permits single cells to be ordered informatically along trajectories reflecting their progression towards differentiation (Trapnell et al., 2014). Cells located in the initial phase of the trajectory (start of pseudotime) correspond to MuSC cQ, as indicated by expression of Pax7 and the Notch target Hes1 (Fig. 3E). Without an apparent hiatus separating them from the initial phase, cells located further on the trajectory initiate Myf5 transcription and upregulate Myod (MuSC eA). PMs were identified as a distinct group of cells expressing transcripts involved in cell cycle progression (cyclin E and G1/S checkpoint), chromatin remodeling complexes, and regulators of translational initiations (eukaryotic initiation factors, eIFs) (Fig. 3E, Table S2). PMs branched off in two distinct trajectories. One trajectory identified cells progressing towards differentiation, as indicated by activation of myogenin and the muscle structural transcript Tnnt2 (troponin T2) (Fig. 3E). Genes encoding cyclins were downregulated and cell cycle inhibitors Cdkn1a/b upregulated in cells located within this branch (Fig. 3E, DM). Another trajectory stemming out of the PM node identified cells that retained Ccnd1 and Ccnd2 without detectable myogenin expression (Fig. 3E). This branch corresponds to proliferating myoblasts that have not entered the differentiation program and may contain a self-renewal subpopulation of described ‘reserve cells’ (Qiu et al., 2017; Cacchiarelli et al., 2018).

To evaluate how the MuSC transcriptome is impacted by muscle injury, we generated and compared scRNA-seq datasets of homeostatic (uninjured) and regenerating MuSCs. Tibialis anterior (TA) and extensor digitorum longus (EDL) muscles of two 3-month-old C56BL/6J mice were injured with the myotoxic agent notexin and MuSCs were FACS-isolated 60 h after injury (Fig. 4A, Fig. S3B). After quality control, we retained 3550 MuSCs for downstream scRNA-seq analysis. On average, we detected 1336 expressed genes in each individual MuSC (248,914 mean reads per cell; Fig. S3C). The two samples were tested for similarity and merged for further analysis (Fig. S3D,E, Table S2 60 h Inj1 versus 60 Inj2 sheet). t-SNE analysis of injured MuSCs identified three clusters, MuSCs 60 h Cl1, Cl2 and Cl3, distinct from those of homeostatic MuSCs (Fig. 4A, Fig. S3G,H).

GO analysis revealed activation of common transcriptional programs in the three injured MuSC clusters (Fig. S3F, Table S2). Translational regulators were enriched in all injured MuSCs (Eif3c, Eif4ebp1), which had also activated Ccnd1 transcription (Fig. 4B, Fig. S3F, Table S2). Cells located in one of the injured MuSC clusters expressed myogenin and were enriched for the terms ‘muscle contraction’ and ‘myotube differentiation’, indicating that they had entered differentiation (Fig. 4B, Table S2). As noted in bulk transcriptome analyses of injured or culture-activated MuSCs (Pallafacchina et al., 2010; Liu et al., 2013; Ryall et al., 2015b; Pala et al., 2018), expression of genes involved in oxidation-reduction processes, glycolysis, tricarboxylic acid (TCA) cycle, mitochondrial genes, serine biosynthesis and pentose phosphate pathways was also enriched in injured MuSCs (Fig. 4C, Fig. S3F, Table S2). Among them, we found profilin 1 (Pfn1), deletion of which induces a switch from glycolysis to mitochondrial respiration in hematopoietic stem cells (Zheng et al., 2014), and Park7, encoding a redox-sensitive chaperone protecting neurons against oxidative stress and involved in early-onset Parkinson disease, deletion of which promotes glycolysis in skeletal muscle (Shi et al., 2015). Expression of genes known to regulate metabolic pathways and RNA metabolism (Rplp1, Rpl26) (Fig. 4B, Fig. S3F) and expression of nucleophosmin 1 (Npm1) (Fig. S3I), which promotes RNA metabolism and is essential for embryonic development (Colombo et al., 2005) (Fig. S3F), were also increased.

Pseudotemporal ordering aligned individual homeostatic, injured MuSCs, and PMs, along a linear trajectory. Homeostatic MuSCs were located at one end and PMs were found at the opposite end of this trajectory. Injured MuScs were temporally distinct from homeostatic MuSCs and extended towards PMs (Fig. 4D, Fig. S4A). Correct cell ordering was confirmed by plotting genes known to be expressed in the three different cell states (Fig. 4D). Next, we generated heatmaps of pseudotemporal expression patterns to visualize metabolic pathways comprehensively in homeostatic MuScs, injured MuSCs, and PMs (Fig. 4E, Fig. S4B). Expression of mitochondrial genes progressively increased from homeostatic through injured MuSCs to PMs, serving as a molecular correlate of the increased mitochondria observed in activated and cultured MuSCs (Tang and Rando, 2014; Ryall et al., 2015b). Concordantly, genes of the TCA cycle and the electron transport chain were increased in injured MuSCs and in PMs. With the exception of Hk2 (Pala et al., 2018), glycolysis and pentose phosphate pathways were also activated in injured MuSCs and in PMs (Fig. 4E, Table S2). As previously observed, expression of a subset of oxidative genes (Fabp4, Pnpla2, Acaca, Lipa, Acsl1) was elevated in homeostatic MuSCs and decreased in injured MuSCs and in PMs (Ryall et al., 2015b) (Fig. 4E).

Analysis of average gene expression confirmed pseudotemporal expression patterns indicating that transcripts of metabolic enzymes were increased in both injured MuSCs and proliferating PMs compared with homeostatic MuSCs (Fig. 5A). To evaluate the transcriptional dynamics of different metabolic pathways in homeostatic MuSCs, injured MuSCs and PMs, we generated heatmaps of genetic interactions based on the use of pairwise matrices with Euclidean distance. In such matrices, gene proximity indicates co-expression whereas gene distance indicates absence of co-expression. Grouping genes into clusters based on the similarity between their expression profiles allows evaluation of the presence of functionally connected and related modules (Stuart et al., 2003; Nguyen and Lió, 2009). Connectivity mapping indicated that genes of the TCA cycle established an interacting module in homeostatic MuSCs. This module was stabilized, as revealed by increased gene interactions in injured MuSCs, and extended to encompass the majority of TCA genes in PMs (Fig. 5B). Mitochondrial, oxidative phosphorylation and fatty acid oxidation genes steadily increased their connectivity, both in strength and number, from homeostatic MuSCs to PMs (Fig. 5C,D). In a similar manner, electron transport chain (ETC) genes progressively increased their interaction (Fig. S4C). Glycolytic genes established a well-coordinated network in homeostatic MuSCs (Fig. 5E, left). In injured MuSCs, the connectivity of the glycolytic interaction map transitioned to a different conformation characterized by newly acquired intermediate gene interactions (∼0.75 on the color key) of Ldha with all the glycolytic genes and strengthened connectivity of Pkm (Fig. 5E, middle). Gene connectivity further increased in PMs (Fig. 5E, right). Interactions among members of the pentose phosphate pathway also progressively expanded and stabilized in injured MuSCs and PMs (Fig. 5F). Connectivity of genes involved in ribosomal biogenesis steadily increased from MuSCs to PMs and that of protein-folding genes underwent dynamic changes as well. In contrast, the network of genes for structural proteins remained relatively constant (Fig. S4C).

This analysis revealed that although expression of metabolic pathway genes is concordantly orchestrated resulting in a progressive enrichment of all the energy-producing pathways from homeostatic and injured MuSCs to PMs, the dynamics and connectivity of genes within the individual metabolic pathways are distinct and identify discrete cell states.

We have employed scRNA-seq to identify different cell populations of adult mouse skeletal muscle with a special emphasis on MuSCs. FACS-isolated MuSCs could be identified with distinct transcriptional transition states. This phenomenon has been previously documented using bulk RNA-seq approaches (Machado et al., 2017; van Velthoven et al., 2017) and scRNA-seq of 21 and 235 MuSCs, respectively (van den Brink et al., 2017; Cho and Doles, 2017). Here, we extend those findings by analyzing thousands of MuSCs obtained from homeostatic and notexin-injured muscles and comparing the scRNA-seq datasets with those obtained from proliferating committed PMs. t-SNE analysis of FACS-isolated cells identified two MuSC clusters with overlapping but distinct transcriptomes. Indicating MuSC activation, one of the clusters contained cells with higher Myod and lower Pax7 transcript levels whereas the other cluster was composed of cells retaining a transcriptome closer to that of PFA-fixed MuSCs (Machado et al., 2017). Although freshly isolated MuSCs express Myod and Myf5 transcripts, their cognate proteins are not detected until cell activation (Zammit et al., 2002). This regulatory step is controlled by microRNA 489 (Cheung et al., 2012) and microRNA 31 (Crist et al., 2012), and by the RNA-binding proteins of the tristetraprolins (TTP) family, which destabilize Myod mRNA (Hausburg et al., 2015), and Staufen1, which binds to and prevents translation of Myod mRNA in quiescent MuSCs (de Morree et al., 2017). Our single cell western blot analysis of freshly FACS-isolated MuSCs indicated that the MyoD protein could be detected in only a small percentage of MuSCs (∼3%), suggesting that although Myod transcripts are increased during FACS isolation of MuSCs, they are not immediately translated.

By aligning the individual transcriptomes of MuSCs derived from homeostatic and injured muscles, and primary committed myoblasts along a pseudotemporal single trajectory, we could order them, in an unsupervised fashion, and obtain a dynamic description of the transcriptional events characterizing the progression of the individual cells towards differentiation. As expected, the transcriptomes of close-to-quiescence and early-activation clusters present in freshly isolated MuSCs partially overlapped but the existing differences were bioinformatically recognized, thus allowing the positioning of the two cell clusters in close yet distinct space-time positions along the trajectory. Sharp transcriptome differences between homeostatic and injured MuSCs resulted in injured MuSCs occupying a pseudotemporal position intermediate between homeostatic MuSCs and PMs. Alternate trajectories identified PMs that successfully entered terminal differentiation, as indicated by expression of myogenin, muscle structural genes, and cell cycle inhibitors, or that retained a transcriptome indicative of a committed yet undifferentiated state.

Several chaperones were enriched in the MuSC cQ cluster. Carbon dating experiments indicate an average age of 15.9 years for cells residing in human skeletal muscle (Spalding et al., 2005). Thus, chaperones may afford a quality control mechanism by maintaining appropriate protein folding (Hartl et al., 2011). A concomitant enriched expression of genes involved in ‘response to unfolded protein’ suggests tight proteostatic regulation in the MuSC cQ cluster (García-Prat et al., 2017). The MuSC eA cluster was enriched for ribosome biogenesis, RNA splicing, and translational initiation factors. These GO terms were also present in regenerative MuSCs and PMs. Thus, in addition to their post-transcriptional modifications (Zismanov et al., 2016), expression of translational initiation factors was increased in MuSCs and proliferating PMs. Besides serving as building blocks for muscle growth and repair, MuSCs cross-talk with their immediate cellular environment via extracellular vescicles containing RNA and proteins, mediating intercellular communication. During this process, MuSCs secrete exosomes containing miRNAs regulating extracellular matrix (ECM) remodeling (Fry et al., 2017). Coordinated activation of an extensive network of genes controlling extracellular exosomes biogenesis was observed in MuSC eA, regenerative MuSCs, and PMs supporting a potential role for MuSC-secreted vesicles in ECM remodeling that occurs during muscle repair. Interestingly, transcripts of Ccnd1 were detected in MuSCs cQ (Figs 2 and 3), suggesting a role unrelated to cell cycle regulation for cyclin D1.

Metabolic reprogramming is associated with distinct states of embryonic and somatic cells (Ryall et al., 2015a). In MuSCs, transcriptional activation of both glycolytic and mitochondrial genes occurs during cell activation and proliferation (Pallafacchina et al., 2010; Rocheteau et al., 2012; Liu et al., 2013; Rodgers et al., 2014; Ryall et al., 2015b; Pala et al., 2018) and interfering with either glycolysis or peroxisomal fatty acid oxidation bears functional consequences on cell activation and differentiation (Ryall et al., 2015b; Pala et al., 2018). We took advantage of scRNA-seq granularity to investigate metabolic reprogramming in MuSCs and PMs. Regenerating MuSCs are enriched for all the main energy-producing pathways, including glycolysis (Liu et al., 2013; Pala et al., 2018). Bioenergetics profiles of MuSCs freshly isolated from adult mice revealed a positive basal oxygen consumption rate (OCR)/extracellular acidification rate (ECAR) ratio (indicative of oxidative phosphorylation) (Pala et al., 2018), and cultured MuSCs displayed increased ECAR (indicative of glycolysis) compared with freshly isolated MuSCs (Ryall et al., 2015b). In addition, comparison of transcripts present in both FACS-isolated MuSCS and in vivo 4-thiouracil-labeled MuSC transcriptomes revealed an enrichment of transcripts for mitochondrial-associated functions (TCA, electron transport chain, oxidative reduction) (van Velthoven et al., 2017), suggesting oxidative phosphorylation in quiescent MuSCs (Rodgers et al., 2014; Ryall et al., 2015b). By analyzing gene networks based on pairwise matrices of Euclidean distance, we detected overall increase of all metabolic pathways investigated from homeostatic through injured MuSCs to PMs. In injured muscles, increased glycolysis may generate the glycolytic intermediates necessary for the production of new biomass and metabolites utilized by histone- and DNA-modifying enzymes whereas oxidative phosphorylation would generate the ATP required for MuSC proliferation (Ryall et al., 2015a; Pala et al., 2018). Indeed, metabolites generated by different metabolic pathways regulate the activity of chromatin-modifying enzymes known to control different aspects of MuSC biology (McKinnell et al., 2008; Juan et al., 2011; Kawabe et al., 2012; Liu et al., 2013; Ryall et al., 2015b; Boonsanay et al., 2016; Faralli et al., 2016; Scionti et al., 2017; Das et al., 2017; Tosic et al., 2018; Puri and Mercola, 2012). For instance, the oscillation of intermediary metabolites of the NAD biosynthetic pathways may be relevant for circadian gene expression in MuSCs and skeletal muscle (Nakahata et al., 2009; Ryall et al., 2015b; Solanas et al., 2017; Andrews et al., 2010; Lowe et al., 2018). Lastly, perturbations of specific metabolic pathways induce state changes in MuSCs (Ryall et al., 2015b; Pala et al., 2018), revealing a central role of metabolism in regulating cell state transitions. scRNA-seq has permitted the identification of cellular states and a granular definition of the individual cell transcriptional developmental trajectories accounting for the transition of MuSCs from quiescence to proliferation and differentiation. Additional studies will advance our understanding of the regulatory mechanisms governing these transitions at the individual cell level and provide further insights into how metabolic pathways are transcriptionally tuned.

Mice were housed in a pathogen-free facility and all experiments were performed according to the National Institutes of Health (NIH) Animal Care and Use regulations.

Hindlimb muscles from 3-month-old C56BL/6J wild-type mice were dissociated into single cells by enzymatic digestion and live cells were isolated by FACS. Skeletal MuSCs were sorted following described methods (Liu et al., 2015). Briefly, hindlimb muscles from 3-month-old adult wild-type mice were minced and digested with collagenase for 1 h and MuSCs were released from muscle fibers by further digesting the muscle slurry with collagenase/dispase for an additional 30 min. After filtering out the debris, cells were incubated with the following primary antibodies: biotin anti-mouse CD106 (anti-VCAM1, BioLegend 105704; 1:75), PE/Cy7 Streptavidin (BioLegend 405206; 1:75), Pacific Blue anti-mouse Ly-6A/E (anti- Sca1, BioLegend 108120; 1:75), APC anti-mouse CD31 (BioLegend 102510; 1:75) and APC anti-mouse CD45 (BioLegend 103112; 1:75). Satellite cells were sorted by gating VCAM1-positive, Pacific Blue-labeled Sca1-negative, and APC-labeled CD31/CD45-negative cells. SYTOX Green (ThermoFisher Scientific S7020; 1:30,000) was used as a counterstain.

PMs derived from the hindlimbs of 3-month-old C56BL/6J wild-type mice were cultured for four passages in HAM's F10 (Invitrogen) supplemented with 2.5 ng/µl of bFGF (Invitrogen), 10% horse serum (Invitrogen), and penicillin and streptomycin (Invitrogen).

C57/BL6 mice (3-6 months old) were anesthetized using isoflurane (1.5-2.5%) via a vaporizer with anesthetic gas-scavenging equipment. A small incision (3-5 mm) was made on the right-side leg to expose the TA muscle. A volume of 100 µl notexin (2 µM, Accurate Chemical & Scientific) was evenly delivered into TA and EDL muscles using a 29 gauge insulin needle. The incision was closed by wound clips. Buprenorphine (100 µl 0.1 mg/kg diluted in saline) was delivered subcutaneously at the completion of surgery to alleviate the pain. The mice were put on a separate clean bedding cage and were monitored until they became fully awake, then moved back to their original cages.

Cells were collected, washed once with 1×PBS and re-suspended in PBS with 0.04% bovine serum albumin. Cellular suspensions were loaded on a Chromium Instrument (10x Genomics) to generate single cell GEMs. Single cell RNA-seq libraries were prepared using a Chromium Single Cell 3′ Library & Gel Bead Kit v2 (P/N 120737, 10x Genomics). GEM-RT was performed in a C1000 Touch Thermal cycler with 96-Deep Well Reaction Module (Bio-Rad; P/N 1851197): 53°C for 45 min, 85°C for 5 min; held at 4°C. Following retrotranscription, GEMs were broken and the single-strand cDNA was purified with DynaBeads MyOne Silane Beads (Thermo Fisher Scientific; P/N 37002D). cDNA was amplified using the C1000 Touch Thermal cycler with 96-Deep Well Reaction Module: 98°C for 3 min; cycled 12 times: 98°C for 15 s, 67°C for 20 s, and 72°C for 1 min; 72°C for 1 min; held at 4°C. Amplified cDNA product was purified with the SPRIselect Reagent Kit (0.6× SPRI). Indexed sequencing libraries were constructed using the reagents in the Chromium Single Cell 3′ Library & Gel Bead Kit v2, following these steps: (1) end repair and A-tailing; (2) adaptor ligation; (3) post-fragmentation, end repair and A-tailing double size selection cleanup with SPRI-select beads; (4) sample index PCR and cleanup. The barcode sequencing libraries were diluted at 3 nM and sequenced on an Illumina HiSeq3000 using the following read length: 26bp for Read1, 8 bp for I7 Index and 98 bp for Read2.

Sequences from the Chromium platform were de-multiplexed and aligned using CellRanger ver. 2.2.0 from 10x Genomics with default parameters (https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/what-is-cell-ranger). For clustering, filtering, variable gene selection and dimensionality reduction were performed using Seurat ver.2.3.4 (Butler et al., 2018) according to the following workflow: (1) Cells with fewer than 200 detected genes or in which a high percentage of UMIs mapped to mitochondrial genes were excluded from further analysis. (2) Cells with a lower percentage (0.06-0.1%) of UMIs mapped to mitochondrial genes were retained for further analysis. (3) The UMI counts per million were log-normalized for each cell using the natural logarithm (Butler et al., 2018). (4) Variable genes were selected using thresholds (0.05-0.1) after calculating the binned values from log average expression and dispersion for each gene (Butler et al., 2018). (5) The expression level of highly variable genes was scaled along each gene, and we regressed out cell-cell variation in gene expression by the number of detected molecules (5000-8000), and mitochondrial gene expression. (6) The scaled data were projected onto a low-dimensional subspace of principal component analysis (PCA) using dimensional reduction. The number of PCA were decided through the assessment of statistical plots (Butler et al., 2018). (7) Cells were clustered using a graph-based clustering approach optimized by the Louvain algorithm with resolution parameters and visualized using two-dimensional tSNE. (8) To define the clusters’ characteristics, we identified markers for all clusters with low minimum percentage of genes detected and reassigned the names of the clusters on the basis of the distribution of the specific markers. In addition, using the top n markers based on log-fold changing values, we generated an expression heatmap for given cells and genes in the clusters.

Trajectory analysis of scRNA-seq was performed using Monocle v.2.8.0 (Qiu et al., 2017). First, datasets containing information of cell differentiation were created using cell IDs selected through Seurat package. Then, we performed dimensional reduction using the DDRTree method to visualize the dataset, ordered the cells by global gene expression levels, and visualized the trajectory of the dataset. To illustrate the trends of gene expression across pseudotime, we generated heatmaps to represent smooth expression curves of the genes.

To study gene interactions during cell differentiation, we used distance matrices of gene expressions. After selecting genes related to metabolic pathways, we collected the values of their gene expressions in each state and calculated the similarities of the gene expressions using Pearson’s correlation coefficient to generate the corresponding heatmaps (Stuart et al., 2003).

Freshly FACS-isolated MuSCs (1×10⁵) were resuspended in 1 ml of loading buffer (ProteinSimple) and loaded on small scWest chip run (Milo System, ProteinSimple). MuSCs were settled for 15 min, followed by 5 s lysis and 1 min separation (Milo System, ProteinSimple). scWest chip was incubated with sheep anti-histone H3 (Novus Biologicals, NB100-747; 1:40) and donkey anti-rabbit IgG Alexa 555 (Invitrogen, A-31527; 1:40). Mouse anti-MyoD (Novus Biologicals, NBP2-32883) (100 µg/ml) was revealed by donkey anti-mouse IgG Alexa 647 (Invitrogen, A-31571; 1:20). scWest chip was scanned using a microarray scanner and analysis was performed with Scout Software (Milo System, ProteinSimple).

We thank Berendia Jackson, Dr Khooshbu Shah and the ProteinSimple team for assistance with single cell western blot analysis and Jim Simone from the Flow Cytometry Section, NIAMS, for assistance with FACS. This study utilized the high-performance computational capabilities of the Helix System at the NIH, Bethesda, MD, USA (https://hpc.nih.gov/systems/helix.html/).