Dissecting Murine Muscle Stem Cell Aging through Regeneration Using Integrative Genomic Analysis.

During aging, there is a progressive loss of volume and function in skeletal muscle that impacts mobility and quality of life. The repair of skeletal muscle is regulated by tissue-resident stem cells called satellite cells (or muscle stem cells [MuSCs]), but in aging, MuSCs decrease in numbers and regenerative capacity. The transcriptional networks and epigenetic changes that confer diminished regenerative function in MuSCs as a result of natural aging are only partially understood. Herein, we use an integrative genomics approach to profile MuSCs from young and aged animals before and after injury. Integration of these datasets reveals aging impacts multiple regulatory changes through significant differences in gene expression, metabolic flux, chromatin accessibility, and patterns of transcription factor (TF) binding activities. Collectively, these datasets facilitate a deeper understanding of the regulation tissue-resident stem cells use during aging and healing.

INTRODUCTION

Physical frailty, with its associated immobility and disability, is a major factor limiting independence and quality of life for the elderly and can be partially derived from skeletal muscle atrophy and weakness (Marcell, 2003). Declines in the health and repair of skeletal muscle can be attributed to a population of resident stem cells called satellite cells (Wang and Rudnicki, 2011) or muscle stem cells (MuSCs). In response to damage, MuSCs undergo dramatic molecular transitions (Porpiglia et al., 2017) to regenerate the tissue as well as replenish the reservoir of stem cells for future regenerative needs. However, during aging, decreasing numbers and function of MuSCs (Blau et al., 2015) result in reductions in the rate and magnitude of recovery following muscle injury, leading to persistent tissue damage (Yin et al., 2013) and potentially contributing to age-associated muscle wasting (Almada and Wagers, 2016).

The molecular mechanisms that govern stem cell aging encompass changes in metabolism (Ren et al., 2017), aberrant chromatin packaging (Cheng, 2015; Singh et al., 2019), accumulation of DNA damage (López-Otín et al., 2013), and loss of proteostasis, all of which collectively converge to drive an imbalance between the maintenance of the quiescent (non-proliferating) state, differentiation, and self-renewal. We lack a mechanistic understanding of how each of these levels of genomic regulation are modified in age for MuSCs or how they influence each other, and this limitation is in part driven by the lack of genome-wide datasets for aged MuSCs, notably during regeneration. Understanding how different types of molecular changes (Liu et al., 2013) from natural aging impact defects and delays in healing through MuSCs is critical for prevention of a senescent (Sousa-Victor et al., 2014) or fibrogenic phenotype (Brack et al., 2007) and therapeutic target discovery (Chakkalakal et al., 2012; Bernet et al., 2014) to maintain healthy muscle into old age (Scharner and Zammit, 2011).

Herein, we use an integrative genomics approach and compare the gene expression programs and chromatin landscape of murine MuSCs from distinct age groups during multiple phases of the regenerative response. Our results describe how natural aging impacts changes in the regulome that manifest in the loss of constitutive heterochromatin and aberrant patterns of transcription factor (TF) binding. We link these alterations to a modulation of metabolic flux and corresponding changes in gene expression. Collectively, these data enable the definition of positive regulatory programs that drive healthy regeneration and how chromatin landscapes evolve to negatively affect regeneration as a result of aging.

RESULTS

Aging Modulates Muscle Regeneration through Differences in Gene Expression

To understand how aging impacts the function of MuSCs, hindlimb muscles (tibialis anterior [TA] and gastrocnemius [Gas]) of wild-type (WT) mice (young: 2–3 months and aged: 22–24 months) were injected with barium chloride (BaCl2; Figure 1A). This injury model yields the destruction of muscle fibers but leaves MuSCs intact, facilitating tissue regeneration (Hardy et al., 2016). Consistent with previous observations (Sousa-Victor et al., 2014), histological and immunofluorescence (IF) staining of injured tissue isolated 7 days post-injury (dpi) from aged and young mice displayed reductions and delays in regeneration for aged muscles (Figure 1B; Figure S1A). To understand the sources of dysregulated repair, MuSCs were purified at multiple time points (0, 1, 3, 5, 7, and 21 dpi) by using fluorescence-activated cell sorting (FACS) (Cerletti et al., 2008) with both negative (Sca-1−, CD45−, Mac-1−, and Ter-119−) and positive surface markers (CXCR4+, β1-integrin+) (Figure 1C) from TA and Gas muscles. Previous studies (Maesner et al., 2016) of young MuSCs have shown that these surface markers enrich a highly purified (>90%) population of Pax7+ MuSCs. In line with these results, 93%–97% of our sorted cells stained positive for Pax7 expression (Figure S1B). We also extracted MuSCs from injured tissue 3 dpi using FACS and immunostained the cells for Pax7 and found >93% were Pax7+ (Figure S1B). To further validate the sorted cells were MuSCs, MuSCs were extracted from uninjured hindlimb muscles of transgenic mice harboring a Cre/LoxP-based system for MuSC lineage tracing Pax7CreERT2/+; Rosa26TdTomato (P7TdT) (Keefe et al., 2015).The transcriptome of FACS-analyzed WT MuSCs with P7TdT MuSCs (Figure S1C) was compared by RNA isolation followed by RNA sequencing (RNA-seq) (Aguilar et al., 2016), and the two transcriptomes yielded highly reproducible correlations (Spearman = 0.97). These results, coupled with the high purity observed from IF imaging, indicate FACS-analyzed WT MuSCs and P7TdT MuSCs exhibit nearly complete overlap.

Figure 1.

Muscle Stem Cells Act Aberrantly as a Result of Aging and Poorly Regenerate Muscle after Injury

(A) Schematic of experiment whereby young (2–3 months) and aged (22–24 months) mice are injured by BaCl2 injection in hindlimb muscles (gastrocnemius and tibialis anterior [TA]) and muscle stem cells (MuSCs) are isolated with FACS, pooled across muscles, and profiled during muscle regeneration by using genome-wide chromatin accessibility and gene expression measurements.

(B) Histological assessment (top) and quantification of mean regenerating fiber cross-sectional area (CSA) size (identified through centrally located nuclei) from TA muscles harvested from young and aged mice 7 days post-injury with hematoxylin and eosin (H&E) staining. *p < 0.05, calculated by two-sided, two-sample Student’s t test assuming equal population variance from n = 4 biological replicates.

(C) Representative FACS plots showing negative (Sca-1, Mac-1, CD45, and Ter-119) and positive (CXCR4 and b1-Integrin) surface markers where numbers within gates indicate percentage of cells within gate.

(D) Heatmap of differential expression for 4,985 genes plotted as Z score for young and aged MuSCs isolated from different days post-injury. Clusters are identified by color on the left side of the heatmap.

(E) Dirichlet Process Gaussian Process (DPGP) mixture model-based clustering of gene expression time series data defined as the Z score of young over aged (gray represents 2× standard deviation and black line is cluster mean), where cluster peaks are corresponding to day of MuSC isolation and color-coded to match clusters in (D).

(F) Enriched GO and KEGG pathways for a subset of clusters from (E) that are color-coded.

(G) Venn diagrams of differentially expressed genes from freshly isolated, uninjured, or activated MuSCs from young and aged mice and human MuSCs24. Activated MuSCs were isolated from 3 dpi from young and aged murine muscle and 7 days in culture for human MuSCs, when both of which have undergone at least one cellular division.

(H) Enriched pathways from unique and shared genes from (G). FI, pathways from common genes among all freshly isolated MuSCs; FI–M, pathways from common genes among mouse freshly isolated MuSCs; FI–A, pathways from genes among aged mouse freshly isolated MuSCs; FI–Y, pathways from genes among young mouse freshly isolated MuSCs; Ac, pathways from common genes among all activated MuSCs; Ac–M, pathways from common genes among mouse activated MuSCs; Ac–A, pathways from common genes among aged mouse activated MuSCs; Ac–Y, pathways from common genes among young mouse activated MuSCs.

Gene expression profiles from young and aged sorted MuSCs were generated, and strong agreement (Spearman > 0.80) was observed for each biological replicate isolated from the different time points. Differential expression (DE) analysis revealed that 4,985 genes underwent a change during the time course or as a result of aging. Time series clustering of DE trajectories using non-parametric-based clustering (McDowell et al., 2017) revealed differences between young and aged MuSCs and also displayed variations in enriched Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (Figures 1D–1F). At 1 dpi, enriched pathways for young MuSCs included cellular adhesion and extracellular matrix (ECM) deposition (Ret, Itgb1, Itga6, Col5a3, Laminin β1, and Col15a1), whereas aged MuSCs overexpressed genes associated with histone deacetylase binding (Rac1, Cbx5, and Nudt21). At 3 dpi, both young and aged samples were positive for many cell cycle genes (Ccna2, Ccnb1, Ccnb2, and Ccne1) and pathways, but young samples were observed to increase the expression of TFs, such as Myf5, Myog, and Runx1 (Umansky et al., 2015), and genes associated with muscle contraction (Myh3, Tnnt2, and Myomaker-Tmem8c; Figure S1D). In contrast, at days 3 and 5, aged MuSCs exhibited upregulation of inflammatory markers (Tnfrsf10b/11a/13bf12a, and Irf7) and nervous system development genes (Ntrk3, Cadm1, and Nrn1). At 5 dpi, young MuSCs upregulated pathways associated with chromatin binding (Top2a, Ezh2, Cbx5, MyoD1, and Setdb1) and maintained increases in the expression of genes associated with myogenic differentiation (Myh3, Tnnt2, and Tmem8c), whereas aged MuSCs maintained increases in the expression of inflammatory response genes (Cxcl12, Il-10, Fas, and Ccl3) and complement activation (C1ra, C1s1, and C3). At 7 dpi, young MuSCs upregulated genes associated with mitogen-activated protein kinase (MAPK) and calcium signaling as well as inhibition of Wnt signaling (Camk2g, Sphk1, Mef2c, Sfrp1, and Sfrp5). In contrast, aged MuSCs continued to be enriched for immune response pathways. At 21 dpi, young MuSCs had more differentially expressed genes than aged MuSCs (765 and 161, respectively) and included genes such as Notch2, Zbtb20, and Laminin β1. Many of the enriched genes in aged MuSCs at 21 dpi were also shared with aged MuSCs from 0 dpi and young MuSCs from 1 dpi. Collectively, the observed changes in expression align with histological observations as well as results of other studies (Liu et al., 2013) showing aging induces impairments in the ability to maintain quiescence and regenerate tissue efficiently, which was in contrast to the responses for young MuSCs that displayed stronger and faster regeneration.

Murine MuSCs Display Similar Activation Programs as Human MuSCs

To glean if young or aged murine MuSCs displayed similar changes in expression as a result of activation with human MuSCs, we contrasted the transcriptomes of freshly isolated (FI), uninjured and activated human MuSCs (Charville et al., 2015) with murine uninjured (0 dpi) and activated (3 dpi) MuSCs. By 3 dpi and 7 days in culture, both murine and human MuSCs have fully activated and undergone at least one cellular division (Rodgers et al., 2014). A total of 135 common genes were observed to exhibit DE between FI murine young and aged and FI human MuSCs (Figure 1G). In all 3 FI populations, we observed enrichment of genes such as Calcr, Apoe, Cebpb, Ndrg2, Chrdl2, and Chodl, which have been previously associated with quiescence (Fukada et al., 2007) and decreased in expression after activation. Common enriched pathways included regulation of proliferation (false discovery rate [FDR] = 1.6e-06) and ECM components (FDR = 3.3e-06). Several ECM proteins unique to FI murine MuSCs included vitronectin (Vtn) and decorin (Dcn) and were more strongly expressed in young than aged MuSCs. For all 3 activated MuSC populations, we observed increases in the expression of a series of myogenic differentiation and muscle contractile genes (MyoG, Tmem8c, Myh3, Cav3, Tnnt2, Acta1, and Des). Pathway annotation of activated murine MuSCs revealed enrichments in cell cycle (cell cycle phase, FDR = 2.55e-35; M phase, FDR = 3.38e-35) and cytoskeleton assembly (cytoskeletal, FDR = 6.48e-15; contractile fiber, FDR = 2.66e-07). Pathways enriched in human activated MuSCs (versus FI human) and young activated murine MuSCs included negative regulation of apoptosis (FDR = 2.12e-1), which were absent in activated aged murine MuSCs. In summary, these results suggest similarities in several genes associated with quiescence and activation for human and mouse MuSCs but distinct genes and pathways for ECM-related proteins.

Aged MuSCs Undergo Switches in Metabolism That Associate with Histone Methylation

Many differentially enriched pathways observed from murine MuSC gene expression changes were related to metabolism, which is consistent with the known regulatory role of metabolism on MuSC gene expression (Berger, 2016). To probe deeper into these effects, we used genome-scale metabolic modeling to assess the relationship between 3,744 metabolic reactions, 2,766 metabolites, 1,496 metabolic genes, and 2,004 metabolic enzymes (Duarte et al., 2007; Figure 2A). The expression data were overlaid onto the metabolic model by maximizing the flux through the metabolic reactions that are associated with upregulated genes, while minimizing flux through those reactions that are downregulated genes. The model predicted that young MuSCs upregulated metabolic flux through vitamin A (VA) (retinol) signaling, glutamate metabolism, and oxidative phosphorylation (OxPhos). Previous measurements (Pala et al., 2018) of the bioenergetics (oxygen consumption rate and extracellular acidification rate) of young and aged MuSCs showed aging induces a shift in metabolic substrate use away from oxidative metabolism, which is consistent with the results from the metabolic model. The model also revealed increased flux through folate metabolism and the one-carbon cycle (1CC) for aged MuSCs.

Figure 2.

Alterations in Metabolism Associate with Global Changes in Histone Methylation of Young and Aged MuSCs

(A) Histogram (top) and jitter plot (bottom) showing significant changes in reaction flux plotted as the Z score for each metabolic reaction between young and aged MuSCs using a paired t test (p < 0.05).

(B) Representative immunofluorescence (IF) staining of total histone levels (H3) and repressive chromatin modifications (H3K27me3 and H3K9me3) for young and aged MuSCs. Scale bar represents 100 µm.

(C) Quantification of stains from (B) show higher histone levels (H3) and constitutive heterochromatin modifications (H3K9me3) for young MuSCs, where aged MuSCs displayed an increase in facultative heterochromatin modifications (H3K27me3). ****p < 0.001, as calculated by two-sided, two-sample Student’s t test. n = 1,765–3,825 cells from each of two young and two aged mice.

(D) Heatmap of gene expression for methyltransferases and chromatin enzymes in young and aged uninjured MuSCs plotted as Z score.

Based on the role of 1CC to generate s-adenosyl-methionine (SAM) substrates, which act as methyl group donors for methylation reactions, we next examined whether changes in 1CC coincided with changes to histone methylation (Mentch et al., 2015). To address this question, single uninjured MuSCs from young and aged mice were immunostained, imaged, and enumerated for global levels of total histone levels (H3) and repressive chromatin modifications (H3K27me3 and H3K9me3) (Figures 2B and 2C). Aged MuSCs displayed increased levels of H3K27me3, which was in contrast to young MuSCs that displayed increases in H3 and H3K9me3. These results are consistent with young MuSC increases in the expression of Suv39h methyltransferase enzymes, which are primarily responsible for the formation and maintenance of constitutive heterochromatin (Peters et al., 2003), as well as other chromatin enzymes that contain the Su(var)3–9, Enhancer-of-zeste, and Trithorax (SET) domain (Ezh1, Utx/Kdm6a, Ash1l, Suv39h2, and Setd3) and interact with Suv39h enzymes. These enzymes were decreased in aged MuSCs, which also overexpressed different factors such as Mll5/Kmt2e, Hp-1γ/ Cbx3, Prdm16, Eed, and Setd6 (Figure 2D). Taken together, these results show that aged MuSCs possess differences in metabolic flux that are associated with alterations of histone modifications and changes in the expression of concomitant sets of enzymes. These findings are consistent with previous studies of cellular aging (Liu et al., 2013; Sen et al., 2016), whereby histone loss (Liu et al., 2013) and constitutive heterochromatin are modified.

Engagement of Retinoic Acid Receptors Restrains MuSC Activation but Is Lost with Age

VA-retinoic acid has been shown to restrain human skeletal muscle progenitors from differentiation (Ryan et al., 2012; El Haddad et al., 2017) and to promote quiescence in hematopoietic stem cells (Cabezas-Wallscheid et al., 2017). A loss of the expression of retinoic acid receptors and retinoid x receptors was observed for aged MuSCs (Figure 3A), mirroring the loss of other quiescence-associated genes (Figure S1D). To probe if VA-retinoic acid, as predicted by our metabolic model, impacted young or aged murine MuSC maintenance of quiescence, we isolated MuSCs from uninjured hindlimb muscles of young and aged mice using FACS and cultured the cells in activating conditions for 3 days with and without all-trans retinoic acid (ATRA). We observed ATRA-treated MuSCs reduced activation (MyoD+) and proliferation (Ki67+) in both young and aged MuSCs, and ATRA-treated young MuSCs increased Pax7, which was in contrast to aged MuSCs that maintained the same level of Pax7 (Figures 3B and 3C). These results confirm that aged MuSCs lose VA signaling and suggest the loss of retinoic acid in age contributes to the inability to restrain MuSCs from activation.

Figure 3.

Retinoic Acid Receptors Contribute to Maintenance of MuSC Quiescence but Are Lost in Age

(A) Heatmap of gene expression for genes encoding retinoic acid receptors and retinoid x receptors as well as several downstream retinoic acid target genes in young and aged uninjured MuSCs, plotted as Z scores of transcrips per million (TPM) values.

(B) Representative IF staining of Pax7 and MyoD for young and aged MuSCs following 3 days of treatment of retinoic acid (+RA) or DMSO alone (−RA). Scale bars represents 50 µm in images of aged MuSCs and 25 µm in images of young MuSCs.

(C) Quantitation of images in (B) using a two-sample Student’s t test shows that treatment with ATRA increased Pax7 in young MuSCs (**p < 0.01), decreased MyoD in both young MuSCs and aged MuSCs (****p < 0.0001 and *p < 0.05, respectively), and decreased Ki67 in both young MuSCs and aged MuSCs (*p < 0.05 and **p < 0.01, respectively). Cells were harvested from muscles of two young and two aged mice. In total, at least 50 cells were stained and analyzed per condition. A two-sided, two-sample t test was used to calculate statistical significance.

MuSC Regulatory Networks Change with Age

To obtain a deeper insight into the regulation of the aberrant expression networks observed in aging and changes in metabolism, an assay for transposase-accessible chromatin followed by sequencing (Buenrostro et al., 2013) (assay for transposase-accessible chromatin using sequencing [ATAC-seq]) was used (Figure 1A) on MuSCs at multiple time points (Figure 4A). A total of 238,563 accessible sites were found, and blacklisted sites were removed (Amemiya et al., 2019). Overall, the datasets were highly reproducible across technical and biological replicates (Figure 4A). Principal-component analysis (PCA) of enriched sites showed young samples migrate along a trajectory defined primarily by the first principal component, whereas aged MuSCs used a trajectory defined by the second principal component (Figure 4B). During the early stages of the regenerative process (1 and 3 dpi), both young and aged samples were observed to cluster with each other, and at later time points (5 and 7 dpi), the samples began to reveal a divergence in their regenerative trajectory back toward the uninjured state (0 dpi). Annotation of the ATAC-seq peaks revealed the largest fraction of sites aligned with transcriptional start sites (TSSs), but the majority of peaks were distal to TSSs (Figure S3C). During the regenerative process (1, 3, and 5 dpi), young MuSCs were observed to increase the number of distal sites compared with aged MuSCs, but these dynamics reversed at 7 dpi and resembled 0 dpi. To probe if sites that were opening up as a result of aging were previously marked by facultative heterochromatin modifications (H3K27me3) in young uninjured MuSCs, enriched sites were intersected with chromatin immunoprecipitation sequencing (ChIP-seq) datasets (Liu et al., 2013) for H3K27me3 in MuSCs (Figure 4C). Of the 23,275 sites with increased accessibility in the aged samples, 274 overlapped H3K27me3-bound regions in aged satellite cells, while 399 overlapped H3K27me3-bound regions in young MuSCs (Figure 4C). Integrating these results shows that aging induces increases in the accessibility of sites that were previously demarcated by facultative heterochromatin.

Figure 4.

Chromatin Accessibility Is Modified during MuSC Regeneration and Exhibits Divergent Regenerative Trajectories in Aging

(A) Heatmap of Spearman correlation coefficients for individual replicates isolated from age and time points showing strong reproducibility. Correlation was computed on asinh (counts per million) reads after removal of the contributions from surrogate variables.

(B) Multi-dimensional scaling (MDS) of ATAC-seq enrichments color-coded by day of isolation; circles represent aged samples, and triangles represent young samples.

(C) Intersection of differentially accessible (young versus aged) ATAC-seq peaks from day 0 with H3K27me3 sites previously derived for young and aged MuSCs.

(D) Statistically enriched (Benjamini-Hochberg corrected, p < 0.01) pathways from different days derived from ATAC-seq enrichments using GREAT analysis.

An assessment of potential drivers of changes in chromatin dynamics was performed by grouping chromatin regions with similar accessibility profiles through time, and differential peaks were analyzed with gene set enrichment analysis using the GREAT toolbox (McLean et al., 2010; Figure 4D). This analysis yielded pathway enrichments, such as fatty-acyl-coenzyme A (CoA) biosynthesis for young MuSCs before injury, interleukin-6 (IL-6)-mediated signaling, and muscle cell development enriched for aged MuSCs. After injury, aged MuSCs remained enriched for IL-6 signaling and upregulated Wnt receptor signaling, the Ras pathway, and regulation of retinoblastoma protein. In contrast, after injury, young MuSCs exhibited temporal activation of the Notch-mediated Hes/Hey network, fibroblast growth factor (FGF) signaling, as well as muscle cell differentiation.

Aging Induces Divergent TF Binding Dynamics

To further determine candidate regulators that specify differences in the chromatin accessibility landscape, enriched regions were assessed for TF binding motifs and corresponding changes in expression at each time point. This analysis revealed young MuSCs contained increases in the number of enriched motifs and higher expression for myogenic factors (Mef2a, Mef2c, and MyoG) (Figure 5A; Figures S4A and S4C). Many of the recovered motifs were central regulators of the myogenic lineage (Myf5 and Myog) or were predicted to interact with MyoD, such as Zbtb18/Rp58, Tcf3, and Xbp1 motifs (Figure S4). These collective observations are well illustrated by chromatin profiles of the MyoG and Mef2a loci (Figure 5B). The MyoG promoter displayed differential accessibility in aged MuSCs compared with young throughout the time course, and these chromatin changes mirrored the expression profile that peaked at 3 and 5 dpi for young MuSCs. The MyoG promoter contains MyoD binding sites as well as two overlapping binding sites for Pbx/Meis that cooperate (Sartorelli and Puri, 2018) with MyoD. Similarly, distal sites for Mef2a exhibited accessibility changes comparing young and aged MuSCs, which also contained MyoD binding sites, and correlated with stronger expression levels in young MuSCs during regeneration.

Figure 5.

Aging Engenders Variations in Transcription Factor Binding Dynamics during Regeneration

(A) Comparison of fold change in transcription factor expression in young and aged (x axis) with the HOMER fold enrichment of the corresponding transcription factor binding motif in young and aged over background (y axis). Motifs enriched at day 0 in either young (y axis > 0) or aged (y axis < 0). Points on the chart are sized proportionally to median TPM for the corresponding day and color-coded by the motif family to which they belong.

(B) Normalized tracks of ATAC-seq datasets (fold change of ATAC-seq signal relative to the mm10 background distribution) around the MyoG locus (left) and Mef2a (right), where differences in enrichments are highlighted in gray and all tracks are scaled to the same level (fold change, 0–20).

(C) Top: comparison of aggregate MyoD footprints in young versus aged for each day of isolation, where the y axis values are are logarithms of reads per motif site and the x axis is the distance away from the center of the footprint (+/− 0.1 kb from motif center).

(D) Schematic of approach to distinguish direct and indirect binding of MyoD transcription factor whereby MyoD ChIP-seq peaks are categorized for ATAC-seq footprints (FP) and position weight matrix (PWM).

(E) Annotated distribution of MyoD ChIP-seq peaks (y axis) in young (Y) and aged (A) MuSCs for each time point of isolation before and after injury.

(F) Gene expression heatmap plotted as Z score for each time point of isolation before and after injury for MyoD co-binding partners identified through indirectly bound MyoD sites and PWM of a second factor.

Given the pivotal role of MyoD in specifying the chromatin landscape of MuSCs, we performed footprinting analysis of MyoD in both young and aged ATAC-seq datasets and found clear footprints across the genome (Figure 5C). On average, we observed stronger and more MyoD footprints for aged MuSCs at 0 and 7 dpi, and at 1, 3, and 5 dpi, MyoD footprints were stronger and increased in number in young MuSCs. These results were consistent with increases in the number of accessible sites at 0 and 7 dpi for aged MuSCs. Given that MyoD functionally cooperates with other TFs and histone acetyltransferases (protein-protein interactions) to activate transcription at target genes, we integrated ChIP-seq datasets for MyoD (Cao et al., 2010) with MyoD footprints and identified ChIP-seq peaks containing a MyoD footprint motif (directly bound MyoD), ChIP-seq peaks lacking the MyoD motif or footprint (indeterminate sites), and ChIP-seq peaks overlaying a motif but lacking a footprint (indirectly bound MyoD) (Figures 5D and 5E). The fraction of ChIP-seq peaks predicted to encompass direct binding versus being indirectly bound changed throughout the time course, with 3 dpi representing the largest number of directly bound predicted sites. Increases in direct MyoD binding were observed for young MuSCs for all time points compared to aged MuSCs, with the greatest differential at 5 dpi, which is the time that exhibited the greatest number of upregulated genes for young MuSCs (Figure S3B). In contrast, aged MuSCs contained more indirect or indeterminate MyoD binding sites for all time points, suggesting more but weaker MyoD interactions with chromatin as a result of age. We next determined the frequency of indirectly bound MyoD sites and queried if the sites coincided with motifs of a second factor. This analysis recovered enrichments for E-protein motifs (Tcf3/E47, Tcf12, and E2f1); AP-1; and related basic helix loop helix (BHLH) proteins, such as MyoG, Myf5, Prdm16, Ddit3, and Tfap4, a transcriptional repressor (Jung et al., 2008). Many of the associated factors exhibited a stronger expression in young MuSCs compared with aged MuSCs during regeneration (Figure 5F), indicating young MuSCs may exhibit a stronger regenerative response due to combinations of TF co-occupancies. We also observed that Prdm16 and Ddit3 increased expression at 21 dpi in aged MuSCs, which was similar to 0 dpi and in contrast to young MuSCs.

Misregulation of DDIT3 from Aging Inhibits Myogenic Differentiation

MyoD footprinting analysis revealed that different TFs coincided with MyoD binding sites, and these TFs and chromatin enzymes displayed DE at different time points. One of the factors that was overexpressed in aged MuSCs was Ddit3/Chop (DNA damage inducible transcript 3) (Oyadomari and Mori, 2004). Ddit3 is a TF that regulates the growth hormone receptor (Ghr) (Zecchini et al., 2019) and protein synthesis through insulin growth factor 1 (Igf1) (Han et al., 2013) and has been shown to act as a transcriptional repressor of MyoD (Alter and Bengal, 2011) through interaction with HDAC1. Ddit3 is also sensitive to autophagy (García-Prat et al., 2016), the protein degradation process through which components of the cytoplasm are digested by lysosomes. The role of Ddit3 in MuSC aging has not been profiled, and to evaluate the impact of reductions of Ddit3, Ddit3 was silenced in myogenic progenitors (C2C12 cells) by delivery of two distinct Dicer-substrate small interfering RNAs (DsiRNAs) (Kim et al., 2005) packaged in lipid nanoparticles, and knockdown efficacy was verified using qPCR (Figure S5). After knockdown, myoblasts were differentiated and fused into myotubes, and increases in fusion and MyoG expression were observed for Ddit3 knockdowns compared to controls (Figures 6A–6C). To determine if Ddit3 knockdown could rescue the differentiation delay and deficiency of aged MuSCs ex vivo, MuSCs were isolated from young and aged mice and derived into myoblasts. Ddit3 was silenced, and myoblast differentiation and fusion in myotubes were compared with those of controls. We observed Ddit3 knockdown enhanced the myotube fusion index (Figure 6B) and increased MyoG expression. These results confirm that coordinative TF binding modified in age can restore transcriptional output that regulates cellular differentiation.

Figure 6.

Silencing Ddit3 Restores Myogenic Differentiation in Aged MuSCs

(A) Representative IF staining of aged MuSCs differentiated into myoblasts or myotubes for 3 days after knockdown (KD) of Ddit3, where myosin heavy chain 3 (green) and DAPI (blue) are shown. Scale bar represents 100 µm)

(B) Enumeration of fusion index for differentiated myoblasts shows KD of Ddit3 induced more fusion. n = 3 independent experiments. *p < 0.05, calculated by two-sample Student’s t test assuming equal population variance.

(C) qPCR of Myogenin (MyoG) following siRNA KD of Ddit3 and 3 days post-differentiation shows upregulation of MyoG. Bars show mean ± standard deviation. n = 3 replicates; ns denotes not significant (p > 0.05), calculated by Student’s t test assuming equal population variance.

DISCUSSION

A healthy MuSC compartment is critical for maintaining muscle homeostasis through aging, and transcriptional and epigenetic regulation (Ermolaeva et al., 2018) protect these cells from deleterious factors present in their systemic milieu that promote alternative cell fates (Brack et al., 2007) and adoption of a senescent (Sousa-Victor et al., 2014) phenotype. Herein, we present an integrative genomics resource showing how the transcriptional and epigenetic landscape of MuSCs is altered during regeneration and aging. We observe that aged MuSCs display significant changes in their gene expression profile and chromatin landscape prior to injury, and we attribute these differences to the use of distinctive metabolic pathways.

Profiling the expression of aged and young MuSCs has shown aging induces a defect in the regenerative abilities of MuSCs (Liu et al., 2013; Bernet et al., 2014; Boonsanay et al., 2016) and an impaired capacity to self-renew (Bernet et al., 2014; Blau et al., 2015). Consistent with these results, we observed decreases in the expression of genes associated with quiescence and negative regulation of cellular proliferation before injury in aged MuSCs. During regeneration, young MuSCs displayed a stronger expression of myogenic TFs and genes associated with muscle contraction than aged MuSCs. These results suggested aging induces restraint of the strength and timing of the myogenic program through impairments in the ability of myogenic TFs to enact changes in expression. To probe further into this effect and the potential role of changes in chromatin, we first observed variations in the expression of chromatin enzymes, with young uninjured MuSCs displaying increases in the expression of oxygen-sensitive factors (Chakraborty et al., 2019) that remove H3K27me3 (Utx) (Faralli et al., 2016; Chakraborty et al., 2019) and interact with MyoD1 and MyoG (Setd3) (Eom et al., 2011). We also observed young MuSCs had increased expression of Kmt5b/Suv4–20h1, which has been shown to regulate heterochromatin (Cheng, 2015). Aged MuSCs were observed to increase the expression of Prdm16, which contributes to cell fate (An et al., 2017) decisions toward brown adipocytes (Seale et al., 2008), and downregulated Lsd1/Kdm1a as well as E2F4 (An et al., 2017), both of which have been shown to prevent brown adipocyte differentiation in MuSCs (Tosic et al., 2018). These results suggest aging induces a loss of chromatin enzymes that deposit constitutive heterochromatin (Peters et al., 2003) and attenuate the ability of myogenic TFs to enact tissue repair.

The regulation of MuSC metabolism (Ryall et al., 2015) has been shown to be disrupted as a result of aging (Zhang et al., 2016; Pala et al., 2018). Leveraging our expression datasets into an unbiased metabolic model, we observed retinoic acid-induced signaling was enriched in young MuSCs but lost in aging. This pathway, involving metabolism of VA, has been shown to restrain human skeletal muscle progenitors from differentiation (El Haddad et al., 2017) and dormancy of hematopoietic stem cells (Cabezas-Wallscheid et al., 2017). Culturing young and aged MuSCs in retinoic acid prevented activation and reduced proliferation, suggesting the loss of VA and associated engagement with retinoic acid receptors with age contribute to a loss of the ability to maintain quiescence and resist activation. The metabolic model also predicted enrichments in OxPhos for young MuSCs and 1CC for aged MuSCs. The coupled loss of OxPhos and gain of SAMs from the 1CC observed in aging has previously been shown to reduce NAD-dependent histone deacetylases (HDACs), such as sirtuins (Verdin, 2015), and increase chromatin accessibility (Ducker and Rabinowitz, 2017). These studies are in line with our observations for aged MuSCs, where Sirt1 was downregulated and decreased in histone levels, and H3K9me3, a marker for constitutive heterochromatin, were detected. These results were also mirrored by decreases in the components of the nuclear envelope (Lamin B1, LaminB2, Lamin B receptor, Sun1, Nesprin1, and Nesprin2), suggesting an untethering of previously condensed heterochromatin. The redistribution of constitutive heterochromatin (H3K9me3) into facultative heterochromatin (H3K27me3) was consistent with increases in the expression of Eed and Ezh2, which form a complex with Polycomb repressive complex 2 (Prc2) to deposit methyl groups onto lysine 27 of histone 3 and display specificity (Margueron et al., 2009) for H3K27me3. Notably, expression of Ezh1, the paralog to Ezh2 that has weaker methyltransferase activity and also complexes with Prc2, was reduced with age. Whether this switch in comparative abundance impacts levels of mono- and di-methylated H3K27 in favor of tri-methylation could be a fruitful direction for further exploration. Overall, the coupled loss of heterochromatin (Pinheiro et al., 2012) and nuclear envelope proteins with age (Wang et al., 2016) implies the hierarchical organization of the MuSC nucleus may be disrupted (Haithcock et al., 2005), which has been observed in aged hematopoietic stem cells (Chambers et al., 2007; Grigoryan et al., 2018) and has also been shown to disrupt gene repositioning schemes and expression during myogenesis (Robson et al., 2016).

The regulatory landscape of young and aged MuSCs showed aging induces changes in the number of accessible chromatin states before and after regeneration. A detrimental result of these differences in the chromatin state were alterations in the expression and manner of binding of TFs (such as MyoD). Footprinting the MyoD motif during the regenerative process (1, 3, and 5 dpi), we observed increases in accessibility for young MyoD motifs compared to that in aged. To further discern the nature of differences in MyoD binding and potential co-binding partners that would impact the epigenetic landscape, we integrated MyoD ChIP-seq datasets and categorized the ChIP peaks and ATAC-seq peaks containing MyoD footprints. This analysis revealed MyoD exhibited enhancements in MyoD direct binding for young MuSCs and stronger interactions with additional binding partners, such as E-proteins. These results were in contrast to aged MuSCs that displayed increases in indirect or indeterminate MyoD binding. Given MyoD’s limited ability to enact transcriptional activation and that MyoD binding correlates with only a small number of DNA elements that change expression, these results highlight how heterodimerization of MyoD with other factors to activate muscle-specific genes may occur from a distance and these distance-dependent changes appear attenuated with age. Another factor that may contribute to weaker interaction of MyoD with DNA and transcriptional activation in aging is a lack of histone acetylation (Puri et al., 1997; Sen et al., 2016). Hypoacetylation at myogenic differentiation genes would be compatible with our results given acetyl-CoA and glucose metabolism regulate the histone acetylation landscape, and both were modified in aged MuSCs. This model is further supported by the fact that the recruitment and complexation of MyoD with histone acetyltransferases (HATs), such as p300/CBP or p/CAF at enhancers (Pliner et al., 2018), are essential for the proper execution of myogenic differentiation, and factors such as AP-1 (Fos/Jun) and Runx1 that cooperatively bind with MyoD were differentially expressed in old age. Integrating these results, we suggest that changes in the cooperative nature of TFs (Chronis et al., 2017) with aging occur at enhancers and that these changes drive reductions in the regenerative potential of aged MuSCs.

A deleterious effect of prolonged MuSC stress from aging is a loss of proteostasis, where impairments in unfolding of proteins cause stress in the endoplasmic reticulum (ER). The accumulation of unfolded or improperly folded proteins induces autophagic clearance, reactive oxygen generation (ROS), and persistent activation of ER stress-mediated cell death signaling, which is regulated by the Protein Kinase R (PKR)-like ER kinase (PERK)-eukaryotic initiation factor 2α (eIF2α)-activating transcriptional factor 4 (Atf4) pathway. Autophagy and the PERK-eIF2α-Atf4 pathway are essential for MuSC survival and function (Xiong et al., 2017), but less information is known about the other signaling arms of the unfolded protein response (UPR) on MuSC function and in aging. Ddit3, which is an autophagic sensitive downstream target of Atf6 (one of the three UPR branches) was found to be upregulated in aged MuSCs at 0 and 21 dpi. Ddit3 is a transcriptional repressor of MyoD (Sartorelli and Puri, 2018), and we observed that indirectly bound MyoD sites were enriched for Ddit3, suggesting reduction of this factor would exert a positive effect on subsequent myogenic differentiation lost in aging. Consistent with this view, transient Ddit3 knockdown in aged MuSCs enhanced myogenic differentiation by increases in MyoG expression. Given Ddit3 also interacts with HDAC1, which associates with MyoD (Mal et al., 2001) in a manner that results in transcriptional repression by deacetylation, these experiments further suggest that imbalances in the combinatorial network of TFs from aging results in inefficient and/or delayed activation of the myogenic program in MuSCs.

Overall, our studies are a valuable resource for understanding how epigenetic regulation affects aging of a tissue-resident stem cell and may be applicable for other types of trauma, such as volumetric muscle loss (Aguilar et al., 2018), and can also be used to contrast anti-aging models. Future work will investigate the complex and combinatorial control of TF binding in aged MuSCs, along with how constitutive heterochromatin domains diminish with age.

STAR★METHODS

RESOURCE AVAILABILITY

Lead Contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Carlos Aguilar (caguilar@umich.edu).

Materials Availability

This study did not generate new unique reagents.

Data and Code Availability

The datasets generated during this study are available at the Gene Expression Omnibus (GEO) repository under GSE121589. Additionally, the bioinformatics code is available on GitHub (https://github.com/annashcherbina/nobel_lab_projects/tree/master/age_V2). Additional datasets used in this analysis are available from the ENCODE portal (encodeproject.org) and GEO portal (https://www.ncbi.nlm.nih.gov/gds) at the identifiers indicated in the ‘‘Other’’ section of the Key Resource Table.

KEY RESOURCES TABLE

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
APC anti-mouse Ly-6A/E (Sca-1), clone: D7, isotype: Rat IgG2a, κ	BioLegend	Cat# 108112; RRID:AB_313349
APC anti-mouse CD45, clone: 30-F11, isotype: Rat IgG2b, κ	BioLegend	Cat# 103112; RRID:AB_312977
APC anti-mouse/human CD11b, clone: M1/70. Isotype: Rat IgG2b, κ	BioLegend	Cat# 101212; RRID:AB_312795
APC anti-mouse TER-119, clone: TER-119, isotype: Rat IgG2b, κ	BioLegend	Cat# 116212; RRID:AB_313713
PE anti-mouse/rat CD29, clone: HMβ1–1, isotype: Armenian Hamster IgG	BioLegend	Cat# 102208; RRID:AB_312885
Biotin Rat Anti-Mouse CD184, clone: 2B11/CXCR4 (RUO), isotype: Rat IgG2b, κ, lot # 6336587	BD Bioscience	Cat# 551968; RRID:AB_394307
Streptavidin PE-Cyanine7, lot # 4290713	eBioscience	Cat# 25–4317-82; RRID:AB_10116480
Mouse anti-human myosin (embryonic), Isotype: Mouse IgG1, supernatant	Developmental Studies Hybridoma Bank	Cat# F1.652 s
Mouse anti-Pax7, Isotype MIgG1, kappa light chain, supernatant	Developmental Studies Hybridoma Bank	Cat# Pax7 s
Anti-Histone H3, rabbit polyclonal	Abcam	Cat# ab1791; RRID:AB_302613
Tri-Methyl-Histone H3 (Lys27), rabbit monoclonal, clone C36B11	Cell Signaling Technology	Cat# 9733; RRID:AB_2616029
Tri-Methyl-Histone H3 (Lys9), rabbit monoclonal, clone D4W1U	Cell Signaling Technology	Cat# 13969 RRID:AB_2798355
Donkey anti-rabbit IgG (H+L), Alexa Fluor 647	Thermo Fisher	Cat# A-31573; RRID:AB_2536183
Rabbit anti-mouse laminin 1+2, Isotype: Polyclonal IgG	Abcam	Cat# ab7463,;RRID:AB_305933
Goat anti-mouse IgG1, Alexa Fluor 488 conjugate	Invitrogen	Cat# A-21121; RRID:AB_2535764
Goat anti-rabbit H+L, Alexa Fluor 555 conjugate	Invitrogen	Cat# A-21428; RRID:AB_141784
Goat anti-mouse H+L, Alexa Fluor 555 conjugate	Invitrogen	Cat# A28180; RRID:AB_2536164
Anti-Pax7 antibody, mouse monoclonal IgG1, Alexa Fluor 488 conjugate	Santa Cruz Biotechnology	Cat# sc-81648; RRID:AB_2159836
Anti-MyoD antibody (G2), mouse monoclonal IgG2b, Alexa Fluor 647 conjugate	Santa Cruz Biotechnology	Cat# sc-377460; RRID:AB_2813894
Anti-MyoD antibody (G2), mouse monoclonal IgG2b, Alexa Fluor 647 conjugate	Santa Cruz Biotechnology	Cat# sc-377460; RRID:AB_2813894
Anti-Ki-67antibody, mouse monoclonal IgG1, PE conjugate	Santa Cruz Biotechnology	Cat# sc-23900; RRID:AB_627859
Chemicals, Peptides, and Recombinant Proteins
Dispase II (activity ≥ 0.5 units/mg solid)	Sigma	D4693–1G
Collagenase Type II (654 U/mg, non-specific proteolytic activity 487 U/mg)	Life Technologies	17101015
Barium Chloride (BaCl2) solution (20% w/v)	Fisher Scientific	SB822–1
Digitonin	Promega	G9441
NP40	Sigma/Roche	11332465001
1M Tris-HCl pH 7.5	Invitrogen	15567–027
1M MgCl2	Thermo Fisher	AM9530G
DMEM, high glucose, pyruvate	Life Technologies	11995065
Ham’s F-10 Nutrient Mix	Life Technologies	11550043
Fetal Bovine Serum	Life Technologies	10437028
HBSS, no calcium, no magnesium, no phenol red	Life Technologies	14175145
Propidium Iodide - 1.0 mg/mL Solution in Water	Life Technologies	P3566
Mouse on Mouse blocking reagent	Vector Labs	MKB-2213
Isofluorane	Vet One	502017
Hematoxylin	Ricca Chemical Company	3530–16
Eosin	EMD-Millipore	588X-75
DDIT3 DsiRNA 13.9	Integrated DNA Technologies	mm.Ri.Ddit3.13.9
DDIT3 DsiRNA 13.1	Integrated DNA Technologies	mm.Ri.Ddit3.13.1
Matrigel	BD Biosciences	356234
Horse Serum	GIBCO-Invitrogen	26050088
0.25% Trypsin EDTA	GIBCO-Invitrogen	25200072
Opti-MEM	Thermo Fisher	31985088
Lipofectamine RNAiMAX	Lifetech	S-006–100
Fibroblast Growth Factor basic	GIBCO-Invitrogen	PHG0263
Penicillin Streptomycin	GIBCO-Invitrogen	15640055
QIAzol	QIAGEN	79306
Hoechst 33342	Thermo Fisher	H3570
Tween-20	Sigma Aldrich	P1379
Triton X-100	Sigma Aldrich	T8787
SYBR Green PCR MasterMix	Thermo Fisher	4309155
PrimeTime Mouse GAPDH Primer	Integrated DNA Technologies	Mm.PT.39a1
PrimeTime Mouse DDIT3 Primer	Integrated DNA Technologies	Mm.PT.58.30882054
TissuePlus O.C.T Compound	Fisher Scientific	23–730-571
Dako Fluorescence mounting media	Agilent	S3023
Hoechst 33342	Thermo Fisher	H3570
4’,6-Diamidino-2-Phenylindole, Dihydrochloride (DAPI), FluoroPureTM grade	Invitrogen	D21490
MitoTracker Deep Red	Thermo Fisher	M22426
Corning CellTak	Fisher Scientific	C354240
Gelatin from porcine skin, type A	Sigma Aldrich	G2500
All-trans Retinoic Acid	Fisher Scientific	AC207341000
Dimethyl sulfoxide, Hybri-Max sterile-filtered BioReagent	Sigma Aldrich	D2650
Critical Commercial Assays
Smart-Seq v4 Ultra Low Input RNA Kit	Clontech	634888
Nextera XT DNA Library Preparation Kit	Illumina	FC-131–1024
NextSeq 500/550 High Output Kit v2 (150 Cycles)	Illumina	FC-404–2002
Nextera DNA Library Preparation Kit	Illumina	FC-121–1030
SuperScript III First-Strand Synthesis Kit	Thermo Fisher	18080051
QIAGEN miRNeasy Micro Kit	QIAGEN	217084
Satellite Cell Isolation Kit, mouse	Miltenyi	130–104-268
Deposited Data
ATACseq and RNaseq Datasets	This Manuscript	https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE121589
Experimental Models: Organisms/Strains
C57BL/6 wild-type female mice (3–4 months)	Charles River Labs	Strain 027
C57BL/6 wild-type female mice (20–24 months)	NIA	Strain 027
Pax7CreER/+;Rosa26mTmG/+ female mice (3 months)	Jackson Labs	Stock #017763 crossed with stock #007676
Pax7CreER/+;Rosa26nTnG/+ female mice (5 months)	Jackson Labs	Stock #017763 crossed with stock #023537
Software and Algorithms
Bowtie2	Langmead and Salzberg, 2012	http://bowtie-bio.sourceforge.net/bowtie2/index.shtml
Samtools	Li and Dewey, 2011	http://samtools.sourceforge.net/
STAR	Dobin et al., 2013	https://pubmed.ncbi.nlm.nih.gov/23104886/
RSEM		https://pubmed.ncbi.nlm.nih.gov/21816040/
limma	Ritchie et al., 2015	https://pubmed.ncbi.nlm.nih.gov/25605792/
Voom		https://pubmed.ncbi.nlm.nih.gov/24485249/
SVA		https://pubmed.ncbi.nlm.nih.gov/22257669/
DPGP clustering		https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005896
DAVID		https://pubmed.ncbi.nlm.nih.gov/19131956/
MISO		https://idp.nature.com/authorize?response_type=cookie&client_id=grover&redirect_uri=https%3A%2F%2Fwww.nature.com%2Farticles%2Fnmeth.1528
ENCODE ATAC-seq processing pipeline	Lee et al., 2019	https://github.com/ENCODE-DCC/atac-seq-pipeline
cutadapt	Martin, 2011	https://github.com/marcelm/cutadapt/blob/master/CITATION
MACS2		https://genomebiology.biomedcentral.com/articles/10.1186/gb-2008–9-9-r137
Bedtools	Quinlan and Hall, 2010	https://pubmed.ncbi.nlm.nih.gov/20110278/
GREAT		https://pubmed.ncbi.nlm.nih.gov/20436461/
ChromHMM	Ernst and Kellis, 2012	https://pubmed.ncbi.nlm.nih.gov/22373907/
HOMER		https://pubmed.ncbi.nlm.nih.gov/20513432/
DeepTools	Ramírez et al., 2014	https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4987876/
R v3.4.1, v.3.5.0	The R Foundation for Statistical Computing	https://www.r-project.org/
Other
Bioinformatics analyses code	This manuscript	https://github.com/annashcherbina/nobel_lab_projects/tree/master/age
H3K4me3 ChIPseq−Satellite Cells	Liu et al., 2013	GEO:GSM1148118, GEO:GSM1148110
H3K4me3 ChIPseq−Myoblasts	ENCODE	ENCODE:ENCFF360QRN
H3K27me3 ChIPseq−Satellite Cells	Liu et al., 2013	GEO:GSM1148119, GEO:GSM1148111
H3K27me3 ChIPseq−Myoblasts	ENCODE	ENCODE:ENCFF569LDY
MyoD1 ChIPseq	ENCODE	ENCODE:ENCFF423NWT
MyoG ChIPseq	ENCODE	ENCODE:ENCFF980DKG
CTCF ChIPseq	ENCODE	ENCODE:ENCFF297NKN

EXPERIMENTAL MODEL AND SUBJECT DETAILS

Animals

Young (3–4 months) and aged (20–24 months) C57BL/6 wild-type female mice were obtained from Charles River Breeding Laboratories or from a breeding colony at the University of Michigan (UM). Male 5 month old Pax7CreER/+;Rosa26nTnG/+ (nTnG) mice were obtained from a breeding colony at the University of Michigan. All mice were fed normal chow ad libitum and housed on a 12:12 hour light-dark cycle under UM veterinary staff supervision. All procedures were approved by the Institutional Animal Care and Use Committee and were in accordance with the U.S. National Institute of Health (NIH).

Myogenic Progenitor Cells (C2C12s)

C2C12s are an immortalized mouse myoblast cell line commercially available at ATCC (CRL-1772). They were maintained in culture at 37°C and 5% CO2 using F10 supplemented with 20% FBS, 1% Pen Strep, and 0.02ug/mL bFGF. Media was replenished every 2–3 days, and cells were passaged using 0.25% Trypsin EDTA once they reached 70%–80% confluence.

Animal and Injury Model

C57BL/6 wild-type female mice were obtained from Charles River Breeding Laboratories or from a breeding colony at the University of Michigan (UM). All mice were fed normal chow ad libitum and housed on a 12:12 hour light-dark cycle under UM veterinary staff supervision. All procedures were approved by the Institutional Animal Care and Use Committee and were in accordance with the U.S. National Institute of Health (NIH). Young female mice (3–4 months) and aged female mice (20–24 months) were randomly assigned to one of five groups: uninjured, day 1, day 3, day 5, and day 7 injured (n = 4 per group). To induce skeletal muscle injury, mice were first anesthetized with 2% isoflurane and bilaterally administered a 1.2% barium chloride (BaCl2) solution injected intramuscularly into several points of the tibialis anterior and both heads of the gastrocnemius muscles for a total of 80 µL per hindlimb. We used an uninjured set of mice as our controls because of the previous observation (Rodgers et al., 2014) that the contralateral control after muscle injuries adopts an activated state (called GAlert) in muscle stem cells that is unique compared to uninjured muscle stem cells.

METHOD DETAILS

Satellite Cell Isolation via Fluorescence-Activated Cell Sorting

For tissue collection, mice were anesthetized with 2% isoflurane, then euthanized with cervical dislocation, bilateral pneumothorax, and removal of their heart. Tissues were quickly excised in either biohazard containment or a surgical room. Hind limb muscles (tibialis anterior and gastrocnemius) of control and experimental mice were dissected using sterile surgical tools and placed in Petri dishes containing ice-cold PBS. To achieve adequate MuSC yield for downstream analysis, TA and Gas from both legs were pooled, then minced using surgical scissors and transferred into 50 mL conical tubes containing 20 mL of digest solution (2.5 U/mL Dispase II and 0.2% [~5,500 U/mL] Collagenase Type II in DMEM) per mouse. Samples were incubated on a rocker placed in a 37oC incubator for 70 min with manual pipetting the solution up and down to break up tissue every 30 minutes using an FBS-coated 10 mL serological pipette. Once the digestion was completed, 20 mL of F10 media containing 20% heat inactivated FBS was added into each sample to inactivate enzyme activity. The solution was then filtered through a 70 µm cell strainer into a new 50 mL conical tube and centrifuged at 350xg for 5 min at 4oC. The supernatant was discarded and the pellets were resuspended in a total of 6 mL of staining media (2% heat inactivated FBS in Hank’s Buffered Salt Solution - HBSS). The single cell suspension was divided into separate FACS tubes and centrifuged at 350xg for 5 min at 4oC. The supernatants from each FACS tube were discarded and the cell pellets were resuspended in 200 µL of staining media and antibody cocktail containing Sca-1:APC (1:400 dilution), CD45:APC (1:400), CD11b:APC (1:400), Ter119:APC (1:400) and CD29/β1-integrin:PE (1:200) and CD184/CXCR-4:Biotin (1:100), then incubated on ice for 30 minutes in the dark. Following incubation in primary antibodies, samples were diluted with 2mL staining solution per sample, centrifuged at 350xg for 5 min at 4°C, and supernatant discarded. Samples were then re-suspended in a staining buffer containing PECy7:Streptavidin (1:100), and incubated on ice for 20 min in the dark. After incubation, samples were again diluted in 2mL staining buffer, centrifuged at 350xg for 5 min at 4°C, supernatant discarded, then re-suspended in 200uL staining solution. Prior to sorting, cells were filtered through 35um cell strainers, and 1ug of propidium iodide (PI) stain was added in each experimental sample. Cell sorting was done using a BD FACSAria III Cell Sorter (BD Biosciences, San Jose, CA) and APC-PI-PE+PECy7+ MuSCs were sorted into either ice-cold staining solution for immediate processing or into Trizol and snap frozen for later use. Purity of the enriched MuSC population was validated using immunofluorescent labeling of Pax7 (Figure S1B). Coverslides were coated with 22.4ug/mL CellTak in PBS for 20 minutes at room temperature, then cells were seeded and allowed to adhere for 45 minutes at room temperature before being fixed with 4% paraformaldehyde in PBS for 20 minutes at room temperature. Cells were stained and imaged at 10X magnification as described below, using mouse anti-Pax7 (1:10), anti-mouse AF555 (1:300), and Hoechst (1ug/mL). Pooling of MuSCs from both TA and Gas muscles was required to achieve adequate amounts of RNA or DNA for subsequent sequencing-based assays, which prevented our ability to discern the muscle of origin.

Chromatin Accessibility Assay and Sequencing Library Preparation

At least 10,000 MuSCs were centrifuged at 500xg for 5 min at 4°C in a fixed angle centrifuge. After removing the supernatant, the cell pellets were resuspended in 50uL ATAC-Resuspension Buffer (RSB) (500uL 1M Tris-HCl, 100uL 5M NaCl, 150uL 1M MgCl2 in 49.25mL sterile water) with 0.1% NP40, 0.1% tween-20, and 0.01% Digitonin. Samples were pipetted up and down three times to mix, incubated 3 minutes on ice, then diluted in 1mL of ice-cold ATAC-RSB containing 0.1% Tween-20 and centrifuged at 500xg for 10 min at 4°C to pellet the nuclei. The supernatants were carefully removed and the pellets were resuspended in 25 µL of transposase reaction mix (12.5 µL 2x TD buffer, 1.25 µL TN5 transposase (Illumina) and 11.25 µL nuclease-free water). Transposase reactions were carried out by incubating samples at 37°C for 30 min under mild agitation (300 rpm on a Thermo-mixer C, Eppendorf). Once the incubation was completed, sample tubes were placed on ice and the transposed DNA fragments from each sample were purified using a QIAGEN MinElute PCR Purification Kit following manufacture’s protocol. Purified DNA fragments were then amplified for 13 cycles using barcoded PCR primers and NEB Next High Fidelity 2x PCR Master Mix (New England Bio Labs) on a thermal cycler. Double concentrated Ampure beads were used to purify transposed DNA amplicons. The molarity of each DNA library was determined (Agilent 2100 Bioanalyzer), pooled into a single tube and sequenced on a NextSeq (Illumina) using 76-bp paired-end reads.

Bulk mRNA Isolation and Sequencing Library Preparation

MuSCs sorted directly into Trizol were thawed at room temperature, and RNA was extracted using a QIAGEN miRNeasy Micro Kit as per manufacturer’s instructions. RNA concentration and integrity were measured with a Nanodrop spectrophotometer (Nanodrop 2000c) and Bioanalyzer (Agilent 2100). 1–10 ng of high-quality RNA (RIN > 8) was used to produce cDNA libraries using the SmartSeq v4 protocol (Clontech) as per the manufacturer’s instructions. cDNAs were prepared into sequencing libraries using 150 pg of full-length cDNA amplicons (Nextera XT DNA Library Preparation Kit, Illumina) with dual index barcodes as per manufacturer’s instructions. Barcoded cDNA libraries were pooled into a single tube and sequenced on a NextSeq (Illumina) using 76-bp single-ended reads.

Histology and Immunohistostaining

A separate cohort of mice was used for histological evaluation of skeletal muscle injury and regeneration. Young (3–4 months) and old (20–24 months) female mice (n = 4 per group) received bilateral intramuscular injection of the tibialis anterior muscles with a 1.2% BaCL2 solution (40 µL per muscle) while under isoflurane anesthesia. At day 7 post BaCl2 injury, mice were euthanized via bilateral pneumothorax and cervical dislocation, and tibialis anterior muscles were rapidly dissected. Tibialis anterior muscles were embedded in optimum cutting temperature (OCT) compound, rapidly frozen in isopentane, cooled on liquid nitrogen and stored at −80°C until further analysis.

Muscle tissue cross sections (10 µm) were cut in a cryostat at −20°, adhered to Superfrost Plus microscope slides (Fisher Scientific) and air-dried at room temperature. Slides were fixed in 100% acetone for 10 minutes at −20°C, air-dried at room temperature and stained with hematoxylin and eosin (H&E) using routine procedures. Sections for immunofluorescence staining were rehydrated in PBS for 5 minutes and blocked overnight at 4°C in M.O.M. blocking reagent. Sections were then incubated overnight at 4°C with a cocktail of primary antibodies against embryonic myosin (1:20) and laminin (1:200). The following morning, slides were washed three times for 5 minutes in PBS and incubated for 1 hour at room temperature with Goat Anti-Mouse IgG1 Alexa Fluor 488 (1:500) and Goat Anti-Rabbit Alexa Fluor 555 (1:500) secondary antibodies. Nuclei were counterstained with 4’,6-Diamidino-2-Phenylindole Dihydrochloride (DAPI) (1 µg/mL). Slides were washed three times for 5 minutes each in PBS and then mounted with coverslips using Dako Fluorescence mounting media. Fluorescent and bright field images were acquired using a Nikon A1 confocal and Olympus BX51 wide field microscope. Regenerating myofibers (central nuclei and/or eMHC+ cytoplasm) were quantified in a single 20x field of view from the core of the injured region using ImageJ software. Differences between young and aged mice for regenerating myofibers and regenerating myofiber cross-sectional areas were determined by Student’s t tests where statistical significance was set at p < 0.05.

Chromatin Immunofluorescence

Black, 96-well flat clear-bottom plates were coated with 22.4ug/mL CellTak in PBS for 20 minutes at room temperature, followed by three quick rinses with distilled water. FACS-enriched MuSCs were then seeded at a density of 10,000 cells/well and allowed to adhere for 45 minutes at room temperature. Next, cells were washed twice with PBS, fixed with ice-cold 100% methanol for 20 minutes at −20°C, and rinsed three times with PBS. Next, cells were blocked in Odyssey blocking buffer (OBB; LI-COR Biosciences) for 1 hour at room temperature and incubated overnight at 4°C with the following primary antibodies: Histone H3; H3K27me3; and H3K9me3. The following day, cells were washed three times with PBS supplemented with 0.1% Tween-20 (PBS-T) and incubated at room temperature with the secondary rabbit antibody Alexa Fluor 647 for 1 hour. Cells were then washed once with PBS-T and once in PBS. Hoechst 33342 was used to stain nuclei, after which cells were washed twice with PBS and imaged with a 20x objective using the ImageExpress Micro Confocal system. A total of 49 sites were imaged per well. Background subtraction was performed using ImageJ, whereas image segmentation was performed with CellProfiler. Single-cell and population-average image analysis and quantification was performed using MATLAB 2017b software. Immunofluorescence images were assembled using ImageJ and Adobe Illustrator.

Treatment with all-trans Retinoic Acid

Aged and young MuSC isolation.

Aged MuSCs were isolated from the hindlimb muscles of two C57BL/6 wild-type mice (26 months, 1 male and 1 female) via FACS as described above. Young MuSCs were isolated from two 5-month-old male Pax7CreER/+;Rosa26nTnG/+ (nTnG) mice following 5 daily intraperitoneal injections of 20 mg/mL tamoxifen in corn oil (75 mg/kg body weight) and 5 days of recovery from injections. The same dissection, digestion, and filtration techniques were used as above. Cells were incubated with 200 µL of 1 µg/mL DAPI in staining solution at room temperature for 10 minute (protected from light) to stain for viability. Cell suspensions were passed through a 35 µm cell strainer prior to sorting on a Sony MA900 cell sorter, and DAPI-/tdTomato-/GFP+ MuSCs were collected into cold staining solution for immediate processing.

Cell culture.

The wells of a 96-well plate were coated with sterile 0.5% gelatin (in ddH2O) solution for an hour at room temperature and then aspirated and allowed to dry. 100mM aliquots of all-trans retinoic acid (ATRA) were prepared by dissolution in cell culture-grade dimethyl sulfoxide (DMSO). ATRA was dissolved in myoblast media (F10 with 20% FBS, 1% Pen Strep, and 0.02ug/mL bFGF) for a final ATRA concentration of 100 nM and 0.1% DMSO (v/v). Control media was prepared by diluting DMSO in myoblast media (0.1% v/v). Aged and young MuSCs were re-suspended in myoblast containing either ATRA (+RA) or DMSO (-RA). Cells were seeded in gelatin-wells at a density of approximately 6,250 cells/cm2. Media containing either ATRA or DMSO was replenished every 24 hours for 3 days.

Immunofluorescence (IF) staining

Cells in wells were fixed with 100% methanol for 5 minutes at room temperature. After washing with PBS, cells were permeabilized and blocked with 0.3% Triton X-100 and 1% BSA in PBS, incubated overnight at 4°C with a combination of a 1:100 dilution of AF488-conjugated anti-Pax7 antibody and a 1:50 dilution of AF647-conjugated anti-MyoD antibody in 0.2% BSA in PBS. Cells were incubated overnight at 4°C with a 1:50 dilution of PE-conjugated anti-Ki67 antibody in 0.2% BSA in PBS. Each staining combination was performed in duplicate for each treatment condition and age. Following overnight incubation with antibodies, cells were washed 3 times with PBS. Nuclei were counterstained with DAPI (1 µg/mL) for 10 minutes at room temperature. Cells were washed in PBS a final 3 times and left covered in 100 µL of PBS during imaging. 20x magnification images were acquired on a Zeiss Axio Vert.A1 inverted microscope with a Colibri 7 LED light source and an AxioCam MRm camera. Images were subsequently analyzed in Fiji. ROIs were generated by thresholding on the DAPI image to identify nuclei. The average fluorescent intensities of each stain within these ROIs were recorded. The average Pax7, MyoD, and Ki67 signals were compared between +RA and -RA cells for each age using a two-sample t test with the significance level set to α = 0.05. Statistical tests were performed in R and plots were generated using the ggplot2 package in R.

In vitro DDIT3 Knockdown

C2C12s.

C2C12s were seeded at a density of 50,000 cells per well in a 12-well plate in myoblast media (F10 with 20% FBS, 1% Pen Strep, and 0.02ug/mL bFGF). After 24 hours, the media was replaced with myoblast media without antibiotics containing 15nM DDIT3 13.1 DsiRNA and 15nM DDIT3 13.9 DsiRNA encapsulated in RNAiMAX lipid droplets. After 72 hours, transfection media was replaced with differentiation media (DMEM with 5% horse serum and 1% Pen Strep). Cells were incubated for 72 hours in differentiation media prior to being fixed in 4% paraformaldehyde for immunostaining or lysed in trizol for qPCR.

Primary cells.

Mice were euthanized and their hindlimb muscles extracted and digested into a single cell suspension. Satellite cells were enriched using the Miltenyi mouse satellite cell isolation kit following the manufacturer’s protocol. Tissue culture dishes were coated with 10% Matrigel in DMEM as previously published (Motohashi et al., 2014). Enriched satellite cells were seeded onto the Matrigel-coated plates in myoblast media (F10 with 20% FBS, 1% Pen Strep and 0.02ug/mL bFGF). Media was replaced every 48 hours until cells expanded to 50% confluency, at which point they were passaged as previously published and seeded onto a Matrigel-coated 12-well plate. 24 hours after seeding, 0.15nM DDIT3 13.1 DsiRNA and 0.15nM DDIT3 13.9 DsiRNA were encapsulated in RNAiMAX lipid droplets according to manufacturer’s protocols and delivered to cells. Knockdown efficiency after 72 hours in transfection media was validated using quantitative real-time PCR (Figure S4). Post DDIT3 knockdown, the transfection media was replaced with differentiation media (DMEM with 5% horse serum and 1% Pen Strep). Cells were incubated in differentiation media for 72 hours prior to lysing or fixing in 4% paraformaldehyde in PBS.

Immunostaining

Immediately prior to fixation, cells were incubated for 45 minutes in warmed DMEM containing 250nM Mitotracker Deep Red, then washed 3 times with warmed PBS and fixed for 15 minutes in warmed 4% paraformaldehyde in PBS at 37°C. After three quick washes with PBS, cells were permeabilized with 0.1% Triton X-100 and blocked with 1% BSA, 0.1% Tween-20 and 22.52mg/mL glycine in PBS. After blocking, cells were incubated with primary antibodies (1:10 dilution of anti-MYH3) overnight at 4°C followed by secondary antibodies (1:500 dilution of AF555 anti-rabbit) overnight at 4°C. Nuclei were stained with Hoechst 33342. Immunolabeled cells were imaged on a Zeiss epifluorescent microscope using a 10X objective. The fusion index was automatically calculated using MATLAB as the ratio of nuclei within myofibers containing more than 2 nuclei divided by the total number of nuclei per image. The knockdown and differentiation were performed in at least biological triplicates and technical duplicates. Statistical comparison between groups was performed using a two-sided, two sample Student’s t test assuming equal variances (n between 25 and 41 for primary cells harvested from aged mice and between 13 and 15 for primary cells harvested from young mice). P values below 0.05 were considered significant.

Real-Time PCR

RNA was extracted from lysed cells using the QIAGEN miRNeasy micro kit according to manufacturer’s protocol. RNA quality and concentration were determined using a Nanodrop. Within one week, cDNA was synthesized using SuperScript III cDNA Synthesis Kit according to manufacturer’s protocol, and quality was determined using a Nanodrop. 80–100ug cDNA template was plated in triplicate along with SYBR Green PCR MasterMix and 500nM PCR primer, then cycled 40 times starting at 95°C for 10 s followed by 60°C for 30 s on a CFX96 Real-Time thermocycler. Gene expression was quantified using the ΔΔCt method. Statistical significance was determined using a two-sided, two-sample t test for MyoG expression assuming equal variances, and p values below 0.05 were considered significant. Each condition was performed in biological and technical triplicates, and technical triplicates were averaged prior to statistical analysis (n = 3 for each condition).

RNA-Seq Data Processing and Analysis

Double expression quantification.

Double-stranded RNA-seq data was aligned to the mm10 reference genome with the STAR algorithm (Dobin et al., 2013; Li and Dewey, 2011)

STAR–genomeLoad NoSharedMemory–outFilterMultimapNmax 20–alignSJoverhangMin 8–alignSJDBoverhangMin 1–outFilterMismatchNmax 999–outFilterMismatchNoverReadLmax 0.04–alignIntronMin 20–alignIntronMax 1000000–alignMatesGapMax 1000000–outSAMunmapped Within–outFilterType BySJout–outSAMattributes NH HI AS NM MD–outSAMtype BAM SortedByCoordinate–quantMode TranscriptomeSAM–sjdbScore 1–limitBAMsortRAM 60000000000–twopassMode Basic–twopass1readsN −1

Differential gene expression.

Differentially expressed genes between Young and Old samples at each time point (days = 0,1,3,5,7) were identified by limma (Ritchie et al., 2015) analysis in R. The Expected Counts from RSEM were transformed to counts per million using the voom (Law et al., 2014) R package with a design formula: Count~Day + Age, with Day = {0,1,3,5,7}, and Age = {Young,Old}. Surrogate variable analysis was performed with the SVA package (Leek et al., 2012) using a null model of voom$E~1, and a design matrix of voom$E~Day+Age. Contributions from the surrogate variables were quantified and removed from the voom$E data matrix. An lmfit analysis was performed to check for statistically significant (p < 0.05) associations between surrogate variables and the protected variables of Day and Age. The removeBatchEffect function in limma was used to remove contributions of surrogate variables from the data, while protecting the variables of Age and Day. Pairwise Pearson and Spearman correlation values were computed between all sva-corrected replicates. Any replicate that had r < 0.9 with other replicates for a given sample was excluded from further analysis.

The resulting data was analyzed with limma, using a design matrix of voom$E ~Sample. All pairwise contrasts were examined to identify differentially expressed genes between Old versus Young at time t, Old at time t versus Old at time (t-1), Young at time t versus Young at time (t-1). Thresholds of p.adjusted < 0.05 and log-fold-change > = 1 were used to call differential genes.

Time-series Clustering of Differential Genes

The differential genes underwent time-series clustering with the DPGP (Dirichlet process Gaussian process mixture models) algorithm (McDowell et al., 2017). To generate inputs for the algorithm, replicates were averaged, and the fold change in TPM of Old versus Young was computed for each set of averaged replicates.

Pathway and GO Term Enrichment Analysis

The clusters generated by the DPGP algorithm were analyzed with DAVID (Huang et al., 2009a, 2009b)

Transcription Factor Expression Analysis

A similar approach was used to identify differentially expressed transcription factors. A list of known GRCm38 transcription factors was downloaded from Animal Transcription Factor Database (Zhang et al., 2012). Log2-fold change in FPKM was calculated for MDX injured versus control, as well as wild-type injured versus control. abs(log2-fold) change cutoff of 1 was used to determine significance.

ATAC-Seq Data Processing and Analysis

Generation of a merged peak set across ATAC-seq samples.

The ATAC-Seq samples were analyzed with the ENCODE ATAC-seq processing pipeline (https://github.com/ENCODE-DCC/atac-seq-pipeline, version 0.3.4) (Maher, 2012). The cutadapt algorithm (Martin, 2011) was used to trim adapters, and the Bowtie2 (Langmead and Salzberg, 2012) aligner was used to align the reads to the mm10 index. Duplicates were then removed from the aligned reads, and the MACS2 (Zhang et al., 2008) peak caller was used to call peaks from the aligned ATAC-seq samples. The naive overlap peak set from all replicates for a given sample was used. Such peak sets from all samples were concatenated and merged using the bedtools merge command to produce a master peak set across all samples. The read counts for each sample at each peak were obtained by running bedtools coverage on the shifted de-deplicated tagAlign files generated by the pipeline.

Differential chromatin accessibility.

The resulting counts matrix (238,563 peaks x all sample replicates) was analyzed with limma voom to produce differential peak sets between aged and young at each time point. The analysis protocol matched that for RNA-seq data (see Differential Gene Expression) section above, with equivalent thresholds for differential peaks.

Time-series Clustering of Differentially Accessible Chromatin Regions

All differential peaks underwent DPGP clustering, as described above for RNA-seq data. The log2 fold change in counts per million in young samples versus ages samples were calculated for each differential peak. These fold change values for days 0, 1, 3, 5, 7 were used as inputs for DPGP clustering. 10 significant clusters were identified.

GO Term and Pathway Enrichment

Differential peaks between aged and young for each time point underwent analysis with the GREAT algorithm (McLean et al., 2010), using default parameters. A parallel approach to identify enriched GO terms and pathways was also utilized. In this approach, results from the GREAT analysis were used to map peaks with associated genes (incorporating both distal and proximal associations returned by GREAT). The gene sets were then analyzed with DAVID to identify enriched GO terms and Kegg pathways (Benjamini-Hochburn corrected P value < 0.01).

In addition to GREAT analysis of differential peaks at each time point, the peaks for each of the 10 clusters from DPGP clustering (described above) underwent analysis with GREAT and GO. The enrichment analysis on peak clusters was performed to allow for analysis of more complex peak trajectories compared to enrichment at a single time point.

Data Visualization

The p value bigWig signal tracks generated by the ENCODE ATAC-seq processing pipeline were visualized with the WashU browser (Zhou and Wang, 2012):

(http://epigenomegateway.wustl.edu/browser/?genome=mm10&datahub=http://mitra.stanford.edu/kundaje/annashch/age/datahub.pval.json&tknamewidth=150).

Differential Chromatin State Distributions

The core mark 15-state chromatin state model (Ernst and Kellis, 2012) built on skeletal muscle tissue (ENCODE cell id E107) was used to identify the chromatin state distributions for differential peaks (Aged versus Young) at each time point. For each time point, the fraction of peaks in each of the 15 chromatin states was determined.

Motif Enrichment Analysis

The HOMER algorithm (Heinz et al., 2010) was used to identify enriched motifs in Aged versus Young (and vice versa) at each time point. For each time point, the foreground to HOMER consisted of differentially accessible peaks in aged versus young (and Young versus Aged in a parallel run), while the background consisted of all differential peaks merged across time points. The purpose of this analysis was to identify motifs that were enriched only at a particular time point.

In a separate run, the foreground was kept the same, but the background was the default reference mm10. The purpose of this analysis was to identify motifs that were enriched in the Aged or Young samples relative to the genome background.

Comparison of TF expression and corresponding motif enrichment across time

Differential TF motifs from HOMER (see above) were cross-referenced with bulk RNA-seq TPM values at each time point in young and aged samples. For days 0, 1, 3, 5, 7, the log2 fold change of young TPM versus aged TPM was compared to the −10log10(corrected P value) from HOMER to generate a scatterplot with R ggplot2 (Wickham, 2016). The TF’s were color-coded on the scatterplot based on TF family membership. The mapping of TF to corresponding TF family was accomplished in accordance with the description below (Heatmap of TF family enrichment in Young versus Aged).

Footprint analysis

Footprint analysis was performed for the Myod transcription factor motifs within the ATAC-seq datasets. The HOMER scanMotifGenomes.pl command was executed on the mm10 genome with the HOMER pwm for HOMER Myod motifs:

>AGCAGCTGCTGC MyoD(bHLH)/Myotube-MyoD-ChIP-Seq(GSE21614)/Homer 7.249827 −2.493974e+04 0 49999.2,13125.0, 15210.6,14377.0,0.00e+00

0.432 0.137 0.291 0.140

0.400 0.102 0.455 0.042

0.001 0.997 0.001 0.001

0.966 0.001 0.001 0.032

0.001 0.014 0.984 0.001

0.001 0.996 0.002 0.001

0.093 0.001 0.001 0.905

0.001 0.001 0.997 0.001

0.003 0.518 0.079 0.400

0.123 0.268 0.093 0.516

0.187 0.234 0.377 0.202

0.165 0.373 0.197 0.265

The found motifs were intersected with ATAC-seq optimal overlap peaks to identify peaks with TF footprints. Motif positions overlapping the peaks were extended to 200 bp. The bam mm10 alignments from the ENCODE ATAC-seq pipeline were corrected for ATAC enzymatic bias using the TOBIAS toolkit (Bentsen et al., 2019) (ATACorrect). The resulting corrected bigwigs were used with the TOBIAS command PlotAggregate to generate log-transformed footprints for the Myod regions.

Comparison with external ChIP-seq datasets

Peak call files from TF-ChIP-seq and histone ChIP-seq datasets were downloaded from studies examining satellite cells and myoblasts as indicated below:

Bedtools intersect was used to determine the number of overlapping peaks between the optimal overlap ATAC-seq peak sets and each of the samples above. A Wilcoxon rank enrichment test was performed to determine whether the change in fraction of overlapping peaks with each dataset in the table above was significantly different between young and aged samples at each of the 5 time points (d0-d7).

Derivation of reaction fluxes using gene expression data and genome-scale metabolic modeling

Gene expression data was used as input to derive reaction flux information from the human genome-scale metabolic model (RECON1) (Duarte et al., 2007; Shlomi et al., 2008) using a modeling approach detailed in Shen et al. (2019a, 2019b) This approach maximizes the flux through the metabolic reactions that are upregulated in a condition while minimizing flux through those reactions that are associated with downregulated genes. All of the data processing steps described henceforth were carried out using MATLAB R2018b (https://www.mathworks.com/products/matlab.html). First, expression data was normalized across each gene. Next, normalized expression data between the old and young groups was compared to attain a list of significantly expressed genes using a threshold of p value < 0.05. Based on this list, up- and downregulated genes were determined for each sample based on a z-score threshold of 1.5 and −1.5, respectively. The lists of up- and downregulated genes for each sample were then overlaid onto the RECON1 model based on gene-protein-reaction annotations in the model. Finally, reaction flux data was generated using a linear optimization version of the iMAT algorithm (Duarte et al., 2007; Shlomi et al., 2008) with the following inputs: the RECON1 model, the list of up- and downregulated genes, and the recommended values for the optional parameters (rho = 1E-3, kappa = 1E-3, epsilon = 1, mode = 0).

Visualization (R)

Histogram and jitter plot.

The visualization methods for histogram and jitter plot generation (Figure 2) were carried out through R software (https://www.r-project.org/). For both plots, the z-score was calculated for each reaction by comparing flux values between aged and young groups using a paired t test. Negative and positive extreme z-scores values were floored to −5 and 5, respectively, and invalid values were defaulted to zero with a p value = 1. Significant changes in reaction flux were defined by a p value < 0.05, and these reactions were presented along with the plots (excluding transport reactions).

QUANTIFICATION AND STATISTICAL ANALYSIS

Statistical details, including sample size (n), what n represents, and statistical test used can be found in the figure legends. In most cases, sample size was great enough to assume normality based on the central limit theorem, and parametric statistical tests were used. Two-sided tests were employed for more conservative calculations of significance. The threshold for statistical significance was set at p < 0.05. Unless otherwise stated, data in bar graphs are expressed as mean ± standard deviation. A combination of MAT-LAB_R2019b, R (v3.6), and Prism GraphPad were used for statistical analysis.

	Satellite Cells	Myoblasts
H3K4me3	Liu et al., 2013	ENCODE
	• d0 Aged (GSM1148118)	• ENCFF360QRN
	• d0 Young (GSM1148110)
H3K27me3	Liu et al., 2013	ENCODE
	• d0 Aged (GSM1148119)	• ENCFF569LDY
	• d0 Young (GSM1148111)
MyoD1		ENCODE
		• ENCFF423NWT
MyoG		ENCODE
		• ENCFF980DKG
CTCF		ENCODE
		• ENCFF297NKN