Multimodal cell atlas of the ageing human skeletal muscle

Muscle biopsy and ethical clearance

Samples were taken during orthopaedic surgery with informed consent from the 18 patients in the European cohort and the 13 patients in the Asian cohort; for one individual below 18 years, the informed consent was obtained from the legally acceptable representative. The study was performed following the Declaration of Helsinki. Ethical approval was granted for the European cohort by the Research Ethics Committee of Hospital Arnau de Vilanova (CEIm 28/2019), and for the Asian cohort by the Institutional Ethics Committee of the First Affiliated Hospital/School of Clinical Medicine of Guangdong Pharmaceutical University, Guangzhou (China) (2020-ICE-90). Exclusion criteria were myopathy, haemiplegia or haemiparesis, rheumatoid arthritis or other autoimmune connective tissue disorders, inability to consent, prior hospital admission in the previous month or major surgery in the previous 3 months. For the European cohort, the individuals’ medical and functional states were assessed according to the Barthel index¹⁰ and Charlson Index¹¹. The Barthel index estimates the grade of dependency of the individual ranging from 0 (totally dependent) to 100 (independent). The Charlson Index indicates the grade of comorbidities associated with the individual and ranged from 0 (without comorbidity) to 6 (the individual with a higher number of comorbidities) in our samples. A list of detailed information for the individuals is provided in Supplementary Table 1.

Animal experiment

C57Bl/6 (wild type) mice were bred and raised until 8–12 weeks of age at the animal facility of the Barcelona Biomedical Research Park (PRBB). They were housed in standard cages under a 12 h–12 h light–dark cycle and given unrestricted access to a standard chow diet. All experiments adhered to the ‘three Rs’ principle—replacement, reduction and refinement—outlined in Directive 63/2010 and its implementation in Member States. Procedures were approved by the PRBB Animal Research Ethics Committee (PRBB-CEEA) and the local government (Generalitat de Catalunya), following European Directive 2010/63/EU and Spanish regulations RD 53/2013. Both male and female mice were used for experiments and were maintained according to the Jackson Laboratory guidelines and protocols. Mice were randomly allocated to experimental or treatment groups. No blinding was used. No statistical methods were used to predetermine the sample size. Muscle injury was induced by intramuscular injection of CTX (Latoxan, L8102, 10 µM) and mice were euthanized at 7 days after injury as previously described¹².

Muscle sample processing

Muscle samples were obtained in all cases by selecting a macroscopically healthy area of muscle, without signs of contusion or haematoma. A small portion of muscle was removed by blunt dissection following the course of the myofibres and avoiding the use of electrocautery. The samples were immediately processed into three groups and stored next to the operating room as follows: (1) fixed with paraformaldehyde before being mounted in OCT compound as described previously (for immunochemistry and immunofluorescence)⁶⁰; (2) immediately frozen in liquid nitrogen (for snRNA-seq and snATAC-seq); and (3) tissue-digested (for scRNA-seq).

Single-cell preparation from skeletal muscle

Before the experiment, the post-operative muscle was immediately transferred in prechilled Dulbecco’s modified Eagle’s medium (DMEM, Corning, 10-017-CVR). For single-cell isolation, adipose and tendon tissues were removed using forceps, the remained muscle chunks were mechanically shredded on ice in a 10 cm plate. Next, prechilled DMEM medium was added to the plate for collecting muscle tissues and transferred into a 50 ml tube. After standing for 3 min, the supernatant containing the remaining adipose tissues was discarded. The remained muscle tissues were transferred to a 15 ml tube for digestion in 5 ml tissue digestion buffer (0.2 mg ml⁻¹ liberase (Roche, 5401119001), 0.4 μM CaCl₂ (Thermo Fisher Scientific, J63122AE), 5 μM MgCl₂ (Thermo Fisher Scientific, R0971), 0.2% BSA (Genview, FA016), 0.025% trypsin-EDTA (Thermo Fisher Scientific, 25300120). The muscles were digested in a shaking metal bath at 1000 rpm, 37 °C for 1 h, and mixed by inversion every 10 min. After all tissue pieces were digested, 3 ml of fetal bovine serum (FBS, Cellcook, CM1002L) was added to the mixture to terminate the digestion. The cell suspension was filtered through a 100 μm strainer, and centrifuged at 700g for 10 min at 4 °C to pellet the cells. The cell pellet was then resuspended in 10 ml wash buffer (DMEM medium supplemented with 10% FBS) and filtered through a 40 μm strainer, then centrifuged at 700g for 10 min at 4 °C to pellet the cells. The resultant single-cell suspensions were washed twice with prechilled PBS supplemented with 0.04% BSA and were used as input for scRNA-seq library construction.

Single-nucleus extraction from skeletal muscle

Single-nucleus isolation was performed as previously described⁶. In brief, tissues were thawed, minced and transferred to a 2 ml Dounce homogenizer (Sigma-Aldrich, D8938) with 1 ml of homogenization buffer A containing 250 mM sucrose (Sigma-Aldrich, S8501), 10 mg ml^–1 BSA, 5 mM MgCl₂, 0.12 U μl^–1 RNasin (Promega, N2115) and 1× cOmplete Protease Inhibitor Cocktail (Roche, 11697498001). Frozen tissues were kept in an ice box and homogenized by 25–50 strokes of the loose pestle (pestle A), after which the mixture was filtered using a 100 µm cell strainer into a 1.5 ml tube. The mixture was then transferred to a clean 1 ml Dounce homogenizer with 750 μl of buffer A containing 1% Igepal (Sigma-Aldrich, CA630), and the tissue was further homogenized by 25 strokes of the tight pestle (pestle B). The mixture was then filtered through a 40 μm strainer into a 1.5 ml tube and centrifuged at 500g for 5 min at 4 °C to pellet the nuclei. The pellet was resuspended in 1 ml of buffer B containing 320 mM sucrose, 10 mg ml⁻¹ BSA, 3 mM CaCl₂, 2 mM magnesium acetate, 0.1 mM EDTA (Thermo Fisher Scientific, 15575020), 10 mM Tris-HCl (Invitrogen, AM9856), 1 mM DTT (Invitrogen, 707265ML), 1× Complete Protease Inhibitor Cocktail and 0.12 U μl⁻¹ RNasin. This was followed by centrifugation at 500g for 5 min at 4 °C to pellet the nuclei. The nuclei were then washed twice with prechilled PBS supplemented with 0.04% BSA and finally resuspended in PBS at a concentration of 1,000 nuclei per μl for library preparation.

Library preparation and sequencing

sc/snRNA-seq library preparation

scRNA-seq libraries were prepared using the DNBelab C Series Single-Cell Library Prep Set (MGI, 1000021082)⁴⁹. In brief, the single-cell/nucleus suspensions were converted to barcoded scRNA-seq libraries through droplet encapsulation, emulsion breakage, mRNA-captured bead collection, reverse transcription, cDNA amplification and purification. Indexed sequencing libraries were constructed according to the manufacturer’s instructions. Library concentrations were quantified using the Qubit ssDNA Assay Kit (Thermo Fisher Scientific, Q10212). Libraries were sequenced using the DIPSEQ T1 sequencer.

snATAC-seq library preparation

snATAC-seq libraries were prepared using the DNBelab C Series Single-Cell ATAC Library Prep Set (MGI, 1000021878)⁴⁹. In brief, nuclei were extracted from tissue using the same protocol described above. After Tn5 tagmentation, transposed single-nucleus suspensions were converted to barcoded snATAC-seq libraries through droplet encapsulation, pre-amplification, emulsion breakage, captured bead collection, DNA amplification and purification. Indexed libraries were prepared according to the manufacturer’s instructions. Concentrations were measured with a Qubit ssDNA Assay Kit. Libraries were sequenced by a BGISEQ-2000 sequencer.

sc/snRNA-seq raw data processing, clustering and cell type annotation

Raw data processing

Raw sequencing reads were filtered, demultiplexed, and aligned to hg38 reference genome using a custom workflow (https://github.com/MGI-tech-bioinformatics/DNBelab_C_Series_HT_scRNA-analysis-software)⁴⁹. For scRNA-seq, reads aligned to gene exons were counted. For snRNA-seq, reads aligned to gene loci, including both exons and introns, were counted. Doublets were identified and filtered by DoubletFinder (v.2.0.3)⁶¹. Ambient RNA for snRNA-seq was reduced using SoupX (v.1.4.8)⁶² with the default settings.

Integration, clustering and cell type annotation

The resulting count matrix for cells/nuclei was filtered by the number of unique molecular identifiers (UMIs) > 1,000, gene > 500 and mitochondria content < 5%. Global clustering was performed using Scanpy (v.1.8.1)⁶³ in Python (v.3.7). Filtered data were normalized to total counts and log-transformed. The top 3,000 highly variable genes were selected, and the number of UMIs and the percentage of mitochondrial genes were regressed out. Each gene was scaled with the default options, followed by dimensionality reduction using principal component analysis. Batch effects between snRNA-seq and scRNA-seq were corrected using Harmony⁶⁴. Next, the batch-effect-corrected top 30 principal components were used for generating the neighbourhood graph with the number of neighbours set at 10. The cell clustering was further performed with the Louvain algorithm and annotated by canonical markers, putative scRNA-seq- and snRNA-seq-derived myofibre fragments were removed from the analysis. For satellite cell, immune cell, vascular cell and stromal cell reclustering, cells/nuclei were subset from the global clustering object and processed according to the same procedure as described above. For the reclustering of myonuclei, data were processed in Seurat (v.4.0.2)⁶⁵, and only snRNA-seq data were retained for further analysis. In brief, myonuclei data were subjected to SCTransform-based normalization, anchor identification between samples, integration, Louvain clustering and projection onto the UMAP space. Clustering results were further annotated by highly expressed genes.

Analysis of cell type composition variation in ageing

A generalized linear mixed model with a Poisson outcome¹⁴ was used to model the effect of age on cell-type-specific counts as previously reported, accounting for the possible biological (sex, ethnicity) and technical (omics, sequencing batch) covariates. The effect of each biological/technical factor on cell type composition was estimated by the interaction term with the cell type. The fold change is relative to the grand mean and adjusted. The statistical significance of the fold change estimation was measured by the LTSR, which is the probability that the estimated direction of the effect is true. As an alternative method, the proportion for each population was estimated over the total number of nuclei/cells for a given dataset (Supplementary Table 9).

Transcriptional and epigenetic heterogeneity analysis

Transcriptional heterogeneity analysis was performed as previously described⁶⁶. In brief, snRNA-seq data for each cell type in each age group were downsampled to 300 nuclei. For cell types with fewer than 300 nuclei, all nuclei were included for analysis. The resultant gene × cell matrix was further downsampled to make an equal number of UMI counts and cells between adult/older adult groups in each cell type. Next, all genes were ranked into ten blocks on the basis of the average expression value, and the 10% genes with the lowest coefficient of variation in each block were used to calculate the Euclidean distance between each cell. This Euclidean distance was used to measure transcriptional heterogeneity for each cell. For epigenetic heterogeneity, we adapted the same analysis method as transcriptional noise but using the rounded gene score matrix as input.

Myonucleus classification

Myonuclei were classified on the basis of previous markers associated with the described pure myofibre types (type I, type IIA and type IIX) and the hybrid myofibres (hybrid I/IIA, hybrid IIA/IIX). A module score was calculated for each myofibre type based on the expression of the following markers⁶⁷: type I (TNNT1, MYH7, MYH7B, TNNC1, TNNI1, and ATP2A2), type II (TNNT3, MYH1, MYH2, TNNC2, TNNI2, ATP2A1), type IIA (MYH2, ANKRD2, NDUFA8, MYOM3, CASQ2, HSPB6, RDH11, AIMP1) and type IIX (MYH1, MYLK2, ACTN3, MYBPC2, PCYOX1, CAPZA1, CD38, PDLIM7, COBL, TMEM159, HNRNPA1, TFRC). On the basis of these scores, myonuclei were first classified as type I, type II or hybrid I/IIa; thereafter, type II myonuclei were further classified as type IIA, type IIX or hybrid IIA/IIX. A residual amount of myonuclei remained unclassified due to the lower expression of these genes.

Differential gene expression and functional enrichment analysis

Seurat was used to compute the DEGs for each population and subpopulations between samples in the younger and older cohorts with the thresholds set at log₂[fold change] > 0.25 and Q < 0.05 (Supplementary Table 10). For myofibre subpopulations, the thresholds were set at: log₂[fold change] > 1 and Q < 0.05. The obtained DEGs for each comparison were used as input in Metascape online tool⁶⁸ to perform functional enrichment analysis, with a Q value threshold set at 0.05 (Supplementary Table 11). Heatmap results were plotted using pheatmap (v.1.0.12) in R.

Identification of coexpressing gene modules

Hotspot (v.1.1.1)¹⁶ was used to compute coexpressing gene modules among myofibre populations. The normalized expression matrix for the top 5,000 variable genes, the RegMyon-¹⁹, sarcomeric-⁶⁷ and atrophy-related²⁴ genes (Supplementary Table 3) were used as input. In brief, the k-nearest neighbour graph was created using the create_knn_graph function with the parameters: n_neighbors = 30, and then genes with significant correlation (Q < 0.05) were retained for further analysis. The modules were identified using the create_modules function with the parameters min_gene_threshold = 10 and fdr_threshold = 0.05.

Pseudotime analysis

For the myofibre degeneration trajectory, DCLK1⁺ (type I), ID1⁺ (type I), ID1⁺ (type II), ENOX1⁺ (type II) and other unperturbed myonuclei were selected for pseudotime analysis using Monocle3⁶⁹. After trajectory construction, myonuclei were ordered by pseudotime, and the corresponding gene expression matrixes were aggregated into 100 bins. The top 4,000 variable genes in type I or type II myonucleus trajectory were selected and visualized by k-means clustering heat map ordered by the pseudotime.

Cell–cell interaction analysis

CellChat (v.1.1.0)⁴⁵ detected ligand–receptor interactions on integrated sc/snRNA-seq data according to the standard procedures. The expression matrix and the cell type information were imported to CellChat. Specialized myonuclei, mast cells and erythrocyte clusters were removed from the analyses due to the insufficient number of cells/nuclei or the disproportionate number of cells/nuclei between the younger and older cohorts. The overall communication probability among the cell clusters was calculated using the computeCommunProb function with a trim set at 0.1.

snATAC-seq data processing

Raw data processing, clustering and cell type annotation

Raw sequencing reads were filtered, demultiplexed and aligned to the hg38 reference genome using PISA (https://github.com/shiquan/PISA)⁷⁰. Fragment files for each library were generated for downstream analysis. The transcription start site enrichment score, number of fragments and doublet score for each nucleus were calculated using ArchR⁷¹. Nuclei with transcription start site enrichment scores < 8 and number of fragments < 1,000 were removed from the analysis. Doublets were filtered out using the filterDoublets function with the settings filterRatio = 2. We next performed latent-semantic-indexing-based dimensionality reduction on the 500 bp tiles across the genome using the addIterativeLSI function of ArchR. Anchors between the scATAC-seq and scRNA-seq/snRNA-seq datasets were identified and used to transfer cell type labels identified from the scRNA-seq/snRNA-seq data. For co-embedding of snRNA-seq/scRNA-seq and snATAC-seq data, an anchor-based integration approach was applied based on the sequencing techniques. Then, data were further subjected to batch correction by Harmony among samples. Pearson’s correlation between snRNA/scRNA-seq and snATAC-seq was performed based on the integrated assay.

Motif enrichment analysis

Before motif enrichment, a reproducible peak set was created in ArchR⁷¹ using the addReproduciblePeakSet function based on cell types/subtypes. Differentially enriched peaks were identified using the getMarkerFeatures function with the thresholds set at log₂[fold change] > 0.5 and Q < 0.1. The motif presence in the peak set was determined with the addMotifAnnotations function using CisBP motif database (v.2)⁷².

TF occupancy

TF occupancy was evaluated by footprinting analysis implemented in ArchR⁷¹. In brief, putative binding sites of selectively enriched motifs were first inferred using the addMotifAnnotations function. Next, footprintings for the putative TF-binding sites were calculated using the getFootprints function, in which the Tn5 insertion bias was taken into account. The results were further plotted using the plotFootprints function.

GRN analysis

Construction of the GRNs was performed using FigR³¹. In brief, we first sampled an equal number of nuclei (20,000) in snRNA-seq and snATAC-seq analysis of myofibre and performed data integration using scOptMatch implemented in FigR. For creating the co-embedding map in these two independent datasets, we first input the variable features taken from the snRNA-seq and snATAC-seq datasets to perform CCA using the RunCCA function in Seurat. After integration, pairs of ATAC–RNA cells were identified by geodesic distance-based pairing using the pairCells function, and unpaired cells were removed from the analysis. Significant (P < 0.05) peak-to-gene associations were then identified among the cell pairs in type I or type II myonuclei. The DORCs were defined as peak-gene associations ≥ 6. For inference of the GRNs, the smoothed DORC score, RNA counts, snATAC-seq peak counts and the significant peak-to-gene associations were fed into runFigRGRN function, generating the GRNs. Next, the activators and repressors were identified by ranking the TFs by average regulation score.

GWAS analysis

Association of GWAS traits with skeletal muscle cell types

To identify trait/disease-relevant cell types, we performed LDSR analysis⁷³, a method for partitioning heritability from GWAS summary statistics. In brief, differentially accessible peaks for each adult/older adult cell type were identified (log₂[fold change] > 1 and Q < 0.01). The LDSC analysis was performed according to the standard workflow (https://github.com/bulik/ldsc/wiki). The summary statistics file for each trait was downloaded from the GWAS catalogue database⁷⁴ or published studies^50,51 (Supplementary Table 8a).

Fine mapping of non-coding variants and predicting the effect of TF binding

Lead SNPs were taken from low-hand-grip strength and lean-body-mass traits^50,51,75. FUMA, a web-based platform for GWAS analysis⁷⁶, was used to identify high-correlation SNPs with an LD r² ≥ 0.8 with lead SNPs. High-correlation SNPs within ±50 bp of the differentially accessible peaks were identified for further analysis. The peak-to-gene associations were determined using addPeak2GeneLinks function in ArchR package in the integrative object. To identify SNPs that affect TF binding, we used two approaches, (1) gkm-SVM⁵⁴ and (2) SNP2TFBS⁷⁷. For gkm-SVM, TF models were used from https://github.com/ren-lab/deltaSVM/tree/master/gkmsvm_models, and effective alleles were identified using the gkmExplain function⁷⁸. For SNP2TFBS tools, the analysis was performed in the SNP2TFBS web interface (https://ccg.epfl.ch/snp2tfbs/) following the tutorial.

Histology and immunofluorescence

Cryostat sections (10 μm thickness) were collected from muscles and stained with haematoxylin and eosin (Sigma-Aldrich, HHS80 and 45235) to assess tissue morphology or SA-β-gal (AppliChem, A1007,0001) for senescence cells with a modified staining protocol as described previously^12,79. Histochemical SDH staining was assayed by placing the slides in a solution containing sodium succinate (Sigma-Aldrich, S2378) as a substrate and nitro-blue tetrazolium (Sigma-Aldrich, N6876) for visualization of the reaction for 1 h at 37 °C. The intensity and pattern of staining were evaluated using light microscopy⁸⁰. Muscle collagen content was quantified after Sirius Red (Sigma-Aldrich, 365548) staining as previously described⁸¹. For immunofluorescence, the sections were air-dried, fixed, washed on PBS, permeabilized with Triton X-100 0.5% (Sigma-Aldrich, 11332481001) and incubated with primary antibodies (diluted as indicated below) after blocking with a high-protein-containing solution (BSA at 5%) (Sigma-Aldrich, A7906-100G) in PBS overnight at 4 °C. Subsequently, the slides were washed with PBS and incubated for 1 h at room temperature with the appropriate secondary antibodies diluted at 1:500; DAPI (Thermo Fisher Scientific, 62248) at 1:1,000 for nuclei; and WGA (Thermo Fisher Scientific, W11261) at 1:200 for cell/myofibre membrane. After washing, the tissue sections were mounted with Mowiol (Sigma-Aldrich, 81381) or Fluoromount-G (SoutherBiotech, 0100-01). Quantitative results for histology and immunofluorescence are listed in Supplementary Table 12. Primary antibodies were as follows: PAX7 (DSHB, PAX7, 1:50), PDGFRa (eBioscience, 17-1401-81, 1:100), perilipin-1 (Cell Signalling, 9349, 1:100), filamin C (MyBiosource, MBS2026155, 1:100), TNNT2 (Bioss, BS-10648R, 1:100), CD11b (eBioscience, 14-0112-85, 1:100), CD3 (Invitrogen, 14-0038-82, 1:100), CD19 (eBioscience, 14-0199-82, 1:100), NCAM1 (Cell Sciences, MON9006-1, 1:100), MYH7 (MyHC type I) (DSHB, A4.840, 1:10), MyHC type IIA/IIX (DSHB, SC-71, 1:70), laminin-647 (Novus Biologicals, NB300-144AF647, 1:200), FOS (Cell Signalling, 2250S, 1:200), ACVR2A (R&D, AF340, 1:100), ITGA7 (BioCell Scientific, 10007, 1:100), dystrophin (Sigma-Aldrich, D8168, 1:100). Secondary antibodies were as follows: goat anti-mouse IgM (DyLight 550, Invitrogen, SA5-10151), goat anti-mouse IgG1 (Alexa Fluor 488, Invitrogen, A21121), goat anti-mouse IgG (Alexa Fluor 488, Invitrogen, A11001), goat anti-mouse IgG (Alexa Fluor 568, Invitrogen, A11004), goat anti-rabbit IgG (Alexa Fluor Plus 488, Invitrogen, A32731TR), goat anti-rabbit IgG (Alexa Fluor Plus 647, Invitrogen, A32733TR), donkey anti-goat IgG (Alexa Fluor Plus 647, Invitrogen, A32849TR), goat anti-rat IgG (Alexa Fluor 568, Invitrogen, A11077).

Digital image acquisition and processing

Immunohistochemistry images were acquired using an upright microscope (Leica DMR6000B) equipped with a DFC300FX camera, and, for immunofluorescence pictures, using a Hamamatsu ORCA-ER camera. Images were acquired using HCX PL Fluotar objectives (×10/0.30 NA, ×20/0.50 NA and ×40/0.75 NA) and LAS AF software (Leica, v.4.0). Immunofluorescence pictures were also obtained using the Nikon Ti2 fluorescence microscope with NIS Elements software (v.4.11.0), and a confocal microscope (Zeiss 980 Airyscan2) with ZenBlue software (v.3.5) and a ×20 air objective. The acquired images were composed, edited and analysed using Fiji (ImageJ, v.2.14.0/1,54f). To reduce background, brightness and contrast adjustments were applied to the entire image. Myofibre size was assessed using the MyoSight tool³⁴, with a manual correction applied after automated outlining, and the cross-sectional area (CSA) was determined using Fiji. Signals of SA-β-gal, PAX7, PDGFRα, perilipin, CD11B, CD3, CD19, TNNT2, NCAM1, filamin C, SDH and FOS staining were manually counted in Fiji. The area of ACVR2A, Sirius Red and ITGA7 staining was calculated by normalizing the positive-signal area to the total imaged area in Fiji.

Statistical analysis

The sample size of each experimental group or number of independent experiments is described in the corresponding figure legend. The calculation method for P values is explained in the figure legends. The number of replicates for each experiment is presented in the figure legends. For Pearson’s correlation, statistical significance for positive or negative correlation (represented as the R value) was set at P < 0.05 and shading represents the 95% confidence interval along the correlation line (Supplementary Table 5). For the box plots, the central line shows the median, the box limits indicate the upper and lower quartiles, and the whiskers indicate 1.5× the interquartile range. Python, R or Prism (v.10) were used for statistical analyses.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.