Next Article in Journal
Antigen Cross-Presentation by Murine Proximal Tubular Epithelial Cells Induces Cytotoxic and Inflammatory CD8+ T Cells
Next Article in Special Issue
Update on the Molecular Aspects and Methods Underlying the Complex Architecture of FSHD
Previous Article in Journal
MicroRNA-145 Impairs Classical Non-Homologous End-Joining in Response to Ionizing Radiation-Induced DNA Double-Strand Breaks via Targeting DNA-PKcs
Previous Article in Special Issue
The Evolution of Complex Muscle Cell In Vitro Models to Study Pathomechanisms and Drug Development of Neuromuscular Disease
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Cells
Volume 11
Issue 9
10.3390/cells11091508
cells-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Unraveling the Molecular Basis of the Dystrophic Process in Limb-Girdle Muscular Dystrophy LGMD-R12 by Differential Gene Expression Profiles in Diseased and Healthy Muscles

by
Christophe E. Depuydt
1,
Veerle Goosens
2,
Rekin’s Janky
3,
Ann D’Hondt
4,
Jan L. De Bleecker
5,
Nathalie Noppe
2,
Stefaan Derveaux
3,
Dietmar R. Thal
6,7 and
Kristl G. Claeys
1,4,*
1
Laboratory for Muscle Diseases and Neuropathies, Department of Neurosciences, KU Leuven, and Leuven Brain Institute (LBI), Herestraat 49, 3000 Leuven, Belgium
2
Department of Radiology, University Hospitals Leuven, Herestraat 49, 3000 Leuven, Belgium
3
VIB Nucleomics Core, Herestraat 49, 3000 Leuven, Belgium
4
Department of Neurology, University Hospitals Leuven, Herestraat 49, 3000 Leuven, Belgium
5
Department of Neurology, University Hospital Gent, Corneel Heymanslaan 10, 9000 Gent, Belgium
6
Department of Pathology, University Hospitals Leuven, Herestraat 49, 3000 Leuven, Belgium
7
Laboratory for Neuropathology, Department of Imaging and Pathology, KU Leuven, and Leuven Brain Institute (LBI), Herestraat 49, 3000 Leuven, Belgium
*
Author to whom correspondence should be addressed.
Cells 2022, 11(9), 1508; https://doi.org/10.3390/cells11091508
Submission received: 27 March 2022 / Revised: 19 April 2022 / Accepted: 29 April 2022 / Published: 30 April 2022
(This article belongs to the Special Issue Cellular and Molecular Mechanisms in Neuromuscular Diseases: From Basic Research to Clinical Approach)
Browse Figures
Review Reports Versions Notes

Abstract

:
Limb-girdle muscular dystrophy R12 (LGMD-R12) is caused by two mutations in anoctamin-5 (ANO5). Our aim was to identify genes and pathways that underlie LGMD-R12 and explain differences in the molecular predisposition and susceptibility between three thigh muscles that are severely (semimembranosus), moderately (vastus lateralis) or mildly (rectus femoris) affected in this disease. We performed transcriptomics on these three muscles in 16 male LGMD-R12 patients and 15 age-matched male controls. Our results showed that LGMD-R12 dystrophic muscle is associated with the expression of genes indicative of fibroblast and adipocyte replacement, such as fibroadipogenic progenitors and immune cell infiltration, while muscle protein synthesis and metabolism were downregulated. Muscle degeneration was associated with an increase in genes involved in muscle injury and inflammation, and muscle repair/regeneration. Baseline differences between muscles in healthy individuals indicated that muscles that are the most affected by LGMD-R12 have the lowest expression of transcription factor networks involved in muscle (re)generation and satellite stem cell activation. Instead, they show relative high levels of fetal/embryonic myosins, all together indicating that muscles differ in their baseline regenerative potential. To conclude, we profiled the gene expression landscape in LGMD-R12, identified baseline differences in expression levels between differently affected muscles and characterized disease-associated changes.

1. Introduction

The limb-girdle muscular dystrophies (LGMDs) constitute a group of rare, progressive and genetic muscle disorders, with weakness and atrophy of mainly pelvic and shoulder girdle muscles [1]. LGMDs are inherited mostly in an autosomal recessive manner (LGMD-R), more frequent than autosomal dominant (LGMD-D). Autosomal recessive pathogenic variants in the anoctamin-5 encoding gene (ANO5) cause LGMD-R12 anoctamin5 related (LGMD-R12) and distal Miyoshi muscular dystrophy type 3 (MMD3) [2,3,4]. In ANO5-related muscular dystrophies, a male predominance and women often showing a less severe phenotype have been reported [4]. Dominant ANO5 mutations cause the bone disorder gnathodiaphyseal dysplasia 1 (GDD1) [5,6]. To date, pathophysiology of ANO5-related muscular dystrophies is largely unknown and disease-specific treatments do not exist.
The ANO5-gene is highly expressed in skeletal and cardiac muscles and in bones [6]. ANO5 or transmembrane 16E protein (ANO5/TMEM16E) belongs to a family of ten transmembrane proteins that have a role either as calcium-activated ion channels, lipid scramblases, or both [7,8]. ANO5 is the only member of this protein family associated with muscular dystrophy. Its lipid scramblase and ion channel activities have been shown to play a role in sarcolemmal repair and myoblast fusion during muscle regeneration [9,10,11,12,13].
Similar to other muscular dystrophies, the muscle biopsy in LGMD-R12 reveals dystrophic changes consisting of necrotic and regenerating muscle fibers, and replacement of muscle cells by connective and fatty tissue as disease progresses [3,14]. In addition, muscle biopsy in LGMD-R12 often reveals inflammatory infiltrates [15,16,17]. Although pathomechanisms of the dystrophic process in LGMD-R12 and other muscular dystrophies are largely unknown, several pathways that might contribute to progressive muscle cell death have been suggested, such as oxidative stress [18], mitochondrial dysfunction [19], defective membrane repair [9], inflammation [20], compromised regeneration [13] and impaired satellite cell activation [21].
Interestingly, although the genetic defect is the same in all muscles in LGMD-R12, specific muscles are selectively involved, which can be seen on muscle Magnetic Resonance Imaging (MRI). LGMD-R12 patients show a predominant affection of the posterior thigh muscles, starting at the semimembranosus, whereas vastus lateralis is affected later in the disease course, and rectus femoris, gracilis, sartorius and biceps femoris short head show fatty infiltration only at more advanced stages of the disease [4,22,23]. A selective muscle involvement also occurs in other muscular dystrophies, but the patterns of affected muscles differ depending on the mutated gene [23]. It is highly relevant from a therapeutic perspective to understand the reasons why some muscles escape the dystrophic effects of disease-causing mutations.
The degree of muscle involvement can be determined semiquantitatively on muscle MRI using various scales, such as the 4-point Mercuri score, which has been shown to correlate with histopathological findings and disease progress [24,25].
The objectives of this study were to identify genes and pathways that underlie LGMD-R12, resulting in new molecular insights into the dystrophic process, and explain differences in the molecular predisposition and susceptibility between different muscles. We used RNA-seq and studied the gene expression profiles in three differentially affected muscles of the thigh (severely affected semimembranosus, intermediately affected vastus lateralis, mildly affected to preserved rectus femoris) in 16 male LGMD-R12 patients and in the same three muscles in 15 healthy male age-matched control individuals. We identified marker genes and pathways for the different muscles in patients and in controls and applied single cell integration and deconvolution analysis to estimate how disease affects different cell types that compose human skeletal muscle tissue.

2. Patients and Methods

2.1. Patients and Controls

We included 16 symptomatic, ambulatory, adult, male patients with genetically proven LGMD-R12, i.e., carrying two pathogenic variants in the anoctamin-5 gene (ANO5) (Table 1), and 15 age-matched healthy male control individuals. Patients were followed at the Neurology Department, Neuromuscular Reference Center, at the University Hospitals Leuven and Gent, Belgium. Wheelchair-bound patients and/or patients with complete fatty infiltration of the three target muscles (semimembranosus, vastus lateralis, rectus femoris) on muscle MRI were excluded, as well as patients with a predominantly distal Miyoshi muscular dystrophy type 3 (MMD3) phenotype or with asymptomatic hyperCKemia. Further exclusion criteria were: presence of a contra-indication for MRI (e.g., MRI non-compatible pacemaker or claustrophobia), a blood coagulation disorder, use of anticoagulation medication or potentially toxic medication on muscles such as steroids or statins at the time of muscle biopsy and active alcohol abuse.
In all participants, the 6-min walk distance (6MWD) and 10-m walk test (10MWT) were performed (Table 1) and a whole-body muscle MRI (1.5 Tesla, Philips Ingenia, Philips Medical Systems, Best, The Netherlands) with axial and coronal T1-weighted scans was carried out, prior to muscle biopsy sampling. The target muscles were scored on the MR images using the Mercuri score, which is a 4-point grading system to categorize disease severity in individual muscles based on visual inspection of fatty tissue infiltration on MRI: normal appearance (score 1/normal), less than 30% affected (score 2/mild), between 30–60% affected (score 3/moderate), between 60–100% affected (score 4/severe) [24].
The study was approved by the Ethics Committee Research UZ/KU Leuven (S-59867). Written informed consent was obtained from all participants.

2.2. Muscle Biopsies

In all patients and controls, three muscle biopsies were taken, one in the semimembranosus (SM), one in the vastus lateralis (VL) and one in the rectus femoris (RF) muscle. We performed vacuum-assisted needle biopsies using the EnCor Enspire (breast) biopsy system with 10G EnCor needles (Bard Benelux, Olen, Belgium). For every muscle biopsy a separate needle was used. Muscle biopsies were taken under ultrasound guidance after injecting local anesthesia (lidocaine 2.0%) in the skin, subcutaneous fat and fascia but not in the muscle itself. Prior to the muscle biopsy procedure, a whole-body muscle MRI including the thigh muscles was carried out to determine the level for biopsy sampling, in order to avoid completely fatty replaced tissue. Since the muscle transcriptome is known to be altered by several variables, the following measures were taken, as well as matching for age and sex, to eliminate possible confounding factors: biopsies were taken in the morning, to avoid differences in circadian rhythm between participants [26], and participants were instructed not to perform intensive physical activities for at least one week before the muscle sampling, in order to avoid bias from differences in exercise or work load [27,28].
Immediately after sampling, the biopsies from each of the three muscles were mounted on a separate cork and snap frozen in isopentane that was cooled with liquid nitrogen. All biopsies were stored at −80 °C. For RNA extraction, 200 slices of 5 µm thickness from each biopsy were cut with a microtome-cryostat and collected in two RNase-free tubes of 2.0 mL (100 slices in each tube), which were kept frozen at −80 °C until RNA extraction was performed.

2.3. RNA Procedures

2.3.1. RNA Extraction

We used the same protocol for RNA extraction from muscle tissue as the one described by Cummings et al. [29]. RNA was extracted from the 93 muscle biopsies using the Direct-zol RNA Miniprep Plus kit ZY-R2073 (Zymo Research, Irvine, CA, USA) according to the manufacturer’s instructions. In order to avoid technical variability, all samples were extracted with the same kit on two consecutive days in the same lab by the same lab technician.

2.3.2. RNA Quality Control

The 93 extracted RNA samples were treated with DNase prior to their submission to the VIB Nucleomics Core (Leuven, Belgium, www.nucleomics.be, accessed on 15 July 2019). The Nanodrop ND-1000 (Nanodrop Technologies, Wilmington, DE, USA) was used to measure RNA concentration and purity spectrophotometrically. A Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA) was applied to assess RNA integrity. The RNA Integrity Number (RIN) scores mostly varied from 6.3 (partially degraded) to 8.4, except for a few samples that were highly degraded.

2.3.3. Library Preparation

Per sample, an amount of 5 ng of total RNA (or 7 µL RNA when 5 ng was not available) was used as input for the SMART-Seq Stranded Kit (Cat. No. 634444; protocol version “022819”; low input option; Takara Bio USA, San Jose, CA, USA). This kit can deal with degraded as well as high-integrity input RNA. Positive and negative controls were included in the experimental design using 5 ng control RNA (included in the kit) and 7 µL RNase-free water, respectively.
First, RNA is converted to cDNA using random priming (scN6 Primer) and SMART (Switching Mechanism At 5′ end of RNA Template) technology (Takara Bio USA, San Jose, CA, USA) and then full-length adapters for Illumina sequencing (including specific barcodes for dual-indexing libraries) are added through PCR using a limited number of cycles (5 cycles). The PCR products are purified and then ribosomal cDNA is selectively depleted by cleaving the ribosomal cDNAs by scZapR in the presence of mammalian-specific scR-Probes, which target nuclear and mitochondrial rRNA sequences. The library fragments derived from non-rRNA molecules remain untouched by this process. The remaining cDNA fragments are further amplified with primers universal to all libraries (14 cycles). Lastly, the PCR products are purified once more to yield the final cDNA library. All libraries were finally quantified using Qubit dsDNA HS kit (Thermo Fisher Scientific, Waltham, MA, USA) and their size distribution was checked using a Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA). Nine of the 93 libraries were excluded from downstream sequencing as there was no yield: control C6 (VL), control C8 (SM), patient 6 (SM), patient 7 (SM), patient 8 (VL), patient 12 (RF) and patient 16 (SM, VL, RF) (Supplementary Tables S1 and S2).

2.3.4. Sequencing (RNA-Seq)

Sequence-libraries of each of the 84 remaining samples were equimolarly pooled and sequenced on two Illumina NovaSeq 6000 S1 100 flow-cells (Xp workflow, Paired Read 51-8-8-51, 1% PhiX v3) (Illumina, San Diego, CA, USA) at the VIB Nucleomics Core (Leuven, Belgium, www.nucleomics.be, accessed on 12 September 2019). RNA-seq was performed twice, in two runs. In the first run, data from 84 samples were analyzed (i.e., 93 samples minus 9 excluded samples because of too low RNA amount). The first run was successful with high quality (~95% of the bases ≥30) and output according to the expectations (~2000 M Passed Filter clusters). With 25M PF reads on average per sample, the output per sample was quite variable, ranging from 9.4M to 47.4M PF reads. We also observed a high percentage of duplicates in some samples and a high percentage of adapters in the sequences, corresponding to samples with very low input material. The mapping of the reads is relatively correct (90%), but the filtering based on the mapping quality has a huge impact (~50% loss of reads for many samples), which should be due to duplicated reads as the filtering discards similar reads. Out of 84 samples, we selected the 38 best samples (<22% No Features, Number of final counts > 12M), whereas 46 samples were pooled for a second sequencing run in order to reach 20M final counts (=assigned reads) for each sample.

2.4. Data Analysis

2.4.1. Preprocessing

Using FastX 0.0.14 and Cutadapt 1.15 (http://hannonlab.cshl.edu/fastx_toolkit/index.html, accessed on 20 September 2019), low quality ends and adapter sequences were trimmed off from the Illumina reads. Next, the following were filtered using FastX 0.0.14 and ShortRead 1.36.1 [30]: small reads (length < 35 bp), polyA-reads (more than 90% of the bases equal A), ambiguous reads (containing N), low-quality reads (more than 50% of the bases < Q25) and artifact reads (all but three bases in the read equal one base type). We then identified and removed reads that align to phix_illumina applying Bowtie2 2.3.3.1 [31].

2.4.2. Mapping

We aligned the preprocessed reads to the reference genome of Homo sapiens (GRCh38) with STAR aligner v2.5.2b [32]. Default STAR aligner parameter settings were used, except for “—outSAMprimaryFlag OneBestScore—twopassMode Basic—alignIntronMin 50–alignIntronMax 500,000—outSAMtype BAM SortedByCoordinate”. Reads with a mapping quality smaller than 20 were removed from the alignments using Samtools 1.5 [33].

2.4.3. Counting

We counted the number of reads in the alignments that overlap with gene features using featureCounts 1.5.3 [34]. We chose the following parameters: -Q 0 -s 2 -t exon -g gene_id, and removed genes for which all samples had less than 1 count per million. We further corrected raw counts within samples for GC-content and between samples using full quantile normalization, with the EDASeq package from Bioconductor [35].

2.5. RNA-Seq Analysis

Raw count matrices were uploaded to the online R-based UniApp data analysis platform (Unicle Biomedical Data Science, Leuven, Belgium) for analysis. We used the R packages plotly and ggplot to generate bar, pie, line and scatter plots.

2.5.1. Quality Control and Data Normalization

Genes expressed at a level of at least 1 count per million reads in at least 10% of samples were filtered and normalized using the EdgeR package. Samples with less than 1 Million reads were filtered out from the final analysis.

2.5.2. Principal Component Analysis (PCA)

We first normalized the raw count RNA-seq data from 84 patients for total read depth. Next, we performed quality filtering for read mapping and the number of detected genes and excluded 5 samples: patient 3 (SM), patient 5 (SM, VL), patient 12 (SM, VL) (Supplementary Table S1). The normalized data were auto-scaled and PCA was performed on the top 2000 most highly variable genes (with the Seurat R package [36]) to build a 2-dimensional representation of the data.

2.5.3. Pair-Wise Differential Analysis

We used limma [37] for differential expression analysis between two specific clusters and visualized using a volcano plot.

2.5.4. Gene Set Enrichment Analysis (GSEA)

We applied gene set enrichment analysis (clusterProfiler R package) to compare gene expression profiles between two groups [38]. We performed GSEA using a subset of gene sets selected from the Molecular Signatures Database (MSigDB version 7.41; http://bioinf.wehi.edu.au/software/MSigDB/, accessed on 23 September 2021), which is a collection of annotated gene sets. GSEA scores were calculated for sets with at least five detected genes, all other parameters were considered default.

2.5.5. Heatmap Analysis

Gene expression heatmaps are based on averaged auto-scaled data. We produced heatmaps by the heatmaply R package.

2.5.6. Deconvolution Analysis

Bulk RNA-seq count data were deconvoluted into cell composition matrices with the MUSIC algorithm [39] on a reference single cell RNA-seq dataset derived from the mononuclear cells of the human vastus lateralis [40] with default parameters.

2.6. Single Cell scRNA-Seq Analysis

We used a publicly available single cell dataset of human muscle to obtain markers for different stromal cell populations. The raw count matrix of human muscle cells as described in the study “Single cell transcriptional profiles in human skeletal muscle” was downloaded from the GEO repository under the accession code GSE130646 [40].

2.6.1. Quality Control and Data Normalization

We applied the following quality control steps for the human freshly-isolated muscle cells: (i) genes expressed by <3 cells were not considered; (ii) cells that had over 2000 expressed genes (possible doublets), or over 4% of unique molecular identifiers (UMIs) derived from the mitochondrial genome were deleted. The data of the remaining cells were natural-log transformed using log1 p and normalized using the Seurat package [36].

2.6.2. Dimension Reduction

After auto-scaling, genes with high variability were identified using the Seurat FindVariableGenes function. The top 2000 most highly variable genes were included in the analysis. Subsequently, the normalized data were first summarized by PCA, and the first 20 PCAs were visualized using t-Distributed Stochastic Neighbor Embedding (t-SNE, Rtsne package) with a perplexity value of 100 and a learning rate of 100. This representation was only used to visualize the data.

2.6.3. Clustering Analysis and Annotation

We performed graph-based clustering to cluster cells according to their gene expression profiles as implemented in Seurat with the first 20 PCA dimensions, number of neighbors set to 10 and at a resolution of 0.8. Cell clusters were annotated based on canonical markers.

2.6.4. Pair-Wise Differential Analysis

We performed differential expression analysis between two specific clusters using limma [37].

2.6.5. Marker Gene Analysis

A two-step approach was used to obtain ranked marker gene lists for each cluster. Marker genes for a given cluster should have the highest expression in that cluster compared to all other clusters and are therefore uniquely assigned to one cluster. Next, marker genes were ranked using a product-based meta-analysis [41]. We performed pair-wise differential analysis of all clusters against all other clusters separately and ranked the results of each pair-wise comparison by log2 fold change. The most downregulated genes received the highest rank number and the most upregulated genes received the lowest rank number (top ranking marker genes). For each cluster, we combined the rank numbers for all genes in all pair-wise comparisons by calculating their product to obtain a final list of ranked marker genes for each cluster. Clusters were annotated based on literature-curated marker genes of canonical muscle cells. Cells that could not be unambiguously assigned to a biologically meaningful phenotype were excluded from the analysis because they might represent low quality cells or doublets.

2.6.6. Heatmap Analysis

To account for cell-to-cell transcriptomic stochastics, all heatmaps are based on cluster-averaged gene expression. For visualization, data were auto-scaled. Using the heatmaply package, heatmaps were produced.

3. Results

3.1. Demographic Data

Mean age at inclusion was 45.7 years (range 26 to 64 years) for the patients (Table 1) and 45.1 years (range 25 to 63 years) for the controls (p > 0.05). Further clinical and genetic data of the patients are detailed in Table 1.

3.2. Mercuri Score Is Correlated with Functional Outcomes in LGMD-R12 Patients

We used the 4-point Mercuri score to categorize disease severity in the three target muscles (SM, VL, RF) based on visual inspection of fatty tissue infiltration on muscle MRI in all patients and controls [24]. Mercuri scores in patients’ muscles varied between 1 and 4 (Table 1, Figure 1A), whereas in controls the Mercuri score for all muscles was equal to 1 (i.e., normal) (Figure 1B). The Mercuri score in the patients was lowest for the rectus femoris with mean 1.2 (95% CI 1.0–1.4) (least affected), intermediate for the vastus lateralis with mean 2.0 (95% CI 1.3–2.7) (intermediately affected) and highest for the semimembranosus with mean 2.8 (95% CI 2.1–3.4) (most affected) (ANOVA, p = 0.001) (Figure 1C), which is consistent with the literature [4,22,23]. We calculated the Mercuri score with years after disease onset and showed that the Mercuri score significantly increased over time with increasing disease duration (Pearson correlation coefficient r = 0.5666, p = 0.0331) (Figure 1D). Two patients (Table 1, patients 7 and 11) showed a longer disease duration with, however, lower Mercuri scores (Figure 1D), which might be explained by the different genetic defects in ANO5 of these patients compared to the others (Table 1). Functionally, the Mercuri score significantly increased with decreasing 6MWD (r = −0.8387, p < 0.001; Figure 1E) and with increasing 10MWT (r = 0.8543, p < 0.001; Figure 1F), which corresponds to the progressive functional decline of the disease over time.

3.3. Unbiased Analysis: Gene Expression Profiles Correlate with the Mercuri Score

We performed principal component analysis (PCA) using highly variable genes as input and correlated gene expression signatures with the Mercuri score (Figure 2). Interestingly, color coding the unbiased PCA plot for the Mercuri score showed that the first component (the one that explains most variability in the data and hence is the strongest signature) was highly correlated with the Mercuri score. The second component was correlated with the mitochondrial read count (Figure 2). To determine whether the gene signature that correlated with the Mercuri score was based on genes involved in muscle fiber degradation and connective and fat tissue deposition, we color coded the PCA plot for expression of a known muscle membrane gene (alpha-sarcoglycan, SGCA) [42], fibroblast marker (collagen I A1, COL1A1) [43] and adipocyte marker (adiponectin, ADIPOQ) [44]. This analysis showed a gradient of decreasing expression of muscle marker genes and a gradual increase in the expression of fibroblast, adipocyte and other stromal cells (Figure 2). Thus, our exploratory PCA analysis was consistent with the replacement of muscle tissue by connective and fat tissue as disease progresses.

3.4. Quantification of Gene Signatures Associated with LGMD-R12

We next aimed to formally quantify changes in gene expression signatures, by performing a differential gene expression analysis of LGMD-R12 patients versus controls (irrespective of the Mercuri score, muscle or other stratification criteria) (Figure 3A). This analysis showed in patients an increase in genes involved in response to muscle injury and inflammation, such as osteopontin (SPP1) [45,46], immunoglobulin heavy constant gamma 1/2 (IGHG1/2), leukocyte receptor tyrosine kinase (LTK), joining chain of multimeric IgA and IgM (JCHAIN), serum amyloid A1 (SAA1), as well as embryonic and fetal myosins, including myosin heavy chain 3 (MYH3) and myosin light chain 4 (MYL4). Embryonic myosins are normally not expressed in adult muscle but can be upregulated in response to injury in regenerating muscle fibers [47]. Interestingly, the LGMD-R12 disease-causing gene ANO5 ranked in the top 0.1% of most downregulated genes, whereas other anoctamin genes showed a trend to being slightly upregulated (ANO3, ANO7, ANO9, ANO7L1) (Figure 3A). These findings correspond to recent hypotheses that anoctamin homologs might compensate for the loss of ANO5 in mice [48].
To quantify gene signature changes at the gene set level we performed a gene set enrichment analysis (GSEA) of the hallmark gene sets based on the results of the differential gene expression analysis of LGMD-R12 patients versus controls (Figure 3B). This analysis showed that LGMD-R12 patients have upregulated gene sets related to immune response, phagocytosis and collagen remodeling while having downregulated gene sets relating to energy metabolism and protein synthesis.

3.5. Pathway Analysis Characterizes Gene Sets for LGMD-R12 Progression

To characterize gene sets for LGMD-R12 progression we performed GSVA analysis to alter the gene-by-sample matrix to a gene set-by-sample matrix and used the latter matrix to perform a differential analysis between semimembranosus muscle samples with Mercuri score 1 versus Mercuri score 4 (Figure 3C,D) [49]. This calculation showed that the cell type signature of skeletal muscle cells was downregulated in Mercuri score 4 samples, while stromal and immune cell signatures were upregulated. Consistently, when evaluating cellular pathways in a similar analysis, we found that pathways associated with skeletal muscle were downregulated (ribosomes, mitochondria and aerobic respiration), while gene sets associated with cell death, inflammation and collagen metabolism were upregulated (Figure 3C,D). To evaluate gradual differences in pathway expression between different Mercuri scores, we performed marker pathway analysis on the GSVA transformed matrix and visualized the results as a heatmap. This assessment showed a progressive downregulation of pathways such as myogenesis, protein metabolism and oxidative phosphorylation, while pathways associated with inflammation and cell death were upregulated (Figure 3E).

3.6. Deconvolution Analysis Shows an Increase in Fibroadipogenic Progenitor (FAP) Cells in LGMD-R12 Muscles

Our analysis of pathways associated with LGMD-R12 were consistent with a decrease in muscle and an increase in stromal cell signatures. However, this analysis did not provide insight into whether these findings reflect an overall change in gene expression in all cells of the muscle, or if the actual cellular composition of the disease muscle is changed. To estimate differences in cell type composition we performed a deconvolution analysis [50]. For this analysis we used marker genes derived from a recently described skeletal muscle single cell RNA-seq dataset, cell types that together make up human skeletal muscle tissue (Supplementary Figure S1) [40], and deconvoluted our bulk RNA-seq dataset. While preliminary, our analysis indicated a progressive increase in FAP lumican+ (LUM+) cells and their signatures but a decrease in (vascular) smooth muscle signatures (vSMC) (Figure 3F,G).

3.7. Distinct Gene Signatures in Different Muscles in Healthy Controls

To investigate the hypothesis that baseline differences in the molecular composition of different muscles influence susceptibility to muscular dystrophy, we performed a detailed analysis of gene signatures of the three biopsied muscles in healthy control individuals. We first performed a PCA analysis, which showed a clear separation of semimembranosus, vastus lateralis and rectus femoris on the third component (this first component highlighted a lowly sequenced outlier sample) (Figure 4A). Color coding PCA plots for marker genes for fast- and slow-twitch muscle showed a clear gradient, with the rectus femoris showing the highest expression of the fast-twitch muscle myosins (myosin heavy chain 4 (MYH4), myosin heavy chain 2 (MYH2), myosin heavy chain 1 (MYH1) and myosin light chain 1 (MYL1)) and the semimembranosus presenting the highest expression of the slow-twitch markers (myosin heavy chain 7 (MYH7), myosin heavy chain 6 (MYH6) and myosin light chain 3 (MYL3)) (Figure 4B,C) [51,52]. The vastus lateralis expressed intermediate levels of fast- and slow-twitch myosins. Differential analysis and gene set analysis between the most affected (semimembranosus) and less affected (rectus femoris) muscles showed differences in genes that regulated anterior–posterior axis planning, fast- and slow-twitch muscle fibers and satellite cell marker genes (Figure 4D,E).

3.8. Differential Gene Expression between Different Muscles: Identification of Genes Associated with Different Muscles

To unbiasedly assess genes enriched in the different biopsied muscles (SM, VL, RF) we performed a marker gene analysis (Figure 4F–H; Table 2 and Table 3). In this analysis, each muscle is compared to the other muscles separately via a differential gene expression analysis. The results of the different analyses are then combined with a product-based meta-analysis based on the log-fold changes [53]. We used a heatmap to visualize the top 20 marker genes. Interestingly, several transcription factors were highly differentially expressed across semimembranosus, such as Heart And Neural Crest Derivatives Expressed 2 (HAND2), HAND2 Antisense RNA 1 (HAND2-AS1), Homeobox D8 (HOXD8) and Homeobox D9 (HOXD9), vastus lateralis, such as T-Box Transcription Factor 5 (TBX5), Homeobox A13 (HOXA13), T-Box Transcription Factor 5 Antisense RNA 1 (TBX5-AS1), and rectus femoris such as Transcription Factor 24 (TCF24) (Figure 4F). In addition, we identified differential expression of genes involved in the regulation of protein synthesis in muscle tissue (METTL21C) [54] and glycogen breakdown were highest in the rectus femoris, while adipogenic and lipid pathways were upregulated in the semimembranosus (Figure 4F–H). A similar analysis of early stage (Mercuri score 1 and 2) and late stage (Mercuri score 3 and 4) muscles showed that most variation in disease muscle was explained by disease processes (Supplementary Figure S2; Supplementary Table S3). Together, these findings may suggest that muscles differ in their baseline regenerative potential.

4. Discussion

The main objective of this study was to identify genes and pathways that underlie LGMD-R12 and explain differences in the molecular predisposition and susceptibility between three different thigh muscles that are severely (SM), moderately (VL), or mildly (RF) affected in this disease. We performed an unbiased analysis of 84 muscle biopsies (79 after quality checks) of these three differently affected muscles in 16 LGMD-R12 patients and 15 age-matched controls (all males). Our results indicated that, at the cellular and molecular level, LGMD-R12 muscular dystrophy is characterized by the expression of genes indicative of fibroblast and adipocyte replacement, such as fibroadipogenic progenitor (FAP) cells and immune cell infiltration, while gene signatures associated with striated muscle (protein synthesis, OXPHOS, glycogen-, glucose- and amino acid metabolism) are downregulated in dystrophic muscle. Furthermore, muscle degeneration as quantified by the radiological Mercuri score is associated with an increase in genes associated with muscle injury and inflammation, as well as genes involved in muscle repair/regeneration (including embryonic and fetal myosins). We also identified interesting differential expression patterns in other anoctamin genes that are upregulated in response to ANO5 loss. Analysis of baseline differences in between muscles in healthy individuals indicated that muscles that are the most affected by LGMD-R12 have the lowest expression of transcription factor networks involved in muscle (re)generation and satellite (stem) cell activation. Instead, they show relatively high levels of fetal/embryonic myosins. This is the first study on transcriptomics in LGMD-R12 and the first study overall in which three different muscles were biopsied in patients and in healthy controls, enabling us to conclude to these important and highly relevant findings.

4.1. Genes Involved in Inflammation Are Upregulated in Dystrophic LGMD-R12 Muscles

Our study showed an increase in genes involved in muscle injury and inflammation in dystrophic muscles of LGMD-R12 patients. The number of different immune cells is largely increased compared with healthy muscle. This has been reported in muscles in Duchenne muscular dystrophy (DMD) patients and mdx-mice, including T-lymphocytes (CD4, CD8 and regulatory T-cells or TRegs), natural killer (NK) cells, neutrophils, eosinophils and macrophages [20,55,56]. Furthermore, histological analyses of muscle biopsies from LGMD-R12 patients often show inflammatory changes, characterized by CD45 and CD8 positive leukocytes as well as CD68 and CD206 positive macrophages accumulating within myofibers showing myophagocytosis [16,17]. Inflammatory changes have been reported in the muscle biopsies of patients with other muscular dystrophies as well, such as DMD [57], facioscapulohumeral dystrophy (FSHD) [58], LGMD-R2 [59], LGMD R3 and LGMD R5 [60]. Moreover, on muscle MRI, there are indications that inflammation plays a role in muscular dystrophies, such as LGMD-R12 [16,17] and FSHD [61,62]. Further evidence comes from animal models such as the ANO5-knock-out rabbit model, where scattered necrotic muscle fibers with inflammatory infiltrates are observed in muscle tissue [48].
The precise role of inflammation in the pathomechanism of muscular dystrophies is not well understood. Moreover, it is unclear if the inflammation is rather a primary feature or secondary to muscle degeneration. In contrast to acute muscle injury, the continuing injuries in dystrophic muscle result in the permanent recruitment of pro-inflammatory monocytes and the presence of pro-inflammatory macrophages in the muscle, leading to chronic inflammation [63]. Because of the asynchronicity of the injuries and of the regenerative cues, macrophages adopt a mixed phenotype making them less efficient in myogenesis, while they stimulate fibroadipogenic (FAP) cells (see below) and fibroblastic cells to produce extracellular matrix, leading to fat infiltrates and fibrosis. Chronic inflammation further worsens the disease, and can thus be considered as a secondary event in muscular dystrophies [55,56,64]. Therapeutic strategies inducing a change in macrophages that are present in the dystrophic muscle towards an anti-inflammatory profile could be beneficial, as has been shown in mdx-mice [64]. Another example of therapeutic intervention could be the depletion of osteopontin, which is an immunomodulator that is highly expressed in dystrophic muscles, as was shown in LGMD-R12 muscular dystrophy in our study. A recent study showed that osteopontin ablation ameliorated muscular dystrophy by shifting macrophages to a pro-regenerative phenotype in DMD [45].
In addition, our results showed an increased expression of serum amyloid A1 (SAA1) in dystrophic muscles in LGMD-R12 patients. Serum amyloid A protein (SAA) is an acute-phase protein, which is normally soluble and shows the highest concentration in plasma during inflammation, but which can be abnormally deposited in AA amyloidosis, secondary to chronic inflammatory processes, as fibers of insoluble protein in the extracellular space of various organs and tissues. Interestingly, muscle biopsies of some LGMD-R12 patients showed amyloid depositions around muscle fibers and within blood vessel walls [65,66]. Amyloid subtyping showed the presence of apoliprotein E, apoliprotein A1, apoliprotein A4, serum Amyloid P component and gelsolin in the deposits, whereas SAA was not typed in the depositions by Liewluck et al. [66]. Perhaps the technology was not sensitive enough to detect SAA in the amyloid deposits, or its absence could be explained by molecular interactions between amyloid precursors and associated conformational changes [67].

4.2. Fibroadipogenic Progenitor (FAP) Cells Are Upregulated in LGMD-R12 Muscles

Our data showed that fibroadipogenic progenitor cells (FAPs) are most highly expressed in the more severely affected muscles of LGMD-R12 patients. FAPs are muscle interstitial mesenchymal cells that interact with myogenic stem cells (satellite cells) to support myogenesis and muscle regeneration [68,69,70]. FAPs also have a role in muscle fibrosis and fatty tissue replacement, as is seen in muscular dystrophies [21,69,71,72]. FAPs undergo a large expansion, followed by their macrophage-mediated clearance and the reestablishment of their steady-state pool, during the first days following muscle injury. During this exact time window, FAPs establish a dynamic network of interactions via chemokine and interleukin signaling that culminate in muscle repair, together with the other cellular components of the muscle stem cell niche [73]. In muscular dystrophies, where dystrophic myofibers undergo repeated rounds of injury, the autocrine/paracrine constraints controlling FAP adipogenesis are released, which leads to fat infiltrates [74]. Duchenne muscular dystrophy (DMD) patients and mdx-mouse models show increased expression of platelet-derived growth factor receptor alpha (PDGFRα), which is an FAP marker. Furthermore, FAP accumulation correlates with increased fibro-fatty degeneration and disease severity [70,75]. High transforming growth factor beta (TGF-β) is an inflammatory factor produced in the dystrophic muscle that contributes to increased FAP accumulation and differentiation into fibrogenic and adipogenic fates [72,75]. Similar FAP involvement has also been shown in other muscular dystrophies, such as limb-girdle muscular dystrophy R2 (LGMD-R2) [76] and facioscapulohumeral muscular dystrophy (FSHD) [62,77]. Several clinical trials evaluating the effects of different drugs that alter FAP fate are being/were performed in DMD and LGMD-R2 [21]. Our deconvolution analysis suggested that, specifically, the recently described FAP lumican-positive (LUM+) subpopulation might be involved in the pathogenesis of LGMD-R12 [40]. Thus, our study indicates that targeting the FAP cell population may be therapeutically explored in LGMD-R12 patients. However, while these findings are interesting and can provide direction for future therapeutic interventions in LGMD-R12 patients, deconvolution analysis should be interpreted carefully and further single cell analyses should be performed to quantitatively assess differences in cell type composition.

4.3. Different Muscles Express Different Gene Profiles in Healthy Controls

Our study demonstrated that three different muscles from the thigh in healthy control individuals show different gene expression profiles. Furthermore, the dystrophic process in LGMD-R12 patients resulted in different pathological effects on these three muscles. These basic molecular differences between muscles might at least partially explain the selective involvement of distinct muscles and the absent involvement of others in muscular dystrophies, such as LGMD-R12. The variation in gene expression profiles is far beyond the differences in the proportion of conventional fiber types I (slow-twitch) and IIA/B (fast-twitch), and occurs even in muscles with similar compositions of fiber types, such as the rectus femoris, vastus lateralis and semimembranosus muscles used in our study, which are all muscles with fiber type II predominance [78,79]. More recent studies showed that although skeletal muscle in humans consists of three distinct fiber types (I, IIa, IIx), muscle contains many more isoforms of myosin heavy and light chains, which are coded by a large number of different genes [51,52]. Stuart and coworkers demonstrated that specific human skeletal muscle fiber types contain different mixtures of myosin heavy and light chains, with fast-twitch (type IIx) fibers consistently containing myosin heavy chains 1 (MYH1), 2 (MYH2) and 4 (MYH4) and myosin light chain 1 (MYL1), and slow-twitch (type I) fibers always containing myosin heavy chains 6 (MYH6) and 7 (MYH7) and myosin light chain 3 (MYL3) [51]. In our study, we showed a clear gradient of fast- and slow-twitch muscle genes across the healthy muscles, with the rectus femoris showing the highest expression of the fast-twitch muscle myosins MYH1, MYH2, MYH4 and MYL1 and the semimembranosus presenting the highest expression of the slow-twitch markers MYH6, MYH7 and MYL3. The vastus lateralis expressed intermediate levels of fast- and slow-twitch myosins. These differences in muscle fiber type isoforms might be one explanation of the differential involvement of muscles in muscular dystrophies, such as LGMD-R12.
Transcriptomic studies on different skeletal muscles in healthy animals also showed differential gene expressions between muscles. Terry et al. applied RNA-seq to profile RNA expression in 13 different skeletal muscles from mice and rats and showed extensive transcriptional diversity, with more than 50% of transcripts differentially expressed among skeletal muscles, an observation that cannot be explained by developmental history or fiber type composition alone [80]. The authors suggested a differential role of different transcription factors in distinct muscles to further explain their findings. Transcription factors recognize specific DNA sequences and control transcription by regulating the expression of genes in certain cells, at a certain time and in the right amount. Our study in humans demonstrated that several transcription factors were highly differentially expressed across the three healthy muscles, which may add to the explanation of the distinct susceptibility of different muscles in muscle diseases, such as LGMD-R12.

4.4. Conclusions and Future Perspectives

Overall, this study profiled the gene expression landscape in LGMD-R12. We identified baseline differences in expression levels between differently affected muscles and characterized the disease-associated changes. Our analysis showed differences in the cell composition of LGMD-R12 diseased and healthy muscle and suggested that transcription factor networks in specific cell populations underlie the differential predisposition of different muscles to degeneration. Because bulk RNA-sequencing only provides information on the average gene expression levels in muscle, we did not have the resolution to characterize these findings in more detail. Current single cell RNA-seq methodology would require a larger amount of muscle tissue per sample point, and was therefore not performed in this study. For further characterization at the cellular level, future studies using novel technologies such as single cell RNA-sequencing and possibly spatial RNA-sequencing need to be undertaken. Furthermore, this dataset can be used as a rich source for further exploration and referencing by the research community.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/cells11091508/s1, Supplementary Table S1: Overview of included samples and quality control in patients; Supplementary Table S2: Overview of included samples and quality control in control individuals; Supplementary Table S3: Top 50 of most differentially expressed genes across early stage patient muscles; Supplementary Figure S1: Re-analysis of scRNA-seq data. (A) Pie chart showing the cell type composition of the re-analyzed muscle single cell scRNA-seq data. (B) t-SNE analysis of the re-analyzed muscle scRNA-seq data, color coded for cell type. (C) Dendrogram visualization of hierarchical clustering and heatmap analysis on gene signature correlations between cell types. (D) Heatmap analysis of the top 10 marker genes of all cell types in the re-analyzed muscle scRNA-seq data. Abbreviations: vascular smooth muscle cells (vSMC), natural killer cells (NK). Supplementary Figure S2: Analysis of patient samples. (A) PCA plot of disease samples color coded for muscle of origin. Note that in contrast to healthy samples (see Figure 4A), disease samples do not group according to muscle. (B) PCA plot of disease samples color coded for Mercuri score. Note that disease samples group according to disease stage rather than muscle. (C) Heatmap analysis of the top 10 marker genes of all disease muscles.

Author Contributions

Conceptualization, C.E.D. and K.G.C.; methodology, C.E.D., V.G., R.J., A.D., N.N., S.D., D.R.T. and K.G.C.; software, C.E.D., R.J., S.D. and K.G.C.; validation, C.E.D., V.G., R.J., S.D., D.R.T. and K.G.C.; formal analysis, C.E.D., R.J., S.D. and K.G.C.; investigation, C.E.D., V.G., R.J., A.D., N.N., S.D. and K.G.C.; resources, J.L.D.B. and K.G.C.; data curation, C.E.D., R.J., S.D. and K.G.C.; writing—original draft preparation, C.E.D., R.J., S.D. and K.G.C.; writing—review and editing, C.E.D., V.G., R.J., A.D., J.L.D.B., N.N., S.D., D.R.T. and K.G.C.; visualization, C.E.D., V.G., R.J., N.N., S.D. and K.G.C.; supervision, S.D., D.R.T. and K.G.C.; project administration, C.E.D., A.D. and K.G.C.; funding acquisition, K.G.C. All authors have read and agreed to the published version of the manuscript.

Funding

This study was funded by the Association Française contre les Myopathies (AFM), trampoline grant number 21408.

Institutional Review Board Statement

The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Ethical Committee Research of UZ/KU Leuven (S-59867; date of approval: 19 April 2017).

Informed Consent Statement

Written informed consent was obtained from all patients and control individuals participating in the study.

Data Availability Statement

Data supporting reported results can be obtained from the corresponding author on reasonable request. RNA-seq datasets will be deposited in the GEO database.

Acknowledgments

We are grateful to the patients and control individuals for their participation in this study. We thank the personnel of the Departments of Pathology, Radiology and VIB Nucleomics Core for their technical support. We are thankful to the VIB Nucleomics Core (Leuven, Belgium, www.nucleomics.be, accessed on 15 July 2019) for library preparation and RNA-seq, and to Unicle Biomedical Data Science (Leuven, Belgium; https://unicle.life/, accessed on 8 September 2021) for data analysis using the online UniApp data analysis platform (v0.2.67). Some authors are member of the European Reference Network for Rare Neuromuscular Diseases (ERN EURO-NMD) and of the European Reference Network for Rare Neurological Diseases (ERN-RND).

Conflicts of Interest

C.E.D, V.G., R.J., A.D., J.D.B., N.N. and S.D. report no conflicts of interest. D.R.T. received speaker honorary or travel reimbursement from Novartis Pharma AG (Switzerland), GE Healthcare (UK) and UCB (Belgium), and collaborated with Novartis Pharma AG (Switzerland), Probiodrug (Germany), GE Healthcare (UK) and Janssen Pharmaceutical Companies (Belgium)—all not related to the current study. KGC participated at advisory boards and/or received speaker honorary from Alnylam, Amicus, ArgenX, Biogen, CSL Behring, Ipsen, Janssen Pharmaceutics, Lupin, Pfizer, Roche, Sanofi-Genzyme and UCB. KGC is Chairholder of the Emil von Behring Chair for Neuromuscular and Neurodegenerative Disorders by CSL Behring—all not related to the current study.

References

  1. Straub, V.; Murphy, A.; Udd, B.; LGMD Workshop Study Group. 229th ENMC international workshop: Limb girdle muscular dystrophies—Nomenclature and reformed classification Naarden, the Netherlands, 17–19 March 2017. Neuromuscul. Disord. 2018, 28, 702–710. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Bolduc, V.; Marlow, G.; Boycott, K.M.; Saleki, K.; Inoue, H.; Kroon, J.; Itakura, M.; Robitaille, Y.; Parent, L.; Baas, F.; et al. Recessive mutations in the putative calcium-activated chloride channel Anoctamin 5 cause proximal LGMD2L and distal MMD3 muscular dystrophies. Am. J. Hum. Genet. 2010, 86, 213–221. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Hicks, D.; Sarkozy, A.; Muelas, N.; Köehler, K.; Huebner, A.; Hudson, G.; Chinnery, P.F.; Barresi, R.; Eagle, M.; Polvikoski, T.; et al. A founder mutation in Anoctamin 5 is a major cause of limb-girdle muscular dystrophy. Brain 2011, 134, 171–182. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  4. Sarkozy, A.; Hicks, D.; Hudson, J.; Laval, S.H.; Barresi, R.; Hilton Jones, D.; Deschauer, M.; Harris, E.; Rufibach, L.; Hwang, E.; et al. ANO5 gene analysis in a large cohort of patients with anoctaminopathy: Confirmation of male prevalence and high occurrence of the common exon 5 gene mutation. Hum. Mutat. 2013, 34, 1111–1118. [Google Scholar] [CrossRef] [Green Version]
  5. Tsutsumi, S.; Kamata, N.; Vokes, T.J.; Maruoka, Y.; Nakakuki, K.; Enomoto, S.; Omura, K.; Amagasa, T.; Nagayama, M.; Saito-Ohara, F.; et al. The novel gene encoding a putative transmembrane protein is mutated in gnathodiaphyseal dysplasia (GDD). Am. J. Hum. Genet. 2004, 74, 1255–1261. [Google Scholar] [CrossRef] [Green Version]
  6. Mizuta, K.; Tsutsumi, S.; Inoue, H.; Sakamoto, Y.; Miyatake, K.; Miyawaki, K.; Noji, S.; Kamata, N.; Itakura, M. Molecular characterization of GDD1/TMEM16E, the gene product responsible for autosomal dominant gnathodiaphyseal dysplasia. Biochem. Biophys. Res. Commun. 2007, 357, 126–132. [Google Scholar] [CrossRef]
  7. Whitlock, J.M.; Hartzell, H.C. Anoctamins/TMEM16 Proteins: Chloride Channels Flirting with Lipids and Extracellular Vesicles. Annu. Rev. Physiol. 2017, 79, 119–143. [Google Scholar] [CrossRef] [Green Version]
  8. Boccaccio, A.; Di Zanni, E.; Gradogna, A.; Scholz-Starke, J. Lifting the veils on TMEM16E function. Channels 2019, 13, 33–35. [Google Scholar] [CrossRef] [Green Version]
  9. Griffin, D.A.; Johnson, R.W.; Whitlock, J.M.; Pozsgai, E.R.; Heller, K.N.; Grose, W.E.; Arnold, W.D.; Sahenk, Z.; Hartzell, H.C.; Rodino-Klapac, R. Defective membrane fusion and repair in Anoctamin5-deficient muscular dystrophy. Hum. Mol. Genet. 2016, 25, 1900–1911. [Google Scholar] [CrossRef] [Green Version]
  10. Whitlock, J.M.; Yu, K.; Cui, Y.Y.; Hartzell, H.C. Anoctamin 5/TMEM16E facilitates muscle precursor cell fusion. J. Gen. Physiol. 2018, 150, 1498–1509. [Google Scholar] [CrossRef] [Green Version]
  11. Chandra, G.; Sreetama, S.C.; Mazala, D.A.; Charton, K.; Van der Meulen, J.H.; Richard, I.; Jaiswal, J.K. Endoplasmic reticulum maintains ion homeostasis required for plasma membrane repair. J. Cell Biol. 2021, 220, e202006035. [Google Scholar] [CrossRef] [PubMed]
  12. Foltz, S.J.; Cui, Y.Y.; Choo, H.J.; Hartzell, H. ANO5 ensures trafficking of annexins in wounded myofibers. J. Cell Biol. 2021, 220, e202007059. [Google Scholar] [CrossRef] [PubMed]
  13. Thiruvengadam, G.; Sreetama, S.C.; Charton, K.; Hogarth, M.; Novak, J.S.; Suel-Petat, L.; Chandra, G.; Allard, B.; Richard, I.; Jaiswal, J.K. Anoctamin 5 Knockout Mouse Model Recapitulates LGMD2L Muscle Pathology and Offers Insight Into in vivo Functional Deficits. J. Neuromuscul. Dis. 2021, 8, S243–S255. [Google Scholar] [CrossRef] [PubMed]
  14. Dubowitz, V.; Sewry, C.A.; Oldfors, A. Muscle Biopsy E-Book: A Practical Approach, 4th ed.; Elsevier Health Sciences: New York, NY, USA, 2020. [Google Scholar]
  15. Silva, A.M.S.; Coimbra-Neto, A.R.; Souza, P.V.S.; Winckler, P.B.; Gonçalves, M.V.M.; Cavalcanti, E.B.U.; Carvalho, A.A.D.S.; Sobreira, C.F.D.R.; Camelo, C.G.; Mendonça, R.D.H.; et al. Clinical and molecular findings in a cohort of ANO5-related myopathy. Ann. Clin. Transl. Neurol. 2019, 6, 1225–1238. [Google Scholar] [CrossRef] [Green Version]
  16. Holm-Yildiz, S.; Witting, N.; de Stricker Borch, J.; Kass, K.; Khawajazada, T.; Krag, T.; Vissing, J. Muscle biopsy and MRI findings in ANO5-related myopathy. Muscle Nerve 2021, 64, 743–748. [Google Scholar] [CrossRef]
  17. Christiansen, J.; Güttsches, A.K.; Schara-Schmidt, U.; Vorgerd, M.; Heute, C.; Preusse, C.; Stenzel, W.; Roos, A. ANO5-related muscle diseases: From clinics and genetics to pathology and research strategies. Genes Diseases 2022. [Google Scholar] [CrossRef]
  18. Petrillo, S.; Pelosi, L.; Piemonte, F.; Travaglini, L.; Forcina, L.; Catteruccia, M.; Petrini, S.; Verardo, M.; D’Amico, A.; Musarò, A.; et al. Oxidative stress in Duchenne muscular dystrophy: Focus on the NRF2 redox pathway. Hum. Mol. Genet. 2017, 26, 2781–2790. [Google Scholar] [CrossRef]
  19. Vila, M.C.; Rayavarapu, S.; Hogarth, M.W.; Van der Meulen, J.H.; Horn, A.; Defour, A.; Takeda, S.; Brown, K.J.; Hathout, Y.; Nagaraju, K.; et al. Mitochondria mediate cell membrane repair and contribute to Duchenne muscular dystrophy. Cell Death Differ. 2017, 24, 330–342. [Google Scholar] [CrossRef] [Green Version]
  20. Tripodi, L.; Villa, C.; Molinaro, D.; Torrente, Y.; Farini, A. The Immune System in Duchenne Muscular Dystrophy Pathogenesis. Biomedicines 2021, 9, 1447. [Google Scholar] [CrossRef]
  21. Hogarth, M.W.; Uapinyoying, P.; Mázala, D.A.G.; Jaiswal, J.K. Pathogenic role and therapeutic potential of fibro-adipogenic progenitors in muscle disease. Trends Mol. Med. 2022, 28, 8–11. [Google Scholar] [CrossRef]
  22. Straub, V.; Carlier, P.G.; Mercuri, E. TREAT-NMD workshop: Pattern recognition in genetic muscle diseases using muscle MRI: 25-26 February 2011, Rome, Italy. Neuromuscul. Disord. 2012, 22, S42–S53. [Google Scholar] [CrossRef]
  23. Ten Dam, L.; van der Kooi, A.J.; Verhamme, C.; Wattjes, M.P.; de Visser, M. Muscle imaging in inherited and acquired muscle diseases. Eur. J. Neurol. 2016, 23, 688–703. [Google Scholar] [CrossRef] [PubMed]
  24. Mercuri, E.; Pichiecchio, A.; Allsop, J.; Messina, S.; Pane, M.; Muntoni, F. Muscle MRI in inherited neuromuscular disorders: Past, present, and future. J. Magn. Reson. Imaging 2007, 25, 433–440. [Google Scholar] [CrossRef] [PubMed]
  25. Kinali, M.; Arechavala-Gomeza, V.; Cirak, S.; Glover, A.; Guglieri, M.; Feng, L.; Hollingsworth, K.G.; Hunt, D.; Jungbluth, H.; Roper, H.P.; et al. Muscle histology vs. MRI in Duchenne muscular dystrophy. Neurology 2011, 76, 346–353. [Google Scholar] [CrossRef] [PubMed]
  26. Scotton, C.; Bovolenta, M.; Schwartz, E.; Falzarano, M.S.; Martoni, E.; Passarelli, C.; Armaroli, A.; Osman, H.; Rodolico, C.; Messina, S.; et al. Deep RNA profiling identified CLOCK and molecular clock genes as pathophysiological signatures in collagen VI myopathy. J. Cell Sci. 2016, 129, 1671–1684. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  27. Carson, J.A.; Nettleton, D.; Reecy, J.M. Differential gene expression in the rat soleus muscle during early work overload-induced hypertrophy. FASEB J. 2002, 16, 207–209. [Google Scholar] [CrossRef] [PubMed]
  28. Schmutz, S.; Däpp, C.; Wittwer, M.; Vogt, M.; Hoppeler, H.; Flück, M. Endurance training modulates the muscular transcriptome response to acute exercise. Pflugers Arch. 2006, 451, 678–687. [Google Scholar] [CrossRef] [Green Version]
  29. Cummings, B.B.; Marshall, J.L.; Tukiainen, T.; Lek, M.; Donkervoort, S.; Foley, A.R.; Bolduc, V.; Waddell, L.B.; Sandaradura, S.A.; O’Grady, G.L.; et al. Improving genetic diagnosis in Mendelian disease with transcriptome sequencing. Sci. Transl. Med. 2017, 9, eaal5209. [Google Scholar] [CrossRef] [Green Version]
  30. Morgan, M.; Anders, S.; Lawrence, M.; Aboyoun, P.; Pagès, H.; Gentleman, R. ShortRead: A bioconductor package for input, quality assessment and exploration of high-throughput sequence data. Bioinformatics 2009, 25, 2607–2608. [Google Scholar] [CrossRef] [Green Version]
  31. Langmead, B.; Salzberg, S.L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 2012, 9, 357–359. [Google Scholar] [CrossRef] [Green Version]
  32. Dobin, A.; Davis, C.A.; Schlesinger, F.; Drenkow, J.; Zaleski, C.; Jha, S.; Batut, P.; Chaisson, M.; Gingeras, T.R. STAR: Ultrafast universal RNA-seq aligner. Bioinformatics 2013, 29, 15–21. [Google Scholar] [CrossRef]
  33. Li, H.; Handsaker, B.; Wysoker, A.; Fennell, T.; Ruan, J.; Homer, N.; Marth, G.; Abecasis, G.; Durbin, R.; 1000 Genome Project Data Processing Subgroup. The Sequence Alignment/Map format and SAMtools. Bioinformatics 2009, 25, 2078–2079. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Liao, Y.; Smyth, G.K.; Shi, W. featureCounts: An efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 2014, 30, 923–930. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  35. Risso, D.; Schwartz, K.; Sherlock, G.; Dudoit, S. GC-content normalization for RNA-Seq data. BMC Bioinformatics 2011, 12, 480. [Google Scholar] [CrossRef] [Green Version]
  36. Satija, R.; Farrell, J.A.; Gennert, D.; Schier, A.F.; Regev, A. Spatial reconstruction of single cell gene expression data. Nat. Biotechnol. 2015, 33, 495–502. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  37. Ritchie, M.E.; Phipson, B.; Wu, D.; Hu, Y.; Law, C.W.; Shi, W.; Smyth, G.K. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015, 43, e47. [Google Scholar] [CrossRef] [PubMed]
  38. Yu, G.; Wang, L.G.; Han, Y.; He, Q.Y. clusterProfiler: An R package for comparing biological themes among gene clusters. OMICS 2012, 16, 284–287. [Google Scholar] [CrossRef] [PubMed]
  39. Wang, X.; Park, J.; Susztak, K.; Zhang, N.R.; Li, M. Bulk tissue cell type deconvolution with multi-subject single-cell expression reference. Nat. Commun. 2019, 10, 380. [Google Scholar] [CrossRef] [Green Version]
  40. Rubenstein, A.B.; Smith, G.R.; Raue, U.; Begue, G.; Minchev, K.; Ruf-Zamojski, F.; Nair, V.D.; Wang, X.; Zhou, L.; Zaslavsky, E.; et al. Single-cell transcriptional profiles in human skeletal muscle. Sci. Rep. 2020, 10, 229. [Google Scholar] [CrossRef] [Green Version]
  41. Hong, F.; Breitling, R.; McEntee, C.W.; Wittner, B.S.; Nemhauser, J.L.; Chory, J. RankProd: A bioconductor package for detecting differentially expressed genes in meta-analysis. Bioinformatics 2006, 22, 2825–2827. [Google Scholar] [CrossRef] [Green Version]
  42. Hack, A.A.; Groh, M.E.; McNally, E.M. Sarcoglycans in muscular dystrophy. Microsc. Res. Tech. 2000, 48, 167–180. [Google Scholar] [CrossRef]
  43. Csapo, R.; Gumpenberger, M.; Wessner, B. Skeletal Muscle Extracellular Matrix—What Do We Know About Its Composition, Regulation, and Physiological Roles? A Narrative Review. Front. Physiol. 2020, 11, 253. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  44. Krause, M.P.; Milne, K.J.; Hawke, T.J. Adiponectin-Consideration for its Role in Skeletal Muscle Health. Int. J. Mol. Sci. 2019, 20, 1528. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Capote, J.; Kramerova, I.; Martinez, L.; Vetrone, S.; Barton, E.R.; Sweeney, H.L.; Miceli, M.C.; Spencer, M.J. Osteopontin ablation ameliorates muscular dystrophy by shifting macrophages to a pro-regenerative phenotype. J. Cell Biol. 2016, 213, 275–288. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  46. Wasgewatte Wijesinghe, D.K.; Mackie, E.J.; Pagel, C.N. Normal inflammation and regeneration of muscle following injury require osteopontin from both muscle and non-muscle cells. Skelet. Muscle 2019, 9, 6. [Google Scholar] [CrossRef] [Green Version]
  47. Schiaffino, S.; Rossi, A.C.; Smerdu, V.; Leinwand, L.A.; Reggiani, C. Developmental myosins: Expression patterns and functional significance. Skelet. Muscle 2015, 5, 22. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  48. Sui, T.; Xu, L.; Lau, Y.S.; Liu, D.; Liu, T.; Gao, Y.; Lai, L.; Han, R.; Li, Z. Development of muscular dystrophy in a CRISPR-engineered mutant rabbit model with frame-disrupting ANO5 mutations. Cell Death Dis. 2018, 9, 609. [Google Scholar] [CrossRef]
  49. Hänzelmann, S.; Castelo, R.; Guinney, J. GSVA: Gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013, 14, 7. [Google Scholar] [CrossRef] [Green Version]
  50. Jin, H.; Liu, Z. A benchmark for RNA-seq deconvolution analysis under dynamic testing environments. Genome Biol. 2021, 22, 102. [Google Scholar] [CrossRef]
  51. Stuart, C.A.; Stone, W.L.; Howell, M.E.; Brannon, M.F.; Hall, H.K.; Gibson, A.L.; Stone, M.H. Myosin content of individual human muscle fibers isolated by laser capture microdissection. Am. J. Physiol. Cell Physiol. 2016, 310, C381–C389. [Google Scholar] [CrossRef] [Green Version]
  52. Lee, L.A.; Karabina, A.; Broadwell, L.J.; Leinwand, L.A. The ancient sarcomeric myosins found in specialized muscles. Skelet Muscle 2019, 9, 7. [Google Scholar] [CrossRef] [Green Version]
  53. Breitling, R.; Armengaud, P.; Amtmann, A.; Herzyk, P. Rank products: A simple, yet powerful, new method to detect differentially regulated genes in replicated microarray experiments. FEBS Lett. 2004, 573, 83–92. [Google Scholar] [CrossRef] [PubMed]
  54. Zoabi, M.; Zhang, L.; Li, T.M.; Elias, J.E.; Carlson, S.M.; Gozani, O. Methyltransferase-like 21C (METTL21C) methylates alanine tRNA synthetase at Lys-943 in muscle tissue. J. Biol. Chem. 2020, 295, 11822–11832. [Google Scholar] [CrossRef] [PubMed]
  55. Juban, G.; Saclier, M.; Yacoub-Youssef, H.; Kernou, A.; Arnold, L.; Boisson, C.; Ben Larbi, S.; Magnan, M.; Cuvellier, S.; Théret, M.; et al. AMPK Activation Regulates LTBP4-Dependent TGF-β1 Secretion by Pro-inflammatory Macrophages and Controls Fibrosis in Duchenne Muscular Dystrophy. Cell Rep. 2018, 25, 2163–2176.e6. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  56. Singh, P.; Chazaud, B. Benefits and pathologies associated with the inflammatory response. Exp. Cell Res. 2021, 409, 112905. [Google Scholar] [CrossRef]
  57. Spencer, M.J.; Tidball, J.G. Do immune cells promote the pathology of dystrophin-deficient myopathies? Neuromuscul. Disord. 2001, 11, 556–564. [Google Scholar] [CrossRef]
  58. Arahata, K.; Ishihara, T.; Fukunaga, H.; Orimo, S.; Lee, J.H.; Goto, K.; Nonaka, I. Inflammatory response in facioscapulohumeral muscular dystrophy (FSHD): Immunocytochemical and genetic analyses. Muscle Nerve 1995, 18, S56–S66. [Google Scholar] [CrossRef]
  59. Gallardo, E.; Rojas-García, R.; De Luna, N.; Pou, A.; Brown, R.H. Inflammation in dysferlin myopathy: Immunohistochemical characterization of 13 patients. Neurology 2001, 57, 2136–2138. [Google Scholar] [CrossRef]
  60. Panicucci, C.; Baratto, S.; Raffaghello, L.; Tonin, P.; D’Amico, A.; Tasca, G.; Traverso, M.; Fiorillo, C.; Minetti, C.; Previtali, S.C.; et al. Muscle inflammatory pattern in alpha- and gamma-sarcoglycanopathies. Clin. Neuropathol. 2021, 40, 310–318. [Google Scholar] [CrossRef]
  61. Dahlqvist, J.R.; Andersen, G.; Khawajazada, T.; Vissing, C.; Thomsen, C.; Vissing, J. Relationship between muscle inflammation and fat replacement assessed by MRI in facioscapulohumeral muscular dystrophy. J. Neurol. 2019, 266, 1127–1135. [Google Scholar] [CrossRef]
  62. van den Heuvel, A.; Lassche, S.; Mul, K.; Greco, A.; San León Granado, D.; Heerschap, A.; Küsters, B.; Tapscott, S.J.; Voermans, N.C.; van Engelen, B.G.M.; et al. Facioscapulohumeral dystrophy transcriptome signatures correlate with different stages of disease and are marked by different MRI biomarkers. Sci. Rep. 2022, 12, 1426. [Google Scholar] [CrossRef] [PubMed]
  63. Tidball, J.G.; Welc, S.S.; Wehling-Henricks, M. Immunobiology of Inherited Muscular Dystrophies. Compr. Physiol. 2018, 8, 1313–1356. [Google Scholar] [CrossRef] [PubMed]
  64. Saclier, M.; Ben Larbi, S.; My Ly, H.; Moulin, E.; Mounier, R.; Chazaud, B.; Juban, G. Interplay between myofibers and pro-inflammatory macrophages controls muscle damage in mdx mice. J. Cell Sci. 2021, 134, jcs258429. [Google Scholar] [CrossRef]
  65. Milone, M.; Liewluck, T.; Winder, T.L.; Pianosi, P.T. Amyloidosis and exercise intolerance in ANO5 muscular dystrophy. Neuromuscul. Disord. 2012, 22, 13–15. [Google Scholar] [CrossRef]
  66. Liewluck, T.; Winder, T.L.; Dimberg, E.L.; Crum, B.A.; Heppelmann, C.J.; Wang, Y.; Bergen, H.R.3rd; Milone, M. ANO5-muscular dystrophy: Clinical, pathological and molecular findings. Eur. J. Neurol. 2013, 20, 1383–1389. [Google Scholar] [CrossRef] [PubMed]
  67. Yang, M.; Liu, Y.; Dai, J.; Li, L.; Ding, X.; Xu, Z.; Mori, M.; Miyahara, H.; Sawashita, J.; Higuchi, K. Apolipoprotein A-II induces acute-phase response associated AA amyloidosis in mice through conformational changes of plasma lipoprotein structure. Sci. Rep. 2018, 8, 5620. [Google Scholar] [CrossRef]
  68. Joe, A.W.; Yi, L.; Natarajan, A.; Le Grand, F.; So, L.; Wang, J.; Rudnicki, M.A.; Rossi, F.M. Muscle injury activates resident fibro/adipogenic progenitors that facilitate myogenesis. Nat. Cell Biol. 2010, 12, 153–163. [Google Scholar] [CrossRef] [Green Version]
  69. Uezumi, A.; Ito, T.; Morikawa, D.; Shimizu, N.; Yoneda, T.; Segawa, M.; Yamaguchi, M.; Ogawa, R.; Matev, M.M.; Miyagoe-Suzuki, Y.; et al. Fibrosis and adipogenesis originate from a common mesenchymal progenitor in skeletal muscle. J. Cell Sci. 2011, 124, 3654–3664. [Google Scholar] [CrossRef] [Green Version]
  70. Uezumi, A.; Ikemoto-Uezumi, M.; Tsuchida, K. Roles of nonmyogenic mesenchymal progenitors in pathogenesis and regeneration of skeletal muscle. Front. Physiol. 2014, 5, 68. [Google Scholar] [CrossRef] [Green Version]
  71. Uezumi, A.; Fukada, S.; Yamamoto, N.; Takeda, S.; Tsuchida, K. Mesenchymal progenitors distinct from satellite cells contribute to ectopic fat cell formation in skeletal muscle. Nat. Cell Biol. 2010, 12, 143–152. [Google Scholar] [CrossRef]
  72. Giuliani, G.; Rosina, M.; Reggio, A. Signaling pathways regulating the fate of fibro/adipogenic progenitors (FAPs) in skeletal muscle regeneration and disease. FEBS J. 2021. [Google Scholar] [CrossRef]
  73. Biferali, B.; Proietti, D.; Mozzetta, C.; Madaro, L. Fibro-Adipogenic Progenitors Cross-Talk in Skeletal Muscle: The Social Network. Front. Physiol. 2019, 10, 1074. [Google Scholar] [CrossRef]
  74. Reggio, A.; Rosina, M.; Palma, A.; Cerquone Perpetuini, A.; Petrilli, L.L.; Gargioli, C.; Fuoco, C.; Micarelli, E.; Giuliani, G.; Cerretani, M.; et al. Adipogenesis of skeletal muscle fibro/adipogenic progenitors is affected by the WNT5a/GSK3/β-catenin axis. Cell Death Differ. 2020, 27, 2921–2941. [Google Scholar] [CrossRef] [PubMed]
  75. Mázala, D.A.; Novak, J.S.; Hogarth, M.W.; Nearing, M.; Adusumalli, P.; Tully, C.B.; Habib, N.F.; Gordish-Dressman, H.; Chen, Y.W.; Jaiswal, J.K.; et al. TGF-β-driven muscle degeneration and failed regeneration underlie disease onset in a DMD mouse model. JCI Insight 2020, 5, e135703. [Google Scholar] [CrossRef] [PubMed]
  76. Hogarth, M.W.; Defour, A.; Lazarski, C.; Gallardo, E.; Diaz Manera, J.; Partridge, T.A.; Nagaraju, K.; Jaiswal, J.K. Fibroadipogenic progenitors are responsible for muscle loss in limb girdle muscular dystrophy 2B. Nat. Commun. 2019, 10, 2430. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  77. Bosnakovski, D.; Shams, A.S.; Yuan, C.; da Silva, M.T.; Ener, E.T.; Baumann, C.W.; Lindsay, A.J.; Verma, M.; Asakura, A.; Lowe, D.A.; et al. Transcriptional and cytopathological hallmarks of FSHD in chronic DUX4-expressing mice. J. Clin. Investig. 2020, 130, 2465–2477. [Google Scholar] [CrossRef] [PubMed]
  78. Polgar, J.; Johnson, M.A.; Weightman, D.; Appleton, D. Data on fibre size in thirty-six human muscles. An autopsy study. J. Neurol. Sci. 1973, 19, 307–318. [Google Scholar] [CrossRef]
  79. Garrett, W.E., Jr.; Califf, J.C.; Bassett, F.H. 3rd. Histochemical correlates of hamstring injuries. Am. J. Sports Med. 1984, 12, 98–103. [Google Scholar] [CrossRef]
  80. Terry, E.E.; Zhang, X.; Hoffmann, C.; Hughes, L.D.; Lewis, S.A.; Li, J.; Wallace, M.J.; Riley, L.A.; Douglas, C.M.; Gutierrez-Monreal, M.A.; et al. Transcriptional profiling reveals extraordinary diversity among skeletal muscle tissues. eLife 2018, 7, e34613. [Google Scholar] [CrossRef]
Cells 11 01508 g001 550
Figure 1. Mercuri score is correlated with functional outcomes in LGMD-R12 patients. (A) Axial MRI section through the left thigh of patient 6, indicating the differential involvement of the target muscles: rectus femoris (RF) with Mercuri score 1, vastus lateralis (VL) with Mercuri score 3 and semimembranosus (SM) with Mercuri score 4. (B) A similar axial MRI slice of a control individual is shown. All muscles in healthy controls are per inclusion criterion normal and thus have a Mercuri score of 1. (C) Mean Mercuri score for all patients and for all controls is shown for the three target muscles. The bars indicate the 95% confidence interval. (D) A significant correlation between the mean Mercuri score per patient (dots) and the disease duration is shown in a scatter diagram, with “r” indicating the Pearson correlation coefficient, and “p” the p-value. (E,F) Scatter diagrams illustrating a significant increase in the Mercuri score with decreasing 6MWD (E) and with increasing 10MWT (F).
Figure 1. Mercuri score is correlated with functional outcomes in LGMD-R12 patients. (A) Axial MRI section through the left thigh of patient 6, indicating the differential involvement of the target muscles: rectus femoris (RF) with Mercuri score 1, vastus lateralis (VL) with Mercuri score 3 and semimembranosus (SM) with Mercuri score 4. (B) A similar axial MRI slice of a control individual is shown. All muscles in healthy controls are per inclusion criterion normal and thus have a Mercuri score of 1. (C) Mean Mercuri score for all patients and for all controls is shown for the three target muscles. The bars indicate the 95% confidence interval. (D) A significant correlation between the mean Mercuri score per patient (dots) and the disease duration is shown in a scatter diagram, with “r” indicating the Pearson correlation coefficient, and “p” the p-value. (E,F) Scatter diagrams illustrating a significant increase in the Mercuri score with decreasing 6MWD (E) and with increasing 10MWT (F).
Cells 11 01508 g001
Cells 11 01508 g002 550
Figure 2. Unbiased transcriptome analysis of all genes. (A) PCA of all samples (n = 84), color coded for Mercuri score. (B) PCA of all samples (n = 84), color coded for read count. Note that samples with the highest Mercuri score had overall lower sequencing depth, which can be explained by the fact that the cellularity in high Mercuri score samples is lower due to replacement of muscle cells with fibrotic and adipose tissue. If there is low cellularity, then we also have lower transcripts and hence fewer reads that map to those transcripts. (C) PCA of all samples included in the study (n = 84), color coded for feature (gene) count. Note that samples with the highest Mercuri score had overall highest feature (gene) count. (D) PCA of all samples included in the study (n = 84), color coded for mitochondrial feature (gene) count. (E) PCA of all samples included in the study (n = 84), color coded for passing (grey, n = 79) or failing (red, n = 5) quality control (QC). All samples with 1M or more reads were considered to have passed QC. (F) PCA of samples that passed QC (n = 79), color coded for read count. (G) PCA of samples that passed QC (n = 79), color coded for feature (gene) count. (H) PCA of samples that passed QC (n = 79), color coded for mitochondrial feature (gene) count. (I) PCA of samples that passed QC (n = 79), color coded for Mercuri score. (J) PCA of samples that passed QC (n = 79), color coded for SGCA gene expression. Note, SGCA is a protein of the muscle membrane and SGCA gene expression (muscle membrane marker) was reduced in samples with higher Mercuri scores. (K) PCA of samples that passed QC (n = 79), color coded for COL1A1 gene expression (fibroblast marker). (L) PCA of samples that passed QC (n = 79), color coded for ADIPOQ gene expression (adipocyte marker). In the PCA plots each sample is represented by a dot. The dots are placed on the graph in function of their gene signatures, the closer the dots are together, the more similar their gene expression profile is.
Figure 2. Unbiased transcriptome analysis of all genes. (A) PCA of all samples (n = 84), color coded for Mercuri score. (B) PCA of all samples (n = 84), color coded for read count. Note that samples with the highest Mercuri score had overall lower sequencing depth, which can be explained by the fact that the cellularity in high Mercuri score samples is lower due to replacement of muscle cells with fibrotic and adipose tissue. If there is low cellularity, then we also have lower transcripts and hence fewer reads that map to those transcripts. (C) PCA of all samples included in the study (n = 84), color coded for feature (gene) count. Note that samples with the highest Mercuri score had overall highest feature (gene) count. (D) PCA of all samples included in the study (n = 84), color coded for mitochondrial feature (gene) count. (E) PCA of all samples included in the study (n = 84), color coded for passing (grey, n = 79) or failing (red, n = 5) quality control (QC). All samples with 1M or more reads were considered to have passed QC. (F) PCA of samples that passed QC (n = 79), color coded for read count. (G) PCA of samples that passed QC (n = 79), color coded for feature (gene) count. (H) PCA of samples that passed QC (n = 79), color coded for mitochondrial feature (gene) count. (I) PCA of samples that passed QC (n = 79), color coded for Mercuri score. (J) PCA of samples that passed QC (n = 79), color coded for SGCA gene expression. Note, SGCA is a protein of the muscle membrane and SGCA gene expression (muscle membrane marker) was reduced in samples with higher Mercuri scores. (K) PCA of samples that passed QC (n = 79), color coded for COL1A1 gene expression (fibroblast marker). (L) PCA of samples that passed QC (n = 79), color coded for ADIPOQ gene expression (adipocyte marker). In the PCA plots each sample is represented by a dot. The dots are placed on the graph in function of their gene signatures, the closer the dots are together, the more similar their gene expression profile is.
Cells 11 01508 g002
Cells 11 01508 g003 550
Figure 3. Gene signatures in limb-girdle muscular dystrophy R12. (A) Volcano plot of a differential gene expression analysis of LGMD-R12 patient versus control samples. Genes mentioned in text are highlighted in red. (B) Bar plot showing top enriched gene sets based on the results of differential gene expression analysis of LGMD-R12 patient versus control samples. (C) Volcano plot of a differential gene expression analysis of Mercuri score 1 versus Mercuri score 4 semimembranosus samples. Genes associated with immune cells, collagen production, ribosomes, cell death, mitochondria and aerobic respiration are color coded. Note that pathways associated with skeletal muscles are downregulated (red dots), whereas stromal (collagen) (green dots) and immune cell signatures (blue dots) are upregulated. (D) Volcano plot of a differential pathway expression analysis of Mercuri score 1 versus Mercuri score 4 semimembranosus samples. Pathways associated with skeletal muscle, stromal and immune cells are color coded. Note that genes associated with skeletal muscles are downregulated, stromal and immune cell signatures are upregulated. (E) Heatmap analysis of pathway activity associated with Mercuri score. Note that pathways associated with myogenesis and typical muscle signatures (oxidative phosphorylation, glycolysis and protein metabolism) are progressively downregulated, while pathways associated with cell death, inflammation and mesenchymal transition are upregulated. (F) Heatmap analysis of the proteoglycan 4+ (PRG4+) fibroadipogenic progenitor (FAP) and lumican+ (LUM+) FAP marker genes’ expression. These are the markers from the single cell analysis that were used as input for deconvolution analysis in panel (G). (G) Bar plot quantifying the deconvolution analysis predictions on the relative percentage of the indicated cell types.
Figure 3. Gene signatures in limb-girdle muscular dystrophy R12. (A) Volcano plot of a differential gene expression analysis of LGMD-R12 patient versus control samples. Genes mentioned in text are highlighted in red. (B) Bar plot showing top enriched gene sets based on the results of differential gene expression analysis of LGMD-R12 patient versus control samples. (C) Volcano plot of a differential gene expression analysis of Mercuri score 1 versus Mercuri score 4 semimembranosus samples. Genes associated with immune cells, collagen production, ribosomes, cell death, mitochondria and aerobic respiration are color coded. Note that pathways associated with skeletal muscles are downregulated (red dots), whereas stromal (collagen) (green dots) and immune cell signatures (blue dots) are upregulated. (D) Volcano plot of a differential pathway expression analysis of Mercuri score 1 versus Mercuri score 4 semimembranosus samples. Pathways associated with skeletal muscle, stromal and immune cells are color coded. Note that genes associated with skeletal muscles are downregulated, stromal and immune cell signatures are upregulated. (E) Heatmap analysis of pathway activity associated with Mercuri score. Note that pathways associated with myogenesis and typical muscle signatures (oxidative phosphorylation, glycolysis and protein metabolism) are progressively downregulated, while pathways associated with cell death, inflammation and mesenchymal transition are upregulated. (F) Heatmap analysis of the proteoglycan 4+ (PRG4+) fibroadipogenic progenitor (FAP) and lumican+ (LUM+) FAP marker genes’ expression. These are the markers from the single cell analysis that were used as input for deconvolution analysis in panel (G). (G) Bar plot quantifying the deconvolution analysis predictions on the relative percentage of the indicated cell types.
Cells 11 01508 g003
Cells 11 01508 g004 550
Figure 4. Gene signatures across muscles in healthy control muscles. (A) PCA of healthy control samples color coded for muscle of origin. (B) PCA of healthy control samples color coded for MYH1 gene expression, which is characteristic for fast-twitch muscle fibers. Note the highest expression of MYH1 in the rectus femoris muscle. (C) PCA of healthy control samples color coded for MYH6 gene expression, which is characteristic for slow-twitch muscle fibers. Note the highest expression of MYH6 in the semimembranosus muscle. (D) Volcano plot of a differential gene expression analysis of rectus femoris versus semimembranosus samples of control individuals. Genes associated with the indicated biological themes are color coded. (E) Bar plot showing top enriched gene sets based on the results of differential gene expression analysis of rectus femoris versus semimembranosus control samples. Note that the epithelial to mesenchymal transition signature is higher expressed in the semimembranosus muscle, while pathways associated with glucose metabolism, DNA repair and protein metabolism are higher in the rectus femoris. (F) Heatmap analysis of the top 10 marker genes of the biopsied muscles of healthy control individuals. (G) Heatmap analysis of the top 20 cell type pathways associated with each of the three muscles in healthy controls. Several interesting cell type signatures have been indicated. (H) Heatmap analysis of the top 20 biological processes associated with each of the three healthy muscles. Several interesting cell type signatures have been indicated.
Figure 4. Gene signatures across muscles in healthy control muscles. (A) PCA of healthy control samples color coded for muscle of origin. (B) PCA of healthy control samples color coded for MYH1 gene expression, which is characteristic for fast-twitch muscle fibers. Note the highest expression of MYH1 in the rectus femoris muscle. (C) PCA of healthy control samples color coded for MYH6 gene expression, which is characteristic for slow-twitch muscle fibers. Note the highest expression of MYH6 in the semimembranosus muscle. (D) Volcano plot of a differential gene expression analysis of rectus femoris versus semimembranosus samples of control individuals. Genes associated with the indicated biological themes are color coded. (E) Bar plot showing top enriched gene sets based on the results of differential gene expression analysis of rectus femoris versus semimembranosus control samples. Note that the epithelial to mesenchymal transition signature is higher expressed in the semimembranosus muscle, while pathways associated with glucose metabolism, DNA repair and protein metabolism are higher in the rectus femoris. (F) Heatmap analysis of the top 10 marker genes of the biopsied muscles of healthy control individuals. (G) Heatmap analysis of the top 20 cell type pathways associated with each of the three muscles in healthy controls. Several interesting cell type signatures have been indicated. (H) Heatmap analysis of the top 20 biological processes associated with each of the three healthy muscles. Several interesting cell type signatures have been indicated.
Cells 11 01508 g004
Table
Table 1. Clinical and genetic features of the LGMD-R12 patients and MRI Mercuri scores of their biopsied muscles.
Table 1. Clinical and genetic features of the LGMD-R12 patients and MRI Mercuri scores of their biopsied muscles.
Patient NumberGenderAge at Symptom Onset
(y)		Age at
Study Inclusion
(y)		Disease Duration at Inclusion
(y)		6MWD (m)10MWT (s)Mercuri
Score
SM		Mercuri Score
VL		Mercuri Score
RF		ANO5
Mutations
1M3063333667.8442c.191dupA (p.Asn64Lysfs*15);
c.2317A>G (p.Met773Val)
2M313327853.5211c.191dupA (p.Asn64Lysfs*15);
c.191dupA (p.Asn64Lysfs*15)
3M343736894.5221c.191dupA (p.Asn64Lysfs*15);
c.191dupA (p.Asn64Lysfs*15)
4M303886324.8111c.191dupA (p.Asn64Lysfs*15);
c.1961G>A (p.Arg654Gln) and c.155A>G (p.Asn52Ser)
5M4759123698.1441c.172C>T (p.Arg58Trp);
c.172C>T (p.Arg58Trp)
6M3048184797.1431c.191dupA (p.Asn64Lysfs*15);
c.692G>T (p.Gly231Val)
7M3855177923.3211c.1213C>T (p.Gln405X);
c.1733T>C (p.Phe578Ser)
8M3964253628.1442c.191dupA (p.Asn64Lysfs*15);
c.191dupA (p.Asn64Lysfs*15)
9M354384527.5211c.1210C>T (p.Arg404X);
c.2387C>T (p.Ser796Leu)
10M3346135167.8432c.191dupA (p.Asn64Lysfs*15);
c.294+1G>A (p.?)
11M1548336305.0211c.649-2A>G (p.?);
c.679-2A>G (p.?)
12M3364312889.8444c.41-1G>A (p.?);
c.752C>T (p.Pro251Leu)
13M1326135705.0211c.191dupA (p.Asn64Lysfs*15);
c.191dupA (p.Asn64Lysfs*15)
14M343626924.9211c.191dupA (p.Asn64Lysfs*15);
c.242A>G (p.Asp81Gly)
15M2840124338.2421c.191dupA (p.Asn64Lysfs*15);
c.1213C>T (p.Gln405X)
16(#)M253166873.1211c.2411G>C (p.Cys804Ser);
c.1627dupA (p.Met543Asnfs*11)
y, years; 6MWD, six-minute walk distance; m, meter; 10MWT, 10-m walk test; s, seconds; SM, semimembranosus muscle; VL, vastus lateralis muscle; RF, rectus femoris muscle; ANO5, anoctamin-5 gene; M, male; 16(#), patient 16 was not included for RNA-seq because there was no RNA yield in the three muscle biopsies of this patient. The reference sequence to which the reported variants are referred is NM_213599.2.
Table
Table 2. Top 50 of most differentially expressed genes across healthy muscles.
Table 2. Top 50 of most differentially expressed genes across healthy muscles.
Marker Gene Analysis across Healthy Muscles
RankSemimembranosusRectus FemorisVastus Lateralis
1HAND2-AS1TBX5LINC02107
2C12orf75NTNG2LINC02119
3HAND2GCNT2SBK2
4HOXD8TYRP1RHOXF1-AS1
5FRMD1ANKRD36BP2C1orf158
6LBPMUC22CALML6
7HOXD-AS2HOXA13MYH1
8METTL21CZNF750LRRC37A7P
9SLC1A2LINC01854FAM184B
10HOXD9WFIKKN1RNA5-8S5
11IL31RATBX5-AS1LINC01886
12KIF1AFNDC10FEZF1-AS1
13ANGPTL8IPCEF1CRNDE
14CFAP57JCHAINCPXM1
15DNAH3CD300LBTCF24
16IL22RA1SATB2-AS1CDH22
17BDNFTUBB1PAX3
18LRRC52DIRAS1SNORD115-30
19C10orf67SPTA1GGT7
20CAPN8LINC01968MYHAS
21CROCC2HMGCS2RGS10
22LAD1TREM1SLITRK3
23ANKRD18BCCDC189PLCH1
24GLYATCOL9A1IGFN1
25SCDP2RX3ATRNL1
26IRX6RRM2ACTN3
27PAQR9-AS1IGHMSHISA2
28LINC01484IGHV3-7GADD45G
29C6orf132SPAG17MIR503HG
30HOXD3RPS27AP9FGF10
31TMC1PAX1GREM2
32FOSIRX4GDA
33SCRT1CXCR2P1OPRD1
34RNY4P10DNAH11RN7SL813P
35SLC26A9SIM2UBASH3B
36CCDC78TMEM163SNORD115-23
37COMPCDH20GDNF
38ADRB1TMEM105NANOS1
39PLPPR1SKA3LINC01773
40GPR39SLC7A11-AS1HSD52
41LINC00877EPB42NPR3
42SLC29A4SLC30A8NME9
43RSPO1HS6ST2CHAD
44BARX2FGD5P1MKRN3
45LOXL1-AS1SEL1L2SCT
46GPA33NLRP12RN7SL267P
47SDR42E2SLC4A10KHDRBS2
48CERS3CD160MYH4
49C2CD4BETF1P2RN7SL541P
50TRPM1RYR3SH2D1B
Table
Table 3. Top 50 of most differentially expressed genes across LGMD-R12 patient muscles.
Table 3. Top 50 of most differentially expressed genes across LGMD-R12 patient muscles.
Marker Gene Analysis across Patient Muscles
RankSemimembranosusRectus FemorisVastus Lateralis
1COMPLINC02107SIM2
2HAND2MYH1P2RX3
3HAND2-AS1LRRC37A7PHMGCS2
4COL20A1AQP4CHAC1
5AQP6LINC01773LINC01854
6ADIPOQLINC02119KCTD8
7PLEKHG4BC1orf158TBX5-AS1
8FRMD1PVALBZNF750
9CIDECACTN3IGLV3-21
10TNMDCALML6IGLC3
11HYDINMYLK4HOXA13
12SCDFBP2IGHD
13SALL1HCN1ADAMTS19-AS1
14LRRC74AUNC13CLAMC3
15LEPB3GALT1MAPT-AS1
16SLC1A6ERBB4FBP2
17PLA2G2AGREM2HAND2-AS1
18KCNQ2ATP2A1NEU4
19LGALS12RHOXF1-AS1SNCB
20CHI3L1LINC01886IL20RA
21GRM5NANOS1AQP4
22SLC5A10SHISA2HOXC12
23CCL18UGT3A1DIRAS1
24CUX2MLF1SLC16A3
25KLBNRG4TSHR
26MUC16SH2D1BGDNF
27GRIN2BPLCH1KLHDC7B
28PIEZO2LRRC3BLINC01018
29SAA1MYHASSLC51A
30MKXMYLK2GHRL
31GRM4FEZF1-AS1TBX1
32TSPEARFAM166BIGHV1-3
33TNCSNORD23TYRP1
34MYEOVENO3CRYM
35PCK1ENSAP2PIANP
36PLIN1IRX3IGKV1-5
37CACNA1ICDH22TMEM26
38USH2AATRNL1MTND4P24
39MUC6LANCL1-AS1RN7SKP276
40COL22A1NEK10FAM166B
41SCUBE1PDE4DIPP1HPN
42SCRG1TMEM266OSCAR
43S100A3FABP7HES7
44KRT7KLHL38PAX1
45OPRM1AGMATASB12
46DUX4L19SMCO1DDX11L2
47MYBL2ASB14ANKRD20A21P
48AMZ1NPSR1-AS1SPAG17
49RASAL1PITX1TNNI3
50MROH4PDDIT4LC1orf105
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Depuydt, C.E.; Goosens, V.; Janky, R.; D’Hondt, A.; De Bleecker, J.L.; Noppe, N.; Derveaux, S.; Thal, D.R.; Claeys, K.G. Unraveling the Molecular Basis of the Dystrophic Process in Limb-Girdle Muscular Dystrophy LGMD-R12 by Differential Gene Expression Profiles in Diseased and Healthy Muscles. Cells 2022, 11, 1508. https://doi.org/10.3390/cells11091508

AMA Style

Depuydt CE, Goosens V, Janky R, D’Hondt A, De Bleecker JL, Noppe N, Derveaux S, Thal DR, Claeys KG. Unraveling the Molecular Basis of the Dystrophic Process in Limb-Girdle Muscular Dystrophy LGMD-R12 by Differential Gene Expression Profiles in Diseased and Healthy Muscles. Cells. 2022; 11(9):1508. https://doi.org/10.3390/cells11091508

Chicago/Turabian Style

Depuydt, Christophe E., Veerle Goosens, Rekin’s Janky, Ann D’Hondt, Jan L. De Bleecker, Nathalie Noppe, Stefaan Derveaux, Dietmar R. Thal, and Kristl G. Claeys. 2022. "Unraveling the Molecular Basis of the Dystrophic Process in Limb-Girdle Muscular Dystrophy LGMD-R12 by Differential Gene Expression Profiles in Diseased and Healthy Muscles" Cells 11, no. 9: 1508. https://doi.org/10.3390/cells11091508

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop