Isolation of pure populations of MuSCs and niche cells by fluorescence-activated cell sorting (FACS)
Hindlimb muscles were dissected from young (4–6 weeks old) or aged (22–26 months old) C57BL/6J (Jackson Laboratory, 000664) mice and minced until no visible tissue chunks were visible. Muscle was then digested in a 15 ml Falcon tube containing 2.4 U/ml Collagenase D (11088882001, Roche), 12 U/ml Dispase II (04942078001, Roche) and 0.5 mM CaCl2 in Ham’s F10 media (11550043, Gibco) for two periods of 30 min each on a shaker in a 37 °C/5.0% CO2 Tissue Culture (TC) incubator. After the first 30 min, the remaining muscle chunks were pelleted by centrifugation at 600 × g for 25 s and the supernatant containing released mononuclear cells was supplemented with 9 ml Fetal Bovine Serum (FBS) (080450, Wisent). The remaining pellet of undigested muscle tissue was triturated, and additional digestion buffer was added for the second 30-min incubation. The combined digested tissue was then filtered through a 40μm cell strainer (C352340, Falcon) and spun down at 600 × g for 18 min. The supernatant was discarded and the cell pellet was dissolved in 500 μl 2% FBS/PBS (v:v) with 0.5 mM EDTA, and filtered once more. For isolation of GFP + MuSCs from Pax7-nGFP mice, 0.5 μg of DAPI (Invitrogen, D3671) was added for cell sorting of the DAPI−/GFP+ population using a FACSAria Fusion cytometer (BD Biosciences). For simultaneous isolation of MuSCs, FAPs and macrophages, the digested tissue was supplemented with FC (fragment crystallizable) block (clone 2.4G2). Antibodies for anti-mouse ITGA7-AF647 (R&D systems, FAB3518R, 1:500), anti-mouse CD31-PE (Invitrogen 12-0311-81, 1:5000), anti-mouse Sca1-FITC (BD Pharmigen, 557405, 1:5000), anti-mouse CD45-BV785 (Biolegend, 103149, 1:5000), anti-mouse CD11b-PerCP-Cy5.5 (Biolegend, 101227, 1:5000), anti-mouse F4/80-APC-Cy7 (Biolegend, 123117, 1:1000), and Hoechst 33342 (H1399, Molecular Probes) were added and the sample was incubated at room temperature (RT) for 15 min with intermittent shaking. The sample was then washed with 5 ml of 2% FBS/PBS (v:v) with 0.5 mM EDTA and spun down at 600 × g for 15 min. The stained cells were then resuspended in 700 μl of 5% FBS/PBS (v:v) with 0.5 mM EDTA and filtered through a 40 μm filter before sorting. Color compensation was done by staining 20 μl of UltraComp eCOMP beads (01-2222-41, Invitrogen) in 200 μl of 2% FBS/PBS (v:v) and calculated using DIVA before sorting on a FACSAria fusion III cytometer (BD Biosciences). In order to maintain the ratios of MuSCs, macrophages and FAPs, all cells were sorted into a single tube and 2500 cells were captured per sample for single-cell RNA-Sequencing using 10X Genomics chromium platform. FACS gating strategies benefited from Pasut et al.56, Farup et al.57, and Low et al.58.
Allogeneic MuSC transplantation into a heterologous niche
NOD-Prkdcem26Cd52Il2rgem26Cd22/NjuCrl (NCG) immunocompromised mice (Charles-River) were irradiated on their hindlimb the day before transplantation. Mice were anesthetized by IP injection of a rodent cocktail composed of ketamine (100 mg/kg)/xylazine (10 mg/kg)/acepromazine (3 mg/kg) and their hindlimb was positioned into the field of irradiation of a Multirad225 irradiation machine (Precision X-Ray) and were given a 18 Gy dose of irradiation as described previously59,60,61. The next day, MuSCs were isolated from Tg(Pax7-EGFP)#Tajb46 mice as described previously, and sorted into a 2% FBS/PBS mix. 10 000-20 000 freshly isolated GFP+ donor MuSCs were immediately injected into the Tibialis Anterior (TA) muscle of the previously irradiated hindlimb of host NCG mice using a Hamilton Syringe, in 2–3 separate injection sites along the TA muscle.
Immunostaining of TA muscles
TA muscles were dissected tendon to tendon and rinsed in cold PBS. The isolated muscles were then fixed in 0.5% Paraformaldehyde in PBS for 2.5 h at 4 °C and then placed in 20% Sucrose (S7903, Sigma) at 4 °C overnight. Muscle samples were frozen in aluminum foil cups of Optimal Cutting Temperature (OCT) compound using liquid nitrogen-chilled isopentane and stored at −80 °C prior to sectioning by the cryostat at 8 μm interval. TA muscle cross-sections on glass slides were circled using a hydrophobic PAP pen. Cross-sections were permeabilized in 0.2% Triton X-100/0.1 M glycine in PBS for 7 min on a shaker, and then washed in PBS. Cross-sections were blocked in 3% (m/v) BSA/10 % Goat Serum solution for at least 1 hat RT in a humid chamber. TA muscle sections were then incubated with primary antibodies for PAX7 (DSHB, AB_528428, 1:100), GFP (Invitrogen, A-11120, 1:1000), and/or Laminin (Sigma, L9393, 1:750) diluted in blocking solution in a humid chamber at 4 °C overnight, then washed 3 times for 10 min with PBT (0.05% Triton X-100 in PBS). Secondary antibodies diluted in blocking solution were added for 1 h at RT in a humid chamber, followed by 3 ×10-min washes with PBT. ProLong Gold Antifade Solution with DAPI (P3695, Invitrogen) was used for mounting prior to imaging.
Analysis of cell purity post FACS isolation
Cells were freshly sorted by FACS into a 1.5 ml Eppendorf tube containing 100 μl of 2% FBS/PBS (v:v). Freshly sorted cells were pelleted by centrifugation at 600 × g and resuspended in 50 μl of PBS. The cells were then plated in a well of a 12-well plate, and excess liquid was allowed to evaporate for around 10 min, or until wells appeared dry. Cells were then fixed using 3.2% PFA/PBS to the well, gently to not dislodge the cells, for 10 min at RT. After washing 3x with PBS, cells were permeabilized with 0.5% Triton X-100 for 30 min at RT. Cells were then washed 2X with 0.3% Triton X-100/PBS and blocked with 0.3%Triton X-100/0.5% BSA/PBS for 1 h at RT. Primary antibodies for PDGFRA (Cell signaling, 31745), F4/80 (Invitrogen, 124801-82), and PAX7 (DSHB, AB_528428) were diluted in blocking solution and added for an overnight incubation at 4 °C in a humid chamber. Next, the cells were washed 3X with 0.3% Triton X-100 before the addition of the secondary antibodies, diluted in blocking buffer for 1 h at RT. Finally, cells were washed 2X with 0.3% Triton X-100 before the addition of PBS with DAPI and visualized on a microscope.
Single-cell RNA-sequencing (scRNA-seq)
MuSCs, macrophages and FAPs were isolated from three young (5 weeks old) and four aged (23.5 months old) mice by FACS, as described previously, and immediately processed for single-cell RNA-sequencing (scRNA-seq) using Chromium Single Cell 3’ Reagent Kits (v3 Chemistry) using 10x genomics technology. In order to maintain ratios, all cells were sorted into a single tube and to aim for a total of 2500 captured cells for sequencing per biological replicate. The viability of cells was 85.5% on average, with a capture range of 2437 to 2835 cells. Quality control and processing steps can be found at https://csglab.github.io/transcriptional_reprogramming_muscle_cells/.
SMART-seq library preparation
MuSCs were sorted by FACS directly into 9 μl of nuclease-free water containing 1 μl SMART-Seq Reaction Buffer, which is composed of 10x lysis Buffer supplemented with RNase inhibitors. The volume was brought to 11.5 μl and libraries were processed using the SMART-Seq HT library preparation kit (634456, Takara Bio) as described47,62. cDNA libraries were then purified using Ampure XP beads (A63880, Beckman Coulter) at a 1:1 (v:v) ratio and quantified using the Quant-iT PicoGreen dsDNA Assay kit (P11496, Invitrogen). 0.15 ng of cDNA in a 1.25 μl volume was used as input for Nextera XT (FC-131-1024, FC-131-2001, Illumina) tagmentation, as described63,64. Sequencing libraries underwent Ampure XP size selection using a 1:0.85 (v:v) ratio, and an aliquot was used to verify size on an agarose gel stained with GelGreen dye (41005, Biotium). Final libraries underwent quality control verifications using a bioanalyzer and were sequenced on an Illumina NextSeq500.
Assay for transposase-accessible chromatin using sequencing (ATAC-seq)
Low input ATAC-seq (Assay for Transposase-Accessible Chromatin using sequencing) was performed according to the OMNI ATAC-seq65 protocol. Briefly, 5000 MuSCs were FACS sorted into the lysis buffer (10 mM Tris-HCl pH 7.5, 10 mM NaCl, 3 mM MgCl2, 0.1% Tween-20, 0.1% NP-40, 0.01% Digitonin) and incubated for 5 min on ice and subsequently 3 min at room temperature (RT). Cells were then washed with 100 μl of wash buffer (10 mM Tris-HCl pH 7.5, 10 mM NaCl, 3 mM MgCl2, 0.1% Tween-20) and spun at 800 g for 10 min. The pellet was resuspended in the transposition mixture at a total volume of 10 μl (5 μl of Tagment DNA (TD) buffer, 3.2 μl PBS, 0.89 μl Tn5, 0.1% Tween-20, 0.01% Digitonin and 0.75 μl water). Transposition was performed for 20 min at 37 °C while shaking the tubes every 5-7 min. The DNA was then purified using column purification as described (Qiagen, QIAquick PCR Purification Kit Cat: 28104). The purified DNA was then PCR-amplified using Q5 DNA polymerase for 12 cycles with the Illumina Nextera XT adapters. The amplified DNA was then size selected and purified with Ampure XP beads at a 1:0.85 (v:v) ratio. The quality control of the final libraries was performed by bioanalyzer and paired-end sequencing was performed on NovaSeq 6000 Sprime paired end (PE 150 bp).
Genomic DNA extraction
Genomic DNA was extracted using phenol/chloroform/isoamyl alcohol solution (Invitrogen #15593-031). Briefly, primary myoblasts derived from in vitro expansion of freshly sorted MuSCs were pelleted at 600 × g for 10 min and snap frozen in liquid nitrogen before use. When all samples were collected, cell pellets were thawed on ice and lysed in 400 μl lysis buffer (40 mM Tris-HCl pH 8.0; 1% (vol/vol) Triton X-100; 0.1% (vol/vol) SDS; 4 mM EDTA (pH 8.0) and 300 mM NaCl). Cell lysates were incubated at 37 °C for 45 min with 10 μl RNase A (Invitrogen, cat#12091-021), followed by a 1 h 15 min incubation at 45 °C with 10 μl of Proteinase K (10 mg/ml Millipore, cat#70663). One volume of phenol/chloroform was then added to each lysate, and samples were vortexed vigorously and centrifuged at 2110 × g for 3 min. The aqueous layer was retained and mixed with 2 volumes of 100% ethanol (EtOH). Samples were snap frozen in liquid nitrogen for 3 min, thawed on ice, and centrifuged for 20 min at 2110 × g. The DNA pellet was washed 3X with 750 μl of 75% EtOH, spinning each time 2110 × g for 12 min. The DNA pellet was then air-dried and resuspended in TE buffer.
Whole-genome bisulfite sequencing (WGBS) library preparation
Genomic DNA was sheared using COVARIS, and libraries were constructed using the NxSeq AmpFREE low DNA Library Kit by Lucigen. Bisulfite conversion was done using the EZ-96 DNA Methylation-GoldTM MagPrep kit (D5042, Zymo Research).
Analysis of differentially methylated regions (DMRs) after whole-genome bisulfate sequencing
Quality check and adaptor trimming of raw reads was performed by trim-galore (Default setting, Phred score = 20)66. Reads were aligned to the mouse reference genome (mm10) using bismark (internal bowtie2 aligner’s default settings and allowing one mismatch)67. Extraction of methylation calls was done by ignoring the first 2 base pairs from the 5’ end of read 2 to avoid known experimental-introduced bias68.
The DMRs between young and aged samples, were identified with DMRcaller’s (R/Bioconductor) function computeDMRsReplicates69. DMRcaller implements a Beta regression to test the difference in methylation levels between conditions in regions of 100 bps with at least 4 CpG sites. Only regions with a minimum methylation proportion difference of 40% are reported. P-values are adjusted for multiple testing using Benjamini and Hochberg’s method (adjusted p-value < 0.01)70.
Bulk RNA-seq data analysis
Transcript-level abundances were imported into R and collapsed into gene counts using the tximport package71. To estimate the effects of the experimental conditions, we used DESeq272 to fit the following gene-wise model: count~ time + age + time:age, where the variable time denotes the effect of engraftment, the variable age denotes the effect of aging, and the interaction variable time:age represents the reprogramming effect on genes in the aged niche (i.e. the residual effect beyond that of engraftment or aging that explains observed expression of a given gene). For the purposes of ranking and visualization, we used the apeglm shrinkage approach73.
Single-cell RNA-seq data analysis
We used Salmon and alevin74,75 to quantify spliced mRNA and unspliced pre-mRNA abundances from scRNA-seq data. To create a spliced-intronic Salmon quantification index, we used eisaR76 to generate a combined FASTA file containing exonic and intronic sequences from the mouse Genocode M24 transcriptome and genome reference files. The introns were defined using the “collapse” option in which transcripts of the same gene are first collapsed by taking the union of all the exonic regions within a gene and labeling the remaining non-exonic parts as introns. In contrast to the alternative “separate” strategy in which intronic regions are extracted for each transcript separately, the collapse strategy treats an ambiguous read to be more likely to have come from an exon of a given transcript than from an intron of alternate transcript that happens to overlap that exonic region76. A 90 bp flanking sequence was added to each side of each extracted intronic region (equal to the length of Read 2 in Chromium Single Cell 3’ Gene Expression library with v3 chemistry) in order to quantify intron-exon reads as unspliced pre-mRNA abundances. To exclude reads coming from intergenic regions of the genome, we adopted the selective alignment (SA) strategy by providing Salmon the complete genome sequence as a decoy sequence77.
The alevin quantifications were imported into R 4.0.0 using the tximeta package. The abundances were split into a spliced and unspliced abundances matrices to be used for downstream analysis such as RNA velocity. We identified and removed low-quality or potentially multiplet cells based on whether they had high percentage of mitochondrial reads (>15%), low library sizes, low number of genes detected, or whether they tend to be observed as outliers with respect to the distributions of spliced reads, unspliced reads, or the fraction of unspliced reads in each cell (see the analysis notebooks for details: https://csglab.github.io/aging_muscle_niche/pages/notebooks.html).
To cluster cells, we used UMAP78 followed by the Walktrap community finding algorithm79 to reduce the dimensionality of the data and pull together similar cell populations into separated clusters which can then be labeled. To annotate the main clusters, we used the expression of Adgre1, Ptprc, and Itgam as markers of macrophages, Pax7 as a marker of muscle stem cells (MuSCs), and Pdgfra as a marker of FAPs. The initial clustering showed the presence of three large main clusters corresponding to the three cell types along with four small populations whose expression profiles suggest them as being Schwan and smooth muscle cell clusters respectively (Supplementary Fig. 1e-g). The four small cell clusters were ignored in downstream analysis. Next, the seven samples were batch normalized and corrected for batch effects using fast mutual nearest neighbors (MNN) correction80. The clustering approach was repeated independently on the MNN corrected space of the MuSCs, FAP, and macrophages respectively in order to determine the internal heterogeneity of each cell type cluster. As a result, MuSCs were partitioned into three clusters, the FAPs into two clusters, and the macrophage cells into one cluster only.
We ran differential gene expression analysis separately for each cell type. We discarded genes that had low counts, that had low mean expression across the subclusters making up the cells of each cell type, and that were characterized by severe batch effects. We used ZINB-WaVE81 to compute gene and cell-specific observational weights that can be used to unlock DESeq272 for the differential analysis of single-cell RNA-seq data. For data of each of the three cell types, we independently fit the following gene-wise model: count~ −1 + cluster_condition + sample, where the variable cluster_condition denotes the mean expression of cells in each cluster and age condition, and the variable sample denotes the batch effect of each sample. We manually created the design matrix such that the sample effects of one of the aged samples and one of the young samples are removed in order to fit a statistical identifiable model. We tested several contrasts to assess several hypotheses of differential mean gene expression across the different subclusters. For the purposes of ranking and visualization, we used the ashr shrinkage approach82 to preserve the size of estimated large LFC and compute s-values (the estimated rate of false sign).
Integration of bulk and single-cell RNA-seq data
Using the moderated LFCs generated by models fit on both RNA-seq data modalities, we are able to examine the pattern of concordance of differential gene expression across the multiple datasets. More specifically, the moderated LFC model effect estimates for the bulk data and the muscle stem cells (MuSCs) single-cell data were joined together such that we retain only genes that were not discarded from the scRNA-seq data. The LFCs were used for plotting Supplementary Fig. 9g.
Geneset Enrichment Analysis (GSEA)
We used the R package hypeR83 to conduct hypergeometric tests to determine which respective set of statistically significant differentially regulated genes are over-represented in the 50 pre-defined hallmark gene sets listed in the Molecular Signatures Database (MSigDB)84. The threshold for significance was set at 10-8 for the s-values across the contrast effects accompanied with an absolute value cutoff for the moderated LFCs between 0.5 and 2 depending on the distribution of effect values for each contrast. The six sets of significant genes for MuSCs are the following: (1) genes that are downregulated in aging both in bulk and single cell, (2) genes that are upregulated in aging both in bulk and single cell (scRNA-seq), (3) genes that are upregulated in MuSC2 with respect to MuSC1 cluster, (4) genes that are upregulated in MuSC1 with respect to MuSC2 cluster, (5) genes that are upregulated in MuSC3 with respect to young cells in both MuSC1 and MuSC2 clusters cluster, (6) genes that are downregulated in MuSC3 with respect to young cells in both MuSC1 and MuSC2 cluster. For FAPs, there were four sets: (1) genes that are downregulated in aging (2) genes that are upregulated in aging (3) genes that are upregulated in FAP 2 with respect to FAP1, (4) genes that are upregulated in FAP1 with respect to FAP2.
ATAC-seq pre-processing analysis
The ATAC-Seq data was processed using the ENCODE ATAC-seq pipeline (https://github.com/ENCODE-DCC/atac-seq-pipeline). The pipeline consists of adapter read trimming, read alignment, post-alignment read filtering and de-duplicating, calling peaks, creating signal tracks, generating quality control reports, and running Irreproducible Discovery Rate (IDR) analysis85. The full specification of the pipeline is detailed on the website (https://github.com/ENCODE-DCC/atac-seq-pipeline). For downstream analysis the generated IDR optimal peaks were used.
Analysis of differentially bound ATAC-seq peaks
Peaks from young and aged conditions were merged from the IDR pipeline. Peaks with summits within 200 bp were combined, a new summit was set in the middle. The peaks’ start and end were redefined with a range of 500 bp from the summit.
We created a count matrix by extracting ATAC-seq read counts from peak and background (+/− 10kbp around the peak) regions for all four samples, peaks and samples correspond to rows and columns respectively86. We used DiffRAC and the count matrix to identify significant differentially accessible peaks (adjusted p-value < 0.05 and LFC > 1)50. Differentially accessible peaks were associated to genes using GREAT87 with the default parameters for the Basal plus extension mode, where a proximal gene regulatory domain is defined as 5kbp upstream and 1 kb downstream of the gene TSS, and a distal domain of 1 Mb.
Analysis of differential mRNA stability
Differential mRNA stability was inferred as previously described50. Briefly, RNA-seq reads were mapped to exonic and intronic regions using annotations acquired from Ensembl GRCh38 version 87. DiffRAC was then used to decouple transcriptional and post-transcriptional effects and infer changes in mRNA stability.
Differential analysis of cell-cell interactions
To identify the cell-cell interactions occurring in young and aged mice, we used the scDiffCom package in R44. To run the software, the Single Cell Experiment object containing the scRNA-seq data was converted to a Seurat object and provided as an argument to the scDiffCom method using the mouse genome setting. We defined the conditions as aged and young mice samples and grouped Emitter (ligand source) and Receiver (receptor source) cells by cell type (FAP, MuSC and Macrophage). The calculated interaction scores were classified under 4 categories in aged samples relative to young samples: UP-regulated for log-fold changes above log2 (1.5) and adjusted p-values below 0.05, DOWN-regulated for log-fold changes below −log2 (1.5) and adjusted p-values below 0.05, FLAT for absolute log-fold changes below log2 (1.5), and NSC (not significant) for absolute log-fold changes above log2 (1.5) with adjusted p-values above 0.05. We then calculated the over-representation scores of the Emitters, Receivers and Emitter-Receiver pairs and plotted them as a network.
Transcription factor footprinting
ATAC-seq read biological replicates were merged per condition (young and aged). Peaks per condition were called with MACS2 (–nomodel,–shift −100, and–extsize 200) and merged88,89. Tn5 bias, prediction of TF binding sites, and visualization of differentially bound TFs was done with TOBIAS90. Repetitive and ENCODE’s blacklisted regions were excluded from the analysis91 (Smit, A., Hubley, R. & Green, P. RepeatMasker. http://www.repeatmasker.org/ (2013). JASPAR motifs (non-redundant, vertebrate) were used as reference92.
Animal care
Mice were housed at a temperature of 21° with 20% humidity, with a light/dark cycle of 14 h light/10 h dark. All animal protocols and procedures carried out were approved by the McGill University Animal Care Committee (UACC).
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.