scRNA-seq of skeletal muscle regeneration
To comprehensively evaluate skeletal muscle homeostasis and regeneration throughout aging, we performed scRNA-seq on 65 mouse skeletal muscle samples with the 10x Chromium v2 and v3 platforms. Muscle damage was induced in young (4–7 months old (mo)), old (20 mo) and geriatric (26 mo) C57BL/6J mice by injecting the TA muscles with notexin. The injured and uninjured muscles were collected at six time points (days 0, 1, 2, 3.5, 5 and 7; Fig. 1a,b). Together, these 65 scRNA-seq samples, newly reported here and from our two prior reports11,26, contained a total of 365,710 cell barcodes before quality control and filtering (Extended Data Fig. 1a and Supplementary Table 1). All samples were processed by aligning sequencing reads to the mm10 mouse reference genome, removing ambient RNA signatures with SoupX27, removing low-quality cells and identifying and removing doublets with DoubletFinder28 (Fig. 1c and Extended Data Fig. 1). The samples were integrated with Harmony29 to correct for batch effects (Fig. 1c). The final dataset contained 273,923 cells. We observed that the cell number was relatively evenly distributed across the three age groups (Fig. 1d,g), and within each age group, the number of cells from each time point was relatively consistent (Fig. 1e).
a,b, Overview of experimental design (a) and summary of replicates (b). 3′ scRNA-seq (10× Chromium v2 and v3) was performed on dissociated TA muscles from young (4–7 mo), old (20 mo) and geriatric (26 mo) mice (both sexes) 0–7 days post injury (dpi) with notexin. n = 3–4 replicates per age and dpi. c, Processing workflow. Each scRNA-seq sample was aligned to the mm10 mouse reference genome, ambient RNA was removed by SoupX, low-quality cells were identified and removed, and doublets were identified and removed. All samples were then integrated with Harmony, resulting in a final dataset containing 273,923 cells from 65 samples. See Extended Data Fig. 1 and Supplementary Table 1 for additional detail. QC, quality control; UMI, unique molecular identifier. d, Fraction of cells from each age group. e, Fraction of cells from each dpi within each age group. f,g, Uniform manifold approximation and projection (UMAP) representations of the final dataset. Cells colored by manually assigned cell-type IDs based on the expression of hallmark skeletal muscle genes (see Extended Data Figs. 2 and 3) (f). Cells are colored by age group, with all other cells in gray (g). DC, dendritic cell.
Multistep identification of cell types by clustering
After data integration, shared nearest-neighbor (SNN) clustering was performed, and canonical marker genes were used to manually identify cell types. The initial clustering resulted in 24 clusters that each received a unique cell-type annotation. Of these, nine were myeloid cell clusters (monocytes, macrophages, dendritic cells and neutrophils) that exhibited similar expression profiles. To further clarify the myeloid cells found in the final dataset, the nine myeloid cell clusters were subset out, re-clustered and re-embedded (Extended Data Fig. 2a–c). This resulted in 15 clusters that each received a unique cell-type annotation based on known myeloid markers (Extended Data Fig. 2c,d). Although this re-clustering did not further clarify the monocyte and macrophage annotations, it did help to identify more specific dendritic cell and T cell subtypes (Extended Data Fig. 2c,d). The dendritic cell and T cell-subtype annotations were transferred from the myeloid subset back to the final dataset based on the cell barcode. The monocyte and macrophage annotations were not changed based on the myeloid subset. With the additional myeloid subset annotations, we identified 29 distinct cell types in the final dataset (Fig. 1f and Extended Data Fig. 3). Compared to prior scRNA-seq and single-nucleus RNA-sequencing skeletal muscle studies, we identified similar broad cell types and we identified more specific endothelial, FAP and immune cell subsets7,9,10,11,12,13,14,15,16,17,18,19,20,21,30.
Lymphoid cell types were 5 of the 17 immune cell clusters. We identified a rare B cell cluster (0.47% of the final dataset) that expressed Cd19, Cd22 and Ms4a1 (Extended Data Fig. 3a). We identified a natural killer (NK) cell cluster that expressed Nkg7, Gzma, Klra4 and Klre1 but, importantly, did not express T cell markers (Extended Data Fig. 3a). We identified three T cell clusters, all of which expressed Cd8a and Cd8b1. One of the T cell populations more strongly expressed Cd4, while the other two T cell populations expressed Cd3e. One of the Cd3e+ clusters also strongly expressed the cycling markers Cdk1 and Hmgb2 and were thus identified as cycling Cd3e+ T cells (Extended Data Fig. 3a). The other Cd3e+ T cell population appears to be non-cycling. The B cells, NK cells and the non-cycling Cd3e+ T cells were identified as unique clusters in the initial clustering of the final dataset. The Cd4+ and cycling Cd3e+ T cells were identified when we subset and re-clustered the myeloid populations (Extended Data Figs. 2c,d and 3a).
Although we incubated the single-cell suspension in erythrocyte lysis buffer, we did see a small (0.42% of the final dataset) erythrocyte cluster that uniquely and strongly expressed a variety of hemoglobin genes, including Hba-a1 and Hbb-bs (Extended Data Fig. 3a). Erythrocytes are not a native cell type in skeletal muscle, so we have excluded them from the cell-type dynamics analysis.
We identified three FAP populations (adipogenic, pro-remodeling and stem), a tenocytes population and a Schwann and neural/glial cell population. All three of the FAP populations expressed Pdgfra and Col3a1 (Extended Data Fig. 3b). The adipogenic FAPs also expressed Adam12, Bmp5, Myoc, Col1a1, Dcn, Mmp2 and Apod. The pro-remodeling FAPs uniquely expressed cycling genes like Cdk1 and Tyms in addition to other FAP markers like Tnfaip6, Il33, Adam12, Bgn and Hdlbp. The stem FAPs also expressed Igfbp5, Dpp4, Cd34, Gsn and Mmp2. The tenocyte population expressed some FAP markers like Col1a1, Dcn and Apod, and they expressed tenocyte-specific markers Tnmd and Scx. The Schwann and neural/glial cell cluster expressed Ptn and Mpz (Extended Data Fig. 3b).
We identified a pericyte and smooth muscle cell cluster that expressed the pericyte-specific genes Rgs5 and Acta2, Myl9 and Myh11 (Extended Data Fig. 3c)31. Four endothelial clusters were identified, and they shared strong expression of Cdh5 and Pecam1 (Extended Data Fig. 3c). The arterial endothelial cells uniquely expressed Alpl and Hey1, the capillary endothelial cells strongly expressed Lpl, and the venous endothelial cells expressed Vwf, Hif1a, Icam1, Lrg1 and Aplnr (Extended Data Fig. 3c). The fourth endothelial cluster expressed endothelial markers like Cdh5 and Pecam1 and myeloid markers like S100a8/S100a9, Csf1 and Itgam (Extended Data Fig. 3c). This cluster was small (0.073% of the final dataset) and was made up of cells from multiple replicates from the three ages at 0–2 dpi (Extended Data Fig. 5e). Because this cluster was not unique to a single replicate, or single age, or single time point, we have maintained it in our analysis.
Comparison of cell-type dynamics across mouse lifespan
For each time point independent of age, we calculated the percentage of cells from each cell type. This calculation represents the cell-type abundance within the scRNA-seq dataset, not within skeletal muscle. Multinucleated myofibers are not well captured in the single-cell isolation protocol, so the abundance of myonuclei is underrepresented within the scRNA-seq dataset. Additionally, scRNA-seq libraries are limited by cell capture number, so the cell-type abundances reported are conditional on the prevalence of other cell types. In agreement with previous scRNA-seq findings, uninjured skeletal muscle (day 0) was mainly composed of endothelial cells (49.2%), FAPs (20.4%) and myonuclei (15.1%), and there were small populations of immune cells (2.7%) and MuSCs and progenitors (2.0%; Extended Data Fig. 4a and Supplementary Table 1a)7,11. Of the endothelial cells at day 0, the most prominent subtype was the capillary endothelial cell (42.4%; Fig. 2k, Extended Data Fig. 4b and Supplementary Table 1b). The FAPs present at day 0 were mainly adipogenic FAPs (15.6%; Fig. 2k, Extended Data Fig. 4b and Supplementary Table 1b). Of the immune cells present at day 0, there were a variety of monocytes and macrophages (47.2%), B cells (14.6%), non-cycling Cd3e+ T cells (11.6%), Cd209a+ dendritic cells (10.6%) and neutrophils (9.3%; Extended Data Fig. 4d and Supplementary Table 1d)7.
a–i, Line plots showing the relative abundance of Ccr2+ monocytes/macrophages (a), Mrc1+ monocytes/macrophages (b), Cx3cr1+ monocytes/macrophages (c), cycling monocytes/macrophages (d), patrolling monocytes/macrophages (e), NK cells (f), non-cycling T cells (g), B cells (h), and MuSCs and progenitors (i) as a fraction of total cells from 0 to 7 dpi. For each sample, the number of cells of the reported type was divided by the total number of cells (excluding erythrocytes). Points are each sample (n = 3–4). Ribbon is the standard deviation. Statistical significance of age-specific cell-type dynamics was evaluated using a one-way analysis of variance (ANOVA) and false discovery rate (FDR)-corrected P values are reported (Extended Data Fig. 6). j, Statistical significance of age-specific cell-type dynamics differences as reported by FDR-corrected P values from a one-way ANOVA. Red line denotes the FDR-corrected P = 0.05 threshold. k, Stacked bar plot of the fraction of cell types across all dpi values. All immune cells were grouped into ‘Immune cells (17 cell types)’. l, Scatterplots of the fraction of all T cells (three cell types) at 0 and 5 dpi from the scRNA-seq data. For each sample, the number of all T cells (three cell types) was divided by the total number of cells (excluding erythrocytes). Points are each sample (n = 3–4). Line is the mean for each age group. Significance was evaluated using a two-sided Welch’s t-test. m, Immunohistochemical analysis of CD3+ T cells at 0 and 5 dpi in TA muscles from young, old and geriatric mice. These are representative images from n = 2 samples per age and time point where four images were taken for each sample. Scale bar, 50 µm. Arrow denotes CD3+ T cell. n, Flow cytometric analysis gating strategy of TCRβ+ T cells. o, Scatterplots of the fraction of TCRβ+ cells of all CD45+ cells at 0 and 5 dpi. Points are each sample (n = 2–7). Line is the mean for each age group. Significance was evaluated using a two-sided Welch’s t-test. p, Violin plot of the T cell exhaustion score in all T cells by age at 5 dpi. Significance was evaluated using a two-sided Mann–Whitney U test. q, Dot plot showing the expression frequency and average expression level of select T cell, cycling, inhibitory and T cell exhaustion markers by age at 5 dpi. Prog TEX, progenitor exhausted T cell; Trans TEX, transitory exhausted T cell; Term TEX, terminal exhausted T cell.
Following injury (days 1, 2, 3.5 and 5), the most abundant general cell type was the immune cells (60.9%, 79.7%, 75.2% and 63.2%, respectively; Fig. 2k, Extended Data Fig. 4b and Supplementary Table 1b). As in previous studies7,11, non-immune cells like endothelial cells, FAPs and myonuclei were present following injury, but at transiently lower relative abundances. As expected, the most abundant immune cells immediately after injury (day 1) were neutrophils (32.0%), Ccr2+ monocytes/macrophages (19.7%), and Ctsa+ patrolling monocytes/macrophages (13.7%; Extended Data Fig. 4c and Supplementary Table 1c)7,11. Immediately following injury (days 1 and 2), there was a more pro-inflammatory environment, as evident by the abundance of Ccr2+ monocytes/macrophages. This was followed by a shift at days 3.5 and 5 to a more anti-inflammatory cell population, as evident by the peak in abundance of Cx3cr1+ monocytes/macrophages at day 5 (28.7%; Extended Data Fig. 4c and Supplementary Table 1c)7,11. By day 7, the cell-type abundances were returning to the abundances observed at day 0, but there was still a substantial immune cell population (32.1%; Extended Data Fig. 4a and Supplementary Table 1a). The immune population at day 7 mainly consists of Cx3cr1+ monocytes/macrophages (24.1%), Ctsa+ patrolling monocytes/macrophages (13.2%) and non-cycling Cd3e+ T cells (17.8%; Extended Data Fig. 4d and Supplementary Table 1d).
We next conducted an analysis of cell-type dynamics across age groups by comparing the abundance of a given cell type over the entire injury time course. Of the 28 cell types, 8 were identified as having significantly different cell-type dynamics across the three ages using a nonlinear modeling approach with multiple-hypotheses correction (Fig. 2j and Supplementary Table 2). In response to injury, we first observed an increase in neutrophils that peaked in abundance at day 1 and returned to day 0 levels by day 3.5 (Extended Data Fig. 5n). We also observed a peak in the abundance of the Ccr2+ monocytes/macrophages at days 1 and 2 (Fig. 2a), while the cycling monocytes/macrophages peaked in abundance at day 3.5 (Fig. 2d) and the Cx3cr1+ monocytes/macrophages peaked in abundance at days 3.5 and 5 (Fig. 2c). The Mrc1+ monocytes/macrophages responded early to injury (day 2) and remained high through day 5 (Fig. 2b). The geriatric Mrc1+ monocytes/macrophages maintained a higher abundance from day 2 to day 7 compared to the young and old cells (Fig. 2b). The young and old Ctsa+ patrolling monocytes/macrophages peaked in abundance at day 2, while the geriatric Ctsa+ patrolling monocytes/macrophages peaked in abundance at day 3.5 (Fig. 2e). Additionally, some lymphoid cell types like NK cells, Cd3e+ non-cycling T cells and B cells started to increase in abundance at day 2, day 3.5 and day 5, respectively (Fig. 2f–h). We observed a similar pattern when looking at all three T cell populations combined (Extended Data Fig. 5t). The geriatric NK cells did not increase in abundance within the 7-day time course, unlike the young and old cells (Fig. 2f).
When looking at all three T cell populations combined, we detected very few T cells and no age-specific differences in abundance at day 0. However, at day 5, we observed a higher abundance of T cells in old samples compared to both the young and the geriatric samples (Fig. 2l). To confirm this, we performed immunohistochemistry on sectioned TA muscles and observed that CD3+ T cells are detected more abundantly at day 5 compared to day 0 in the TAs of young, old and geriatric mice (Fig. 2m). We further used flow cytometry to quantify CD45+CD11c−CD11b−TCRβ+ T cells at days 0 and 5 from dissociated TA muscles of young, old and geriatric mice (Fig. 2n). We detected a low abundance of TCRβ+ T cells of all CD45+ hematopoietic cells at day 0, but still detected a significantly higher T cell abundance in old mice compared to young (Fig. 2o; Welch’s t-test, *P value = 0.001). Further, we observed an increase in T cell abundance from day 0 to day 5 and a significantly higher T cell abundance in old mice compared to young mice at 5 dpi (Fig. 2o; Welch’s t-test, *P value = 0.02). Together, these results suggest the abundance of the T cell pool is elevated specifically in older muscles (20 months of age), but not preserved in geriatric ages. To investigate the differences between T cell populations by age at 5 dpi, we calculated a T cell exhaustion score and looked at the expression of exhausted T cell markers (Fig. 2p,q). We found that the old and geriatric T cells have a higher T cell exhaustion score compared to the young (Fig. 2p; Mann–Whitney U test, *P = 2.2 × 10−16 and *P = 1.1 × 10−12, respectively). Genes that have been shown to mark progenitors of exhausted T cells were more highly expressed in old T cells (Bach2, Irf4, Batf, Slamf6, Tcf7). Additionally, genes that have been shown to mark transitory (Tbx21, Runx1/Runx2/Runx3) and terminally (Eomes, Tox, Nr4a1/Nr4a2/Nr4a3, Cd38, Entpd1) exhausted T cells are more highly expressed in both the old and geriatric T cells (Fig. 2q)32,33. Further, we observed that T cell inhibitory genes (Ctla4, Pdcd1, Lag3, Havcr2) were also elevated in the geriatric T cells. We reason that the increase in T cell abundance observed only in old muscles but not in geriatric muscles can be explained by T cells shifting to an inhibited, exhausted state in older muscles and, therefore, being unable to proliferate productively by geriatric ages.
Independent of age, we detected few MuSCs and progenitors immediately following injury (days 1 and 2; 2.1% and 1.1%, respectively) and the abundance of MuSCs and progenitors peaked at day 5 (9.7%; Extended Data Fig. 4a and Supplementary Table 1a). This was in agreement with previous studies7,11. The peak in abundance of MuSCs and progenitors did vary by age with the young cells peaking at day 5, the old cells peaking at day 3.5 and the geriatric cells peaking at day 7 (Fig. 2i). This difference in peak abundance did not result in a statistically significant difference in the MuSCs and progenitors dynamics by age, but it did demonstrate a delayed response by the geriatric MuSCs and progenitors. Independent of age, we detected the most myonuclei at days 0 and 7 (15.1% and 10.6%, respectively; Extended Data Fig. 4a and Supplementary Table 1a). The myonuclei dynamics were very similar between the three ages, but we detected more myonuclei in old and geriatric samples at both days 0 and 7 compared to young (Fig. 2j and Extended Data Fig. 5m).
Senescence scoring based on single-cell transcriptomics
Next, we sought to investigate age-specific differences in senescence within skeletal muscle regeneration. Hallmarks of mammalian aging include stem cell exhaustion, altered cellular communication and cellular senescence34. Identifying senescent cells in scRNA-seq data is challenging because the markers traditionally used to identify senescent cells are lowly expressed, expressed in select cell types in single-cell data and/or assayed in terms of cellular localization and enzymatic function35. For example, senescent cells are commonly identified by persistent expression of cell cycle regulators p16 (Cdkn2a), p21 (Cdkn1a), p53 and/or Rb36. Senescent cells are also marked by the senescence-associated secretory phenotype (SASP), which includes pro-inflammatory cytokines and chemokines, growth modulators, angiogenic factors and matrix metalloproteinases (for example, Mmp3)37. To examine individual gene signatures of senescence, we quantified the abundance of Cdkn2a (encodes p16), Cdkn1a (encodes p21), Mmp3 (a senescence-associated matrix metalloproteinase) and Glb1 (encodes senescence-associated β-galactosidase) across all cell types and ages (Fig. 3a).
a, Dot plot showing the expression frequency and average expression level of select senescence markers. b–d, Scatterplots of the fraction of Cdkn2a+ cells within young, old and geriatric MuSCs and progenitors (n = 20, 21 and 24, respectively) (b), young, old and geriatric T cells (cycling; Cd3e+) (n = 15, 16, 19, respectively) (c), and young, old and geriatric FAPs (adipogenic; n = 20, 21 and 24, respectively) (d). Points are the fraction for each sample. Horizontal line is the median for each age group. Significance was evaluated using a two-sided Mann–Whitney U test. e, Line plot of the fraction of MuSCs and progenitors that coexpress Cdkn2a and Cdkn1a from 0 to 7 dpi. Points are the fraction for each sample (n = 4). Ribbon is the standard deviation. Statistical significance of age-specific dynamics was evaluated using a one-way ANOVA, and the FDR-corrected P value is reported. f, Scatterplot of the fraction of MuSCs and progenitors at 3.5 dpi that coexpress Cdkn2a and Cdkn1a by age. Points are the fraction for each sample (n = 4). Horizontal line is the median for each age group. Significance was evaluated using a two-sided Student’s t-test. g, Scatterplot of the normalized expression level of Cdkn2a and Cdkn1a in MuSCs and progenitors at 3.5 dpi. The density is shown on the top and to the right. Blue line represents the linear trend. Ribbon is the confidence interval. The inset contains the Pearson correlation coefficient and its statistical significance, calculated with the R function cor.test(alternative = ‘two.sided’). h, Significantly upregulated or downregulated GO terms between Cdkn2a+ and Cdkn1a+ (double positive) and all other (other) MuSCs and progenitors at 3.5 dpi. The normalized enrichment score and the FDR-corrected q values were obtained from GSEA. i,j, Feature plots of the final dataset with cells colored by the rescaled two-way FBR Sen Score (i) and the rescaled one-way FBR Sen Score (j). Cells are randomly plotted. k,l, Violin plots of the two-way (k) and one-way (l) FBR Sen Scores in select cell types. m–p, ROC curves based on the coexpression of Cdkn2a and Cdkn1a for the six two-way in vitro Sen Scores (m), the four two-way in vivo Sen Scores (n), the six one-way Sen Scores (o) and select two-way Sen Scores split by age (p). The AUC is reported for each ROC curve.
Expression of these genes depended more on cell type than age, when considering all time points and samples together (Fig. 3a). Given that Cdkn1a was widely expressed across many cell types, we focused on the common senescence hallmark gene Cdkn2a. Cdkn2a transcripts were rarely detected, in agreement with previous observations in the Tabular Muris Senis project8, and primarily observed in the MuSCs and progenitors and cycling T cells. We observed age-associated changes in the relative abundance of Cdkn2a+ cells within any given cell type infrequently significant (Fig. 3b–d). There was no significant difference in the fraction of Cdkn2a+ MuSCs and progenitors by age (Fig. 3b). We did observe a significant difference in the fraction of Cdkn2a+ cycling T cells between the young and old ages (Fig. 3c; Mann–Whitney U test, *P value = 0.04) and in the fraction of Cdkn2a+ adipogenic FAPs between the young and geriatric ages (Fig. 3d; Mann–Whitney U test, *P value = 0.02*).
Within the MuSCs and progenitors, we observed an increase in the fraction of cells that coexpressed Cdkn2a and Cdkn1a from day 0 to day 3.5, after which it returned to near day 0 levels (Fig. 3e). Although there was no significant difference by age in the abundance of these double-positive cells across all time points, there was a significant difference at day 3.5 between the young and old ages (Student’s t-test, *P value = 0.03) and the young and geriatric ages (Student’s t-test, *P value = 0.01; Fig. 3e,f). Although Cdkn2a and Cdkn1a are both cell cycle inhibitors and senescence markers, their expression was not correlated in MuSCs and progenitors at day 3.5 on the individual cell level in this dataset, possibly due to transcript detection dropout (Fig. 3g). We considered the double-positive Cdkn2a+ and Cdkn1a+ cells as candidate senescent MuSCs and progenitors and performed gene-set enrichment analysis (GSEA) on day 3.5 at their peak abundance. GSEA found that double-positive MuSCs and progenitors are enriched for gene sets associated with muscle weakness and various mitosis-related processes but diminished in muscle contraction and cytoskeletal processes (Fig. 3h). Collectively, these GSEA results suggest that double-positive MuSCs and progenitors have signatures of dysregulated muscle function and stalled cell cycle-related gene expression.
We then used the Cdkn2a+ Cdkn1a+ MuSCs and progenitors as a candidate cell population to evaluate broader senescence signatures at the single-cell level. We tested two senescence scoring methods38,39 and 14 senescence-signature (SenSig) gene lists38,40,41,42,43,44,45 (Supplementary Table 2). We refer to the first method as the two-way senescence score (Sen Score) because it calculates a score based on a list of upregulated and downregulated genes. Within this method, we tested nine gene lists that were experimentally derived38,40,44. We refer to the second method as the one-way Sen Score because it calculates a score based on a list of upregulated genes. Within this method, we tested six gene lists, two of which were taken from Gene Ontology (GO) databases42,43, one was curated in another report41 and three were experimentally derived38,45. We refer to the Methods, Fig. 4 and Supplementary Table 2 for more details on these two methods and the gene lists. We note that two of these gene lists are derived from bulk RNA-seq differential expression analyses of p16+ and p16− cells selected based on transgenic reporter status. One gene list (foreign body response, or FBR) was generated by Cherry et al. from p16+ versus p16− CD29+ cells isolated from an FBR-driven skeletal muscle fibrosis model in adult p16-CreERT2;Ai14 reporter mice38. A second gene list (‘aged chondrocytes’) was generated from p16+ versus p16− Aggrecan+ chondrocytes isolated from 20-mo p16-tdTom;Aggrecan-CreERT2;Ai6 mice (B.O. Diekman, personal communication)46,47,48. To compare how the choice of gene list and method impacted senescence scoring across different cell types, we examined the two-way and one-way FBR Sen Scores in the final scRNA-seq dataset. The two-way FBR Sen Scores were more consistently low across most cell clusters and exhibited high scores most notably in the FAPs and MuSCs and progenitors (Fig. 3i,k). The one-way FBR Sen Scores had broader distribution, with many more cell types exhibiting high scores, including FAPs, MuSCs and progenitors, endothelial cells, pericytes and smooth muscle cells and monocytes/macrophages (Fig. 3j,l). Whereas the two-way FBR Sen Scores were mean-centered around zero for each cell-type cluster due to their z-scored counts (Fig. 3k), one-way FBR Sen Scores had varied cell-type averages (Fig. 3l). These differences complicate establishing a threshold for senescence between cell types in the one-way Sen Scores.
a–f, ROC curves based on the coexpression of Cdkn2a and Cdkn1a for the nine two-way Sen Scores (a–c) and the six one-way Sen Scores (d–f) split into young (a and d), old (b and e) and geriatric (c and f) MuSCs and progenitors. g–j, ROC curves based on the coexpression of Cdkn2a and Cdkn1a for the nine two-way Sen Scores (g and h) and the six one-way Sen Scores (i and j) split into G1 (g and i) and non-G1 (S, G2, and M) (h and j) MuSCs and progenitors. k, Table of the AUC for each ROC curve. l–n, Venn diagrams of the unique and shared upregulated genes found in the two-way in vitro FBR SenSig gene lists (l), the two-way in vivo FBR SenSig gene lists (m) and the one-way in vivo FBR SenSig gene lists (n).
We established a scoring approach calibrated for sensitivity and specificity in discriminating Cdkn2a+ Cdkn1a+ MuSCs and progenitors across all ages and time points and present the results in receiver operating characteristic (ROC) curves with performance reported using an area under the curve (AUC) metric (Fig. 3m,p). Between the two senescence scoring methods and the 14 SenSig gene sets, the one-way FBR method performed the best of any experimentally derived approach (AUC = 0.86) and was comparable to the ontology-curated one-way GO:SASP approach (AUC = 0.88; Fig. 3m–o). Notably, it performed far better than the recently described SenMayo list while using the same single-sample GSEA method (AUC = 0.63). Moreover, the one-way FBR Sen Score accurately discriminated double-positive senescent-like MuSCs and progenitors across all three ages (AUCs = 0.85–0.86), suggesting it captures common features of senescence irrespective of age (Fig. 4d–f,k). We concluded that the one-way FBR method was able to accurately identify senescent-like cells in a manner that is not biased by a highly curated gene list.
While evaluating the SenSig gene lists, we posited that gene lists derived from aged mice would perform better when scoring aged cells. However, this is not the case as the SPIDER+/− young MuSCs and the SPIDER+/− geriatric MuSCs gene lists did not perform well in young (AUC = 0.47) and geriatric (AUC = 0.30) cells, respectively (Fig. 3p).
Refined analysis of myogenic subsets
In the final dataset, we identified two broad myogenic clusters (Fig. 1f). We observed a cluster of MuSCs and progenitors that expressed the myogenic transcription factor Pax7 (ref. 3) and a cluster of myonuclei that expressed Acta1, Myh1 and Myh4, genes critical for the contractile function of mature skeletal muscle cells49 (Extended Data Fig. 3c). We subsetted out these myogenic clusters and re-clustered and re-embedded the cells, resulting in nine distinct subclusters for refined annotation (Fig. 5a,b). We identified four progenitor populations that expressed Pax7, Myf5, Myod1, Myog, Mymk and Mymx and three myonuclei subtypes (IIx, IIb, IIx/IIb) that expressed Acta1, Ckm and Tnnt3 (ref. 49). We observed two transcriptomically variant clusters, which expressed both myogenic markers and either endothelial cell markers like Cd34, Cdh5 and Pecam1 or monocyte and macrophage markers like Ccr2 and C1qa (Extended Data Fig. 7). We suspected these clusters were dominated by doublets. The first doublet subcluster contained cells that coexpressed Pax7 and Pecam1, suggesting that these were MuSCs and progenitors and endothelial cell doublets (Extended Data Fig. 7a). The second doublet subcluster contained cells that coexpressed Acta1 and C1qa, suggesting that these were myonuclei and monocyte/macrophage doublets (Extended Data Fig. 7b). We designated these subclusters as ‘doublets 1’ and ‘doublets 2’, respectively, and excluded them from subsequent analyses involving the myogenic subset.
a–c, Pseudotime organization of the myogenic cell subset after re-clustering. Myogenic cells were re-embedded with UMAP (a) and with PHATE (b) and are colored by general myogenic IDs. PHATE embeddings were used by Monocle 3 to organize the cells in pseudotime where the earliest pseudotime is in the upper-left corner and the latest pseudotime is in the upper-right corner. Myogenic cells were organized into 25 approximately equal bins of increasing pseudotime values (c). d, Each cell was assigned an S phase score and a G2M phase score using Seurat’s Cell Cycle Scoring method. Scatterplot of the two scores with the cells colored by the predicted cell cycle phase based on the two scores. e–g, Polar coordinates (d) were converted to Cartesian coordinates and rescaled to fit in a range from 0 to 1. Cells are colored by the predicted cell cycle phase based on the two scores (e) and by the log-normalized expression of Cdk1 (f) and Cdk4 (g). h–k, Dot plots of the average log-normalized expression of Pax7 (h), Myog (i) and Acta1 (j) and by the fraction of non-G1 (S, G2 and M) cells (k) in each dpi and pseudotime-based myogenic cell-state bin. The size of the circle is the percentage of cells in each pseudotime-based myogenic cell-state bin for each age and dpi combination. l, Top 40 differentially upregulated and downregulated genes between cells in MuSCs 1 and MuSCs 2–6 at day 3.5. All genes highlighted here have an FDR-corrected q value < 0.05.
The type IIx cluster expressed Myh1, the type IIb cluster expressed Myh4, and the type IIx/IIb cluster expressed both Myh1 and Myh4 (Fig. 5b and Extended Data Fig. 8a–g)50. A previous study also identified a myonuclei cluster that expressed both Myh1 and Myh4, but this study concluded that this cluster represented cells with high metabolic activity, not cells that represent a transitional state between types IIx and IIb9. We interpret that the type IIx/IIb cluster identified here represents a transitional state between types IIx and IIb because 39% of the cells in this cluster coexpress Myh1 and Myh4 and because this cluster does not differentially express markers of high metabolic activity (Tnnc2, Tnni2, Mb, Cox6a2, Cox6c, Atp5e, Atp5g1; Extended Data Fig. 8f–h)9,51,52. Additionally, this transitional fiber type is the most common transitional fiber type in rat and mouse muscle fibers53,54. Notably, the type IIx/IIb cluster had a lower percentage of mitochondrial reads than the type IIx and type IIb clusters, indicating that these cells were not clustering together due to being lower quality (Extended Data Fig. 8I). We did not identify a cluster of neuromuscular junction or myotendinous junction cells due to a lack of Chrne or Col22a1 expression (Extended Data Fig. 8h)9,20.
We then re-embedded the myogenic subset using Potential of Heat-diffusion for Affinity-based Trajectory Embedding (PHATE)55, and these embeddings were used by Monocle 3 (refs. 56,57,58) to organize the cells in pseudotime (Fig. 5b,c). The pseudotime values were grouped into 25 bins that contain approximately equal numbers of cells. As the cells progressed through pseudotime, the order in which myogenic markers were expressed followed a typical trajectory of myogenesis. Early pseudotime bins had predominant expression of Pax7 and Myf5, but no strong expression of activation markers. Further in the pseudotime progression, cells still expressed Pax7 and Myf5, but they also expressed cycling markers such as Cdk1/Cdk4. In later pseudotime bins, cells expressed Myod1, Myog, Mymx and Mymk, markers of committed and fusing progenitors. In the latest pseudotime bins, cells expressed Acta1, Ckm, Myh1 and Myh4, markers of myogenic maturation (Extended Data Fig. 9a–h).
Pseudotemporal analysis of myogenesis progression
To directly compare myogenesis in regeneration responses, we assembled an annotated ‘cartography’ of myogenic progression arrayed across both day after injury and myogenic pseudotime. We first examined the percentage of cells that fell within each of the 25 initial pseudotime bins by time point and age, and then used expression frequency of myogenic marker genes to inform myogenic cell-state annotations (Extended Data Fig. 9a–h). Pseudotime bins 1–6 exhibited age-specific cell abundances, with bins 1–2 predominantly containing cells from young mice, bin 3 from young and old mice, bin 4 contained cells from old mice only, bin 5 from old and geriatric mice, and bin 6 from geriatric mice only. Cells in bins 1–6 expressed Pax7 and Myf5, which we annotated as ‘MuSCs’ (with sub-stages 1–6 preserved). Cells in bins 7–13 at 0 dpi expressed Pax7, Myf5 and Myod1 and lowly expressed cycling markers Cdk1 and Cdk4, which we annotated as ‘non-cycling MPCs’ (MPC, myogenic progenitor cell). Cells in bin 7 and dpi 1–3.5 and cells in bins 8–11 and dpi 1–7 expressed Myod1 and cycling markers Cdk1 and Cdk4, which we annotated as ‘cycling MPCs’. Cells in bin 12 at dpi 0–7 and cells in bin 13 at dpi 1–2 have diminishing expression of Myf5, Cdk1 and Cdk4, and increasing expression of Myog and Mymk, which we annotated as ‘committing MPCs’. Cells in pseudotime bin 13 and dpi 3.5–7 highly expressed Myog and Mymk and lowly expressed Cdk1 and Cdk4, which we annotated as ‘fusing myocytes’. Cells in pseudotime bins 14–25 and dpi 0–7 expressed Acta1, which we annotated as ‘myonuclei’. These pseudotime-informed myogenic cell-state aggregates (summarized in Extended Data Fig. 9i) were used in subsequent analyses.
To infer the cell cycle phases, we first assigned each myogenic cell S phase and G2M phase scores using Seurat’s standard Cell Cycle Scoring method59 (Fig. 5d). We treated these scores as polar coordinates, which were converted to Cartesian coordinates and normalized to be within a range from 0 to 1 (‘normalized theta’; Fig. 5e–g and Extended Data Fig. 10a,b). We assessed all myogenic cells within this cell cycle progression from a normalized theta of 0 to 1, corresponding to the continuum of G1–S–G2M stages (Fig. 5e and Extended Data Fig. 10a,b). When considering all MuSCs and progenitors, the distribution of normalized theta values increased from day 0 to day 3.5 and nearly returned to day 0 levels by day 7 in all age groups, suggesting a return to quiescence as expected (Extended Data Fig. 10e). We observed a shift to higher normalized theta values at 1 dpi in the geriatric samples compared to the young and old myogenic cells, suggesting an age-skewed cell cycle induction in early injury response that may represent a precocious activation phenotype. Differences by age group were minimal after 3.5 dpi. When examining all myogenic cells by the 25 pseudotime bins, we observed a shift in normalized theta values at bin 7 persisting through bin 13 (Extended Data Fig. 10f). Notably, pseudotime bins 7–13 also highly expressed the cycling markers Cdk1 and Cdk4 (Extended Data Fig. 9g,h).
We then found that Cdk1 and Cdk4 are more highly expressed in cells predicted to be in S or G2M than in cells predicted to be in G1 (Fig. 5f,g). Seurat’s standard G1 cutoff is at a normalized theta value of 0.25. Based on the expression of cycling markers Cdk1 and Cdk4 and the distribution of cells across the normalized theta values, we extended the G1 cutoff to 0.375 for this analysis (Fig. 5f,g and Extended Data Fig. 10c,d). For simplicity, cells with a normalized theta below and above 0.375 were classified as ‘G1’ and ‘non-G1’ (S/G2/M), respectively (Extended Data Fig. 10c,d).
We then calculated the percentage of cells within each pseudotime-informed myogenic cell state and time point by age group. Within MuSCs 1–6 at 0–1 and 5–7 dpi, we observed high levels of Pax7 (Fig. 5h). We also detected a high fraction of non-G1 cells in MuSCs 1–6, especially at 1–3.5 dpi in the old and geriatric mice (Fig. 5k). This suggested that more of the old and geriatric MuSCs were actively cycling after injury compared to the young MuSCs. We compared the 3.5 dpi cells in MuSCs 1 and MuSCs 2–6 by differential gene expression and found that the quiescence-associated genes Cdkn1c (encoding p57Kip2) and Socs3 were upregulated in MuSCs 1 and numerous translation-associated genes such as Rps29 were upregulated in MuSCs 2–6 (Fig. 5l). These expression profiles suggest that MuSCs 1 cells are in a less activated state MuSCs 2–6 cells.
We observed an inverse relationship between the average Myog expression and the fraction of non-G1 cells in the cycling MPCs, committing MPCs and fusing myocytes, as expected for differentiating myogenic cells (Fig. 5i,k). In all ages, we detected the highest average expression of Myog in the fusing myocytes population (Fig. 5i). We detected myonuclei in all ages and at every dpi, but there were fewer myonuclei with lower average expression of Acta1 at dpi 1–3.5 in the geriatric mice compared to the young and old mice (Fig. 5j). These results present an integrated cellular cartography of myogenic trajectories through regeneration, which exhibits age-associated cellular trajectories, particularly in the MuSC pool.
Senescence scoring across the myogenic cell cycle
To explore how cellular senescence manifests within this organized cartography of myogenesis, we focused on the one-way FBR Sen Score, which performed well across ages (Fig. 3o). To identify cells with senescent-like identities, we defined a threshold within the one-way FBR Sen Score based on its relationship with Cdkn2a and Cdkn1a expression (Extended Data Fig. 10g,h). We set a one-way FBR Sen Score threshold at 2,412, where 50% of cells above this value coexpressed Cdkn2a and Cdkn1a, and classified cells above this threshold as ‘Sen Score-high’ and senescent-like (Fig. 6a). We further observed through ROC analysis that the one-way FBR Sen Score could accurately identify double-positive Cdkn2a and Cdkn1a cells from both the G1 and non-G1 fractions of MuSCs and progenitors (Fig. 6b).
a, Scatterplot of the normalized Cartesian coordinate cell cycle scores and the one-way FBR senescence score (Sen Score) in all MuSCs and progenitors. Cells colored in red coexpress Cdkn2a and Cdkn1a. All other cells are colored in gray. The vertical line is the G1 cutoff, and the horizontal line is where 50% of the cells above this line coexpress Cdkn2a and Cdkn1a. b, ROC curves based on the coexpression of Cdkn2a and Cdkn1a in the MuSCs and progenitors (all ages and dpi) in G1 and non-G1 (S, G2 and M) when using the one-way FBR Sen Score. AUC is reported for each ROC curve. c,d, Dot plots of the fraction of cells that coexpress Cdkn2a and Cdkn1a (c) and the fraction of one-way FBR Sen Score-high cells (d) in each dpi and pseudotime-based myogenic cell-state bin. The size of the circle is the percentage of cells in each pseudotime-based myogenic cell-state bin for each age and dpi combination. e–g, Scatterplots of the fraction of MuSCs 1 (e), MuSCs 2–6 (f) and all MPCs (g) at 3.5 dpi that coexpress Cdkn2a and Cdkn1a by age. Points are the fraction for each sample (n = 2–4) and the horizontal line is the median for each age group. Significance was evaluated using a two-sided Welch’s t-test. h–j, Scatterplots of the fraction of MuSCs 1 (h), MuSCs 2–6 (i) and all MPCs (j) at 3.5 dpi that have a high one-way FBR Sen Score by age. Points are the fraction for each sample (n = 2–4), and the horizontal line is the median for each age group. Significance was evaluated using a two-sided Welch’s t-test.
We aimed to quantify the prevalence of cellular senescence within the cartography of myogenesis by age using both double-positive expression status and Sen Score, reasoning that the scores might capture a more expansive set of senescent-like cells. We observed a correspondence between the fraction of cells that coexpress Cdkn2a and Cdkn1a and the fraction of cells that have a high Sen Score in most ages, cell states and time points (Fig. 6c,d). Focusing on 3.5 dpi, we found a significantly higher fraction of Cdkn2a+ and Cdkn1a+ cells in geriatric MuSCs 1 compared to young MuSCs 1 (Fig. 6e; Welch’s t-test, *P value = 0.007). We did not observe age-specific differences in the fraction of Cdkn2a+ and Cdkn1a+ in MuSCs 2–6 (Fig. 6f). We found a significantly higher fraction of Cdkn2a+ and Cdkn1a+ in both old and geriatric MPCs (from both cycling and committing states) compared to young MPCs (Fig. 6g; Welch’s t-tests, *P values = 0.04 and 0.04, respectively). We observed a similar pattern in the fraction of cells that have a high one-way FBR Sen Score, with significantly higher frequencies in the old compared to young MuSCs 1 (Fig. 6h; Welch’s t-test, *P value = 0.006). We did not observe age-specific differences in the fraction of Sen Score high MuSCs 2–6 or MPCs (Fig. 6i,j). Together, these observations point to a transitory senescent-like cell population that is abundant at the self-renewing MuSC 1 stage across all ages, but increases in older mice, potentially underlying a stalled stem cell self-renewal in mouse muscle aging.
Evaluating muscle regeneration with spatial transcriptomics
To further investigate the aging-associated accumulation of senescent-like myogenic cells, we generated Curio Seeker spatial transcriptomic datasets from young and geriatric TAs at 5 dpi. Deconvolution was performed using the aging skeletal muscle regeneration scRNA-seq dataset presented here, and 15 cell types were identified in the spatial transcriptomic datasets (Fig. 7a,c). We delineated the spatial transcriptomic maps into an injury zone and a noninjury zone (Fig. 7b,d) based on immune cell density (Methods). We predominantly detected myonuclei, especially within the noninjury zone, in both the young and geriatric samples, as expected (Fig. 7e,f). Within the injury zone, we detected a greater diversity of cell types, including MuSCs and progenitors, fusing myocytes, monocytes/macrophages, B cells, T cells, dendritic cells, FAPs and endothelial cells (Fig. 7e,f). We observed some age-specific cell-type composition differences including a higher fraction of fusing myocytes and FAPs in the young injury zone and a higher fraction of monocytes/macrophages, T cells and myonuclei in the geriatric injury zone (Fig. 7e,f). Also, as expected, we detected a higher fraction of MuSCs and progenitors and fusing myocytes in the injury zone than the noninjury zone (Fig. 7e–l). We calculated a one-way FBR Sen Score based on a refined FBR SenSig gene list in the spatial datasets (Fig. 7m,n). There were significantly more FBR score-high senescent-like geriatric MuSCs and progenitors in the whole spatial dataset and specifically in the injury zone (Fig. 7o; Mann–Whitney U test, *P = 0.02 and *P = 0.03, respectively). When considering all other cell types, we found that there were significantly more geriatric senescent-like non-MuSC and progenitor cells in the whole spatial dataset and specifically within the noninjury zone (Fig. 7o; Mann–Whitney U test, *P = 7.53 × 10−38 and *P = 1.35 × 10−36, respectively). These high-resolution spatial transcriptomics maps provide further support for the conclusions from the scRNA-seq results that a senescent-like MuSC and progenitor cell population is more abundant in the geriatric mice, suggesting stalled stem cell self-renewal in aging.
a–d, Spatial transcriptomic maps of young (a and b) and geriatric (c and d) TA muscles 5 days after injury. Beads are colored by cell type (a and c) and injury zone (b and d). Scale bars, 500 µm. e,f, Percentage of cells of select cell types in the injury zone (e) and the noninjury zone (f) in the young and geriatric samples. g–n, Zoomed-in section of the young (g, i, k and m) and geriatric (h, j, l and n) spatial maps (as outlined in a–d) with beads colored by all cell types (g and h), MuSCs and progenitors (i and j), fusing myocytes (k and l) and one-way FBR Sen Score (m and n). o, Box plots of the one-way FBR Sen Score in the MuSCs and progenitors and all other cell types within the full transcriptomic map and within the injury and noninjury zones. The center line is the median, box limits are the upper and lower quartiles, whiskers are the 1.5 times the interquartile range and diamonds are the outliers. From left to right, n = 168, 99, 120, 88, 48, 11, 31,936, 2,2821, 11,361, 9,609, 20,575 and 13,212. Significance was evaluated using a two-sided Mann–Whitney U test.
- SEO Powered Content & PR Distribution. Get Amplified Today.
- PlatoData.Network Vertical Generative Ai. Empower Yourself. Access Here.
- PlatoAiStream. Web3 Intelligence. Knowledge Amplified. Access Here.
- PlatoESG. Carbon, CleanTech, Energy, Environment, Solar, Waste Management. Access Here.
- PlatoHealth. Biotech and Clinical Trials Intelligence. Access Here.
- Source: https://www.nature.com/articles/s43587-024-00756-3