Three healthy, moderately active individuals (Table 1) performed three 30 s all-out sprint exercises with 2 min of recovery between sprints. Peak power output in the three sprints was on average 8.87 ± 2 W/kg, corresponding to 10.3 metabolic equivalents (METs). Total energy expenditure was 14.5 ± 6.9 KJ.
Cellular yield, composition, and biological processes in different cell populations
The tissue retrieval was 500 ± 150 mg per sample and a total of 39,015 cells from pre- and 24,332 cells from post-exercise were analyzed with a mean read count of 19,507,451 and 19,328 882, respectively. To investigate reproducibility and increase sequence-depth, cells isolated from the third subject were sequenced twice.
Cell-type annotation and reproducibility
The six samples from the three subjects, in addition to the two re-sequenced samples from the third subject, were integrated, anchored, and clustered (Fig. 1a). Inter-sample reproducibility with regards to cell-type composition was assessed by a correlation matrix using inter-sample anchoring scores17. Agreement both between and within subject was consistent, with anchoring scores ranging from 0.6 to 0.7 (Fig. 1b). All subjects contained cells of all cell types detected (Fig. 1c), with comparable proportions between samples (Fig. 1d). Thus, the cell-type annotations and anchoring were highly reproducible across different subjects and biopsies. Sample distribution across clusters is found in Supplementary Figure 1. Following Leiden-based clustering and cell-type analysis, six major cell types were identified: myogenic cells, endothelial cells, pericyte cells, mesenchymal cells, lymphoid cells, and monocyte cells. Cell-type annotation was based on a combination of marker gene expression compared to the CellMatch database, the scCatch pipeline18, and an enrichment test of cell types by using highly expressed/marker genes in each cluster relative to the transcriptomic profiles of cell types provided by the Human Gene Atlas19,20. All populations, except the pericyte cluster, were consistently annotated by ≥2 of the methods and were therefore considered unambiguously annotated with respect to cellular origin. Composition tables for individual samples can be found in Supplementary Data 1.
Comparison with whole-tissue RNA-sequencing
An obvious advantage of scRNA-seq over whole-tissue RNA-sequencing (RNA-seq) is the ability to deconvolute gene expression profiles into their respective cell-type-based compartments. In addition, it is also possible to provide better coverage of gene expression of less abundant cell types, such as immune and stem cells. To address and objectify this notion, genes detected in the present single-cell experiment were compared with the comprehensive profile of the skeletal muscle transcriptome available through the Genotype-Tissue Expression data (GTEx)21. In the 564 skeletal muscle samples from the GTEx study, 5844 unique transcripts were detected readably, while in the current single-cell experiment a total of 1946 unique transcripts were detected. Of the transcripts detected by scRNA-seq, 504 transcripts were not readably detected in the GTEx data, representing 26% of all the transcripts identified by scRNA-seq (Fig. 1e). Among transcripts detected exclusively in the single-cell experiment, gene and cell ontology analysis revealed strong enrichment of genes derived from immune, endothelial, mesenchymal, and smooth muscle cells (Fig. 1f). Among the differentially enriched biological processes associated with these cell populations were “positive regulation of T cell proliferation” (fdr < 0.01), “positive regulation of angiogenesis” (fdr < 0.01), and “positive regulation of endothelial cell proliferation” (fdr < 0.05), indicative of cell type-specific enrichments (Fig. 1f).
Cellular composition and cell-type characteristics
Endothelial cells accounted for 44% of the total number of cells and were distributed in four distinct but neighboring clusters (ESAM+, VWF+/EPS8+, ACKR1+, and FABP4+/NEAT1−) (Fig. 2a). Overall, endothelial cells were characterized by high expression of known endothelial marker genes, including VWF and ESAM (Fig. 2c, d). In comparison with the other cell types, endothelial cells had a significantly higher expression of 165 genes and enriched for 263 gene ontologies such as “regulation vasculature development” (fdr < 0.05), “actin filament organization” (fdr < 0.01), and “response to wounding” (fdr < 0.001) (Fig. 2b). More detailed tables of cell-type and subpopulation specific differential expression together with ontology enrichment can be found in the Supplemental Information.
Mesenchymal cells represented 26% of the total number of cells and showed significantly increased expression of DCN, CFD, and GSN. These cells were further divided into three distinct subpopulations (i.e., DCN+/CFD+, GSN+/LUM+, S100A8+/S100A9+) (Fig. 2a, c, d). The ontologies enriched in the mesenchymal cell population were generic in nature, including “autophagy” (fdr < 0.01), “regulation of anatomical structure morphogenesis” (fdr < 0.01), “oxidative phosphorylation” (fdr < 0.01), “electron transport chain” (fdr < 0.01), and “cell activation” (fdr < 0.05) (Fig. 2b). The S100A8+/S100A9+ subpopulation was considered ambiguously annotated. However, due to the high expression of DCN, VIM, COL6A3, and GSN marker genes the cluster was deemed to be of mesenchymal origin.
Myogenic cells represented 18% of total cells and were further divided into three distinct subpopulations (i.e., PAX7+, TNNI1+, TNNI2+) (Fig. 2a). In addition to PAX7, the first cluster expressed MYF5 and in a smaller frequency MYOD1, indicating this cluster is a mix of undifferentiated satellite cells and early myoblasts. The remaining myogenic subpopulations expressed genes associated with more mature characteristics such as DES, MYL2, and MYOG (Fig. 2c, d), and perhaps terminally differentiated muscle fibers including ACTA1 and MYH6. A subset of these cells also expressed MYOD1. These two subpopulations were distinguished from each other by the expression of slow-twitch (TNNI1) and fast-twitch (TNNI2) troponins. At the ontology level, the myogenic subpopulations exhibited a gene expression profile dominated by key skeletal muscle functions such as “oxidative phosphorylation” (fdr < 0.001) and “mitochondrial respiratory chain” (fdr < 0.001) (Fig. 2b). Consistent with muscle-specific gene expression, the TNNI1 + and TNNI2 + cells also showed enrichment for programs involved in skeletal muscle protein synthesis and degradation (e.g., TRIM63, SYNPO2).
Pericytes accounted for 6% of the total cells, with two distinct clusters (TAGLN+ and NDUFA4L2+) (Fig. 2a). The TAGLN+ cluster was characterized by the high expression of smooth muscle markers such as ACTA2 and pericyte markers (e.g., PDGFR). In addition, this cluster was also annotated as pericyte in origin by scCatch, and therefore pericyte was considered the most probable cellular origin. Significantly enriched ontologies included “mitochondrial electron transport cytochrome c to oxygen” (fdr < 0.001), “nucleoside triphosphate metabolic process” (fdr < 0.001), and “oxidative phosphorylation” (fdr < 0.001) (Fig. 2b).
Two lymphocyte clusters (NAMPT+ and HCST+) were also identified, which accounted for 4% of the total number of cells (Fig. 2a). Considering these clusters presented gene expression profiles characterizing T-, B-, and NK-cells, they were collectively considered as lymphocytes. Genes that distinguished these clusters from the other cell types included CCL5 which is considered as lymphocyte-specific markers (Fig. 2c, d). Ontological enrichment showed “T-cell immunity” (fdr < 0.01), “cell killing” (fdr < 0.05), and “adaptive immune response” (fdr < 0.05).
One monocyte cluster, which accounted for 2% of the total number of cells, was classified as myeloid in origin, based on CD14 expression (Fig. 2a). The cluster differentially expressed genes including CXCL8, AIF, and TYROBP (Fig. 2c, d). HLA-DRA was also highly expressed in this cluster, which is a gene associated with antigen-presenting cells, due to its structural involvement in the formation of human leukocyte antigen (HLA) class II proteins. These CD14+ monocytes enriched for ontologies including “multi organism metabolic process” (fdr < 0.001), “response to corticosterone” (fdr < 0.05), “cellular response to calcium ion” (fdr < 0.05), “response to mineralocorticoid” (fdr < 0.05), “translation initiation” (fdr < 0.001), “Ribosome assembly” (fdr < 0.001), and “Nuclear transcribed mRNA catabolic process nonsense-mediated decay” (fdr < 0.001) (Fig. 2b). Gene-ontology enrichment and marker gene expression for all cell-types can be found in Supplementary Data 2, 3 and 4 and in Supplementary Fig. 2.
Effects of exercise on cellular composition
The effects of high intensity exercise were first examined with respect to cellular composition in skeletal muscle (Fig. 3a). Three hours after exercise, circulating cells increased substantially, with lymphocytes increasing from 4 to 9% (p = 0.05) and monocytes from 2 to 4% (p < 0.05), with a corresponding decrease in the relative contribution of resident cells, i.e., endothelial cells decreased from 44 to 37% and pericytes decreased from 6 to 5% after exercise. The proportion of myogenic cells (18%) remained unchanged after exercise.
Transcriptional response to exercise
The transcriptional response to exercise was assessed by comparing pre- vs. post-exercise for each cell type separately. A total of 874 (535 unique) genes were differentially expressed (fdr < 0.05) across all different cell types. In terms of number of differentially expressed genes, the mesenchymal cells showed the greatest exercise-related response, with 304 genes differentially expressed (Fig. 3b). In ontological terms, the mesenchymal cells enriched for biological functions involved in tissue regeneration and remodeling, such as “regeneration” (fdr < 0.01), “organ regeneration” (fdr < 0.05), and “wound healing” (fdr < 0.05) (Fig. 3c). Genes driving the enrichment for regeneration biological function in the mesenchymal cells included VIM, UBC, GPX4, and AVCRL1 (Fig. 3d). Several genes involved in cytoskeletal reorganization and cell-cycle activation, such as RHOBTB3, TPM1, and RGCC, were also robustly upregulated after exercise in this cell type.
The endothelial cells showed the second greatest response to exercise with a total of 281 differentially expressed genes (Fig. 3b). The main ontological characterization included cell activation and stress reactions, such as “cell cycle G2 M-phase transition” (fdr < 0.05), “energy reserve metabolic response” (fdr < 0.01), and “tissue regeneration” (fdr < 0.05) (Fig. 3c). Differentially expressed genes such as TIMP3, ACTB, UBC, and CALM1 (Fig. 3e) indicate endothelial re-composition and stress response after exercise.
The myogenic cell populations differentially expressed 111 genes (fdr < 0.05) after exercise, including genes involved in differentiation along myogenic lineage such as MYOD1, and MYF6 (Fig. 4). In the undifferentiated PAX7+ cluster gene ontologies related primarily to cellular stress response, such as “negative regulation of cell death” (fdr < 0.001) and “regulation of growth” (fdr < 0.05) (Fig. 3c). The genes driving the stress-related enrichments for the PAX7+ cluster included UBC, HSP90AB1, LMNA, NCL, and SOD2. Muscle-related functions, such as “actin-mediated cell contraction” (fdr < 0.001), and “muscle system process” (fdr < 0.001), were enriched in the more mature clusters (TNNI1+, and TNNI2+) (Fig. 3c) and mainly driven by gene expression such as DES, ACTA, MYL2, and MYOZ1.
The monocyte cell population differentially expressed 120 genes after exercise (Fig. 3b). Ontologically, the monocytes presented a generic response to exercise, with terms such as “amide biosynthetic process” (fdr < 0.001), “protein targeting” (fdr < 0.001), and “peptide metabolic process” (fdr < 0.001) (Fig. 3c). When considering the differential expression of single genes, the monocyte population significantly expressed NFKBIA, CTSS, HLA-C, and SLC25A6 after exercise, suggesting an inflammatory response and cellular stress reactions to exercise (Fig. 3f). Lymphocytes differentially expressed 9 genes after exercise (Fig. 3b). Ontologically, the lymphocyte cells solely enriched for generic terms, suggesting an absent exercise-specific response. The genes regulating the ontological enrichment were overall generic and lacked established biological functions in the literature (e.g., FTH1, IFITM2, CALM1, and HLA-E) (Fig. 3g).
Pericytes differentially expressed 49 genes (Fig. 3b). In terms of ontological enrichment, the pericytes showed a response indicative of cellular activation, stress response, and regeneration, enriching for 38 ontologies such as “response to wounding” (fdr < 0.01), “actin mediated cell contraction” (fdr < 0.05), and “cell cycle G2 M phase transition” (fdr < 0.05) (Fig. 3c). Genes driving the regenerative response included ACTB, TPM1, TIMP3, and CD36. In terms of the greatest exercise-related response, the pericytes differentially expressed ACTB and ACTA2 after exercise (Fig. 3h). A complete list of differentially expressed genes and ontologies across cell-populations can be found in Supplementary Data 5–8.
Comparison of single-cell with whole-tissue sequencing in relation to exercise
We also examined whether the transcriptional responses to exercise among different cell populations using scRNA-seq were consistent with RNA-seq findings. To this end, the transcriptional responses to exercise in each cell type were compared with findings from a recent RNA-seq study22 investigating the same time-points and with a similar exercise protocol as in the current study. Of the 874 transcripts (535 unique) that were significantly regulated in ≥1 cell-type in the current scRNA-seq experiment, 187 transcripts (129 unique) were also regulated by exercise in RNA-seq study (Fig. 3i). There were no differences between cell populations in terms of coverage by RNA-seq, i.e., all cell types shared ~25% of transcripts regulated by exercise with RNA-seq.
Trajectory analysis of myogenic cells
A major advantage of scRNA-seq is the ability to use the global transcriptome of each cell to classify populations of cells of common origin into a continuum of only 1 to 2 dimensions, thereby visualizing and identifying successive, stepwise changes in gene expression from cell to cell. This technique is often referred to as trajectory analysis. Trajectory analysis has proven to be a powerful tool to deconvolute successive transcription-driven cellular processes in developmental biology and stem cell differentiation. Here, we use trajectory analysis to test whether there is evidence of a continuous transition from undifferentiated myogenic cells into increasingly mature myogenic cells (Fig. 4a). Three distinct clusters corresponding to undifferentiated PAX7+ satellite/myoblast cells, TNNI2+ fast-twitch, and TNNI1+ slow-twitch myogenic cells were selected, and two different trajectories with the undifferentiated cluster as a starting point were calculated using principal curve pseudotime analysis. Genetic markers that determined position along the first principal component included PAX7, APOE, IGFBP5, MYF5, and NCAM1 (Fig. 4b), which are canonical myogenic stem/progenitor cell markers. Among the more differentiated cells, ENO3, TNNI2, TNNT3, MYL2, TNNT1, and TNNC1 genes (Fig. 4c) drove the partitioning into distinct clusters along the second principal component. The transition of undifferentiated PAX7+ satellite/myoblast cells towards a higher degree of differentiation was proportional to expression levels of marker genes along the given trajectories. Accordingly, the expression of TNNC1 and TNNC2 increased proportionally to pseudotime, in slow- and fast-twitch myogenic cells, respectively (r2 = 0.28 and r2 = 0.32, Fig. 4d, e). A corresponding decrease in the expression of MYOD1 (r2 = −0.39 and r2 = 0.13) and MYF6 (r2 = −0.29 and r2 = −0.22) was observed along the trajectories for both fast- and slow-twitch cells (Fig. 4d, e).
Furthermore, we examined the effect of exercise on the transition from PAX7+ undifferentiated satellite/myoblast cells to fast- and slow-twitch expressing cells and the overall transcriptional effect of exercise in these subpopulations. Three hours after a single bout of exercise, there was a small (Δslow-twitch = 5.4%, p < 0.001; Δfast-twitch = 9.0%, p < 0.001) but statistically significant incremental shift in pseudotime toward a higher degree of differentiation in both the slow- and fast-twitch cell populations (Fig. 5a, b). Apart from the effects on pseudotime, there was a common transcriptional response to a single bout of exercise for 24 genes in all three myogenic subpopulations. These genes included mostly mitochondrial and ribosomal genes. More genes were regulated by exercise in the undifferentiated PAX7+ (255 genes) and slow-twitch (138 genes) subpopulations, compared to the fast-twitch subpopulation (51 genes) (Fig. 5c). Genes regulated by exercise in the undifferentiated PAX7+ myogenic cells included NEAT1, NNMT, CXXC5, MT2A, and SQSTM1 (Fig. 5d). The slow-twitch TNNI1+ cells responded to exercise by regulating genes involved in “muscle organ development” (fdr < 0.001) including MYLPF, ACTA1, MYL2, and MYH7 (Fig. 5d). In the fast-twitch TNNI2+ cells, the exercise response included upregulation of changes in MYBPC1, MYBCP2 (“muscle contraction”, fdr < 0.001), TXNIP, UBC, and OPTN (Fig. 5d). Gene-ontology enrichment in the myogenic cells following exercise can be found in Supplementary Data 9.
In vitro validation of myogenic subpopulations
The observation of a regulated expression of contractile elements such as TNNC1 and TNNC2 in myogenic cells, and that such transcriptional upregulation is associated with the initiation of differentiation towards a more mature myofiber-like phenotype was validated in vitro. Primary myoblasts were isolated from human muscle biopsy and kept in proliferation media until confluence. Differentiation towards myotube formation was initiated according to standard protocols and the expression of key-marker genes from the myoblast subpopulations identified as undergoing differentiation were analyzed with qPCR on day 0, day 4, and day 9 of the differentiation process (Fig. 5e). Slow-twitch troponin (TNNC1) increased from 9.3 ± 2 a.u on day 0 to 331 ± 150 after 4 days of differentiation (p < 0.001). It remained elevated after 9 days of differentiation (117 ± 20) (p < 0.001). Gene-expression of fast-twitch troponin (TNNC2) was 2.2 ± 1 a.u on day 0, increased to 123 ± 10 a.u after 4 days of differentiation (p < 0.001), and remained elevated after 9 days of differentiation (155 ± 30) (p < 0.001).
Finally, we utilized a publicly available microarray experiment conducted in a mouse myoblast cell-line (C2C12-cells) evaluation gene-expression in cells undergoing differentiation (n = 3) in relation to cells in proliferation-media (n = 3) (Fig. 5f). TNNC1 and TNNC2 were elevated with a log2FC of 5.7 and 7.7 respectively (fdr < 0.001) in differentiating versus proliferating C2C12-cells.