Thank you for visiting nature.com. The browser version you are using has limited CSS support. For the best experience, we recommend using the latest browser version (or disabling compatibility mode in Internet Explorer). Additionally, to ensure continued support, this site will be free of styles and JavaScript.
Skeletal muscle is a heterogeneous tissue composed predominantly of myofibrils, which in humans are typically classified into three types: one “slow” (type 1) and two “fast” (types 2A and 2X). However, heterogeneity between and within traditional myofibril types remains poorly understood. We applied transcriptomic and proteomic approaches to 1050 and 1038 individual myofibrils from human vastus lateralis, respectively. The proteomic study included men, and the transcriptomic study included 10 men and 2 women. In addition to myosin heavy chain isoforms, we identified metabolic proteins, ribosomal proteins, and cellular junctional proteins as sources of multidimensional intermyofibril variability. Furthermore, despite the identification of clusters of slow and fast fibers, our data suggest that type 2X fibers are phenotypically indistinguishable from other fast-twitch fibers. Furthermore, myosin heavy chain-based classification is insufficient to describe the myofiber phenotype in nemaline myopathies. Overall, our data suggest multidimensional myofiber heterogeneity, with sources of variation extending beyond myosin heavy chain isoforms.
Cellular heterogeneity is an inherent feature of all biological systems, allowing cells to specialize to meet the different needs of tissues and cells.1 The traditional view of skeletal muscle fiber heterogeneity has been that motor neurons define the fiber type within a motor unit, and that fiber type (i.e., type 1, type 2A, and type 2X in humans) is determined by the characteristics of myosin heavy chain (MYH) isoforms.2 This was initially based on their pH ATPase instability,3,4 and later on their molecular expression of MYH.5 However, with the identification and subsequent acceptance of “mixed” fibers that co-express multiple MYHs in varying proportions, skeletal muscle fibers are increasingly viewed as a continuum rather than as distinct fiber types.6 Despite this, the field still relies heavily on MYH as the primary classifier for myofiber classification, a view likely influenced by the limitations and significant biases of early rodent studies whose MYH expression profiles and range of fiber types differ from those in humans.2 The situation is further complicated by the fact that different human skeletal muscles exhibit a diverse range of fiber types.7 The vastus lateralis is a mixed muscle with an intermediate (and therefore representative) MYH expression profile.7 Furthermore, its ease of sampling makes it the best-studied muscle in humans.
Thus, unbiased investigation of skeletal muscle fiber diversity using powerful “omics” tools is critical but also challenging, in part due to the multinucleated nature of skeletal muscle fibers. However, transcriptomics8,9 and proteomics10 technologies have undergone a revolution in sensitivity in recent years due to various technological advances, allowing the analysis of skeletal muscle at single-fiber resolution. As a result, significant progress has been made in characterizing single-fiber diversity and their response to atrophic stimuli and aging11,12,13,14,15,16,17,18. Importantly, these technological advances have clinical applications, allowing for more detailed and precise characterization of disease-associated dysregulation. For example, the pathophysiology of nemaline myopathy, one of the most common inherited muscle diseases (MIM 605355 and MIM 161800), is complex and confusing.19,20 Therefore, a better characterization of the dysregulation of skeletal muscle fibers could lead to significant advances in our understanding of this disease.
We developed methods for transcriptomic and proteomic analysis of single skeletal muscle fibers manually isolated from human biopsy specimens and applied them to thousands of fibers, allowing us to investigate the cellular heterogeneity of human skeletal muscle fibers. In the course of this work, we demonstrated the power of transcriptomic and proteomic phenotyping of muscle fibers and identified metabolic, ribosomal, and cellular junctional proteins as significant sources of interfiber variability. Furthermore, using this proteomic workflow, we characterized the clinical relevance of nematode myopathy in single skeletal muscle fibers, revealing a coordinated shift toward non-oxidative fibers independent of fiber type based on MYH.
To investigate the heterogeneity of human skeletal muscle fibers, we developed two workflows to enable transcriptome and proteome analysis of single skeletal muscle fibers (Figure 1A and Supplementary Figure 1A). We developed and optimized several methodological steps, from sample storage and preservation of RNA and protein integrity to optimizing throughput for each approach. For transcriptome analysis, this was achieved by inserting sample-specific molecular barcodes at the initial step of reverse transcription, allowing 96 fibers to be pooled for efficient downstream processing. Deeper sequencing (±1 million reads per fiber) compared to traditional single-cell approaches further enriched the transcriptome data. 21 For proteomics, we used a short chromatographic gradient (21 minutes) combined with DIA-PASEF data acquisition on a timsTOF mass spectrometer to optimize proteome depth while maintaining high throughput. 22,23 To investigate the heterogeneity of healthy skeletal muscle fibers, we characterized the transcriptomes of 1,050 individual fibers from 14 healthy adult donors and the proteomes of 1,038 fibers from 5 healthy adult donors (Supplementary Table 1). In this paper, these datasets are referred to as the 1,000-fiber transcriptomes and proteomes, respectively. Our approach detected a total of 27,237 transcripts and 2,983 proteins in the 1,000-fiber transcriptomic and proteomic analyses (Figure 1A, Supplementary Datasets 1–2). After filtering the transcriptomic and proteomic datasets for >1,000 detected genes and 50% valid values per fiber, subsequent bioinformatics analyses were performed for 925 and 974 fibers in the transcriptome and proteome, respectively. After filtering, an average of 4257 ± 1557 genes and 2015 ± 234 proteins (mean ± SD) were detected per fiber, with limited inter-individual variability (Supplementary Figures 1B–C, Supplementary Datasets 3–4). However, within-subject variability was more pronounced across participants, likely due to differences in RNA/protein yield between fibers of different lengths and cross-sectional areas. For most proteins (>2000), the coefficient of variation was below 20% (Supplementary Figure 1D). Both methods allowed capturing a wide dynamic range of transcripts and proteins with highly expressed signatures important for muscle contraction (e.g., ACTA1, MYH2, MYH7, TNNT1, TNNT3) (Supplementary Figures 1E–F). Most of the identified features were common between the transcriptomic and proteomic datasets (Supplementary Figure 1G), and the mean UMI/LFQ intensities of these features were reasonably well correlated (r = 0.52) (Supplementary Figure 1H).
Transcriptomics and proteomics workflow (created with BioRender.com). BD Dynamic range curves for MYH7, MYH2, and MYH1, and calculated thresholds for fiber type assignment. E, F Distribution of MYH expression across fibers in transcriptomics and proteomics datasets. G, H Uniform Diversity Approximation and Projection (UMAP) plots for transcriptomics and proteomics colored by MYH-based fiber type. I, J Feature plots showing MYH7, MYH2, and MYH1 expression in transcriptomics and proteomics datasets.
We initially set out to assign MYH-based fiber type to each fiber using an optimized approach that takes advantage of the high sensitivity and dynamic range of MYH expression in omics datasets. Previous studies have used arbitrary thresholds to label fibers as pure type 1, type 2A, type 2X, or mixed based on a fixed percentage of expression of different MYHs11,14,24. We used a different approach in which the expression of each fiber was ranked by the MYHs we used to type the fibers: MYH7, MYH2, and MYH1, corresponding to type 1, type 2A, and type 2X fibers, respectively. We then mathematically calculated the lower inflection point of each resulting curve and used it as a threshold to assign fibers as positive (above the threshold) or negative (below the threshold) for each MYH (Figure 1B–D). These data demonstrate that MYH7 (Figure 1B) and MYH2 (Figure 1C) have more distinct on/off expression profiles at the RNA level compared to the protein level. Indeed, at the protein level, very few fibers did not express MYH7, and no fiber had 100% MYH2 expression. We next used predetermined expression thresholds to assign MYH-based fiber types to all fibers in each dataset. For example, MYH7+/MYH2-/MYH1- fibers were assigned to type 1, while MYH7-/MYH2+/MYH1+ fibers were assigned to mixed type 2A/2X (see Supplementary Table 2 for full description). Pooling all fibers, we observed a remarkably similar distribution of MYH-based fiber types at both the RNA (Figure 1E) and protein (Figure 1F) levels, while the relative composition of MYH-based fiber types varied across individuals, as expected (Supplementary Figure 2A). Most fibers were classified as either pure type 1 (34–35%) or type 2A (36–38%), although a significant number of mixed type 2A/2X fibers were also detected (16–19%). A striking difference is that pure type 2X fibers could only be detected at the RNA level, but not at the protein level, suggesting that rapid MYH expression is at least partly regulated post-transcriptionally.
We validated our proteomics-based MYH fiber typing method using antibody-based dot blotting, and both methods achieved 100% agreement in identifying pure type 1 and type 2A fibers (see Supplementary Figure 2B). However, the proteomics-based approach was more sensitive, more efficient in identifying mixed fibers, and quantifying the proportion of each MYH gene in each fiber. These data demonstrate the effectiveness of using an objective, highly sensitive omics-based approach to characterize skeletal muscle fiber types.
We then used the combined information provided by transcriptomics and proteomics to objectively classify myofibers based on their complete transcriptome or proteome. Using the uniform manifold approximation and projection (UMAP) method to reduce the dimensionality to six principal components (Supplementary Figures 3A–B), we were able to visualize myofiber variability in the transcriptome (Figure 1G) and proteome (Figure 1H). Notably, myofibers were not grouped by participants (Supplementary Figures 3C–D) or testing days (Supplementary Figure 3E) in either the transcriptomics or proteomics datasets, suggesting that within-subject variability in skeletal muscle fibers is higher than between-subject variability. In the UMAP plot, two distinct clusters representing “fast” and “slow” myofibers emerged (Figures 1G–H). MYH7+ (slow) myofibers were clustered at the positive pole of UMAP1, whereas MYH2+ and MYH1+ (fast) myofibers were clustered at the negative pole of UMAP1 (Figures 1I–J). However, no distinction was made between fast-twitch fiber types (i.e., type 2A, type 2X, or mixed 2A/2X) based on MYH expression, suggesting that MYH1 (Figure 1I–J) or other classical 2X myofiber markers such as ACTN3 or MYLK2 (Supplementary Figures 4A–B) expression does not differentiate between different myofiber types when considering the entire transcriptome or proteome. Moreover, compared with MYH2 and MYH7, few transcripts or proteins were positively correlated with MYH1 (Supplementary Figs. 4C–H), suggesting that MYH1 abundance does not fully reflect the myofiber transcriptome/proteome. Similar conclusions were reached when assessing the mixed expression of the three MYH isoforms at the UMAP level (Supplementary Figs. 4I–J). Thus, while 2X fibers can be identified at the transcript level based on MYH quantification alone, MYH1+ fibers are not distinguishable from other fast fibers when considering the entire transcriptome or proteome.
As an initial exploration of slow fiber heterogeneity beyond MYH, we assessed four established slow fiber type-specific proteins: TPM3, TNNT1, MYL3, and ATP2A22. Slow fiber subtypes showed high, although not perfect, Pearson correlations with MYH7 in both transcriptomics (Supplementary Figure 5A) and proteomics (Supplementary Figure 5B). Approximately 25% and 33% of slow fibers were not classified as pure slow fibers by all gene/protein subtypes in transcriptomics (Supplementary Figure 5C) and proteomics (Supplementary Figure 5D), respectively. Therefore, slow fiber classification based on multiple gene/protein subtypes introduces additional complexity, even for proteins known to be fiber type-specific. This suggests that fiber classification based on isoforms of a single gene/protein family may not adequately reflect the true heterogeneity of skeletal muscle fibers.
To further explore the phenotypic variability of human skeletal muscle fibers at the scale of the whole omics model, we performed unbiased dimensionality reduction of the data using principal component analysis (PCA) (Figure 2A). Similar to UMAP plots, neither participant nor testing day influenced fiber clustering at the PCA level (Supplementary Figures 6A–C). In both datasets, MYH-based fiber type was explained by PC2, which showed a cluster of slow-twitch type 1 fibers and a second cluster containing fast-twitch type 2A, type 2X, and mixed 2A/2X fibers (Figure 2A). In both datasets, these two clusters were connected by a small number of mixed type 1/2A fibers. As expected, overrepresentation analysis of the main PC drivers confirmed that PC2 was driven by contractile and metabolic signatures (Figure 2B and Supplementary Figures 6D–E, Supplementary Datasets 5–6). Overall, MYH-based fiber type was found to be sufficient to explain continuous variation along PC2, with the exception of so-called 2X fibers that were distributed throughout the transcriptome within the fast cluster.
A. Principal component analysis (PCA) plots of transcriptome and proteome datasets colored according to fiber type based on MYH. B. Enrichment analysis of transcript and protein drivers in PC2 and PC1. Statistical analysis was performed using the clusterProfiler package and Benjamini-Hochberg adjusted p-values. C, D. PCA plots colored according to intercellular adhesion gene ontology (GO) terms in the transcriptome and costamere GO terms in the proteome. Arrows represent transcript and protein drivers and their directions. E, F. Uniform manifold approximation and projection (UMAP) feature plots of clinically relevant features showing expression gradients independent of slow/fast fiber type. G, H. Correlations between PC2 and PC1 drivers in transcriptomes and proteomes.
Unexpectedly, MYH-based myofiber type explained only the second highest degree of variability (PC2), suggesting that other biological factors unrelated to MYH-based myofiber type (PC1) play an important role in regulating skeletal muscle fiber heterogeneity. Overrepresentation analysis of the top drivers in PC1 revealed that variability in PC1 was primarily determined by cell-cell adhesion and ribosome content in the transcriptome, and costameres and ribosomal proteins in the proteome (Figure 2B and Supplementary Figures 6D–E, Supplementary Data Set 7). In skeletal muscle, costameres connect the Z-disc to the sarcolemma and are involved in force transmission and signaling. 25 Annotated PCA plots using cell-cell adhesion (transcriptome, Figure 2C) and costamere (proteome, Figure 2D) features revealed a strong left shift in PC1, indicating that these features are enriched in certain fibers.
A more detailed examination of myofiber clustering at the UMAP level revealed that most features exhibited a myofiber type-independent MYH-based expression gradient rather than myofiber subcluster-specific. This continuity was observed for several genes associated with pathological conditions (Figure 2E), such as CHCHD10 (neuromuscular disease), SLIT3 (muscle atrophy), CTDNEP1 (muscle disease). This continuity was also observed across the proteome, including proteins associated with neurological disorders (UGDH), insulin signaling (PHIP), and transcription (HIST1H2AB) (Figure 2F). Collectively, these data indicate continuity in fiber type-independent slow/fast twitch heterogeneity across different myofibers.
Interestingly, driver genes in PC2 showed good transcriptome-proteome correlation (r = 0.663) (Figure 2G), suggesting that slow- and fast-twitch fiber types, and in particular the contractile and metabolic properties of skeletal muscle fibers, are transcriptionally regulated. However, driver genes in PC1 showed no transcriptome-proteome correlation (r = -0.027) (Figure 2H), suggesting that variations unrelated to slow/fast-twitch fiber types are largely regulated post-transcriptionally. Because variations in PC1 were primarily explained by ribosomal gene ontology terms, and given that ribosomes play a crucial and specialized role in the cell by actively participating in and influencing protein translation,31 we next set out to investigate this unexpected ribosomal heterogeneity.
We first colored the proteomics principal component analysis plot according to the relative abundance of proteins in the GOCC term “cytoplasmic ribosome” (Figure 3A). Although this term is enriched on the positive side of PC1, resulting in a small gradient, ribosomal proteins drive partitioning in both directions of PC1 (Figure 3A). Ribosomal proteins enriched on the negative side of PC1 included RPL18, RPS18, and RPS13 (Figure 3B), while RPL31, RPL35, and RPL38 (Figure 3C) were the main drivers on the positive side of PC1. Interestingly, RPL38 and RPS13 were highly expressed in skeletal muscle compared to other tissues (Supplementary Figure 7A). These distinctive ribosomal signatures in PC1 were not observed in the transcriptome (Supplementary Figure 7B), indicating post-transcriptional regulation.
A. Principal component analysis (PCA) plot colored according to cytoplasmic ribosomal gene ontology (GO) terms across the proteome. Arrows indicate the direction of protein-mediated variation in the PCA plot. Line length corresponds to the principal component score for a given protein. B, C. PCA feature plots for RPS13 and RPL38. D. Unsupervised hierarchical clustering analysis of cytoplasmic ribosomal proteins. E. Structural model of the 80S ribosome (PDB: 4V6X) highlighting ribosomal proteins with different abundances in skeletal muscle fibers. F. Ribosomal proteins with different stoichiometry localized near the mRNA exit channel.
The concepts of ribosomal heterogeneity and specialization have been proposed previously, whereby the presence of distinct ribosome subpopulations (ribosomal heterogeneity) can directly influence protein translation in different tissues32 and cells33 through the selective translation of specific mRNA transcript pools34 (ribosome specialization). To identify subpopulations of ribosomal proteins co-expressed in skeletal muscle fibers, we performed an unsupervised hierarchical clustering analysis of ribosomal proteins in the proteome (Figure 3D, Supplementary Data Set 8). As expected, ribosomal proteins did not cluster by fiber type based on MYH. However, we identified three distinct clusters of ribosomal proteins; the first cluster (ribosomal_cluster_1) is coregulated with RPL38 and hence has increased expression in fibers with a positive PC1 profile. The second cluster (ribosomal_cluster_2) is coregulated with RPS13 and is elevated in fibers with a negative PC1 profile. The third cluster (ribosomal_cluster_3) does not show coordinated differential expression in skeletal muscle fibers and can be considered the “core” skeletal muscle ribosomal protein. Both ribosomal clusters 1 and 2 contain ribosomal proteins that have previously been shown to regulate alternative translation (e.g., RPL10A, RPL38, RPS19, and RPS25) and functionally influence development (e.g., RPL10A, RPL38).34,35,36,37,38 Consistent with the PCA results, the observed heterogeneous representation of these ribosomal proteins across fibers also showed continuity (Supplementary Figure 7C).
To visualize the location of heterogeneous ribosomal proteins within the ribosome, we used a structural model of the human 80S ribosome (Protein Data Bank: 4V6X) (Figure 3E). After isolating ribosomal proteins belonging to different ribosomal clusters, their locations were not closely aligned, suggesting that our approach failed to provide enrichment for certain regions/fractions of the ribosome. Interestingly, however, the proportion of large subunit proteins in cluster 2 was lower than in clusters 1 and 3 (Supplementary Figure 7D). We observed that proteins with altered stoichiometry in skeletal muscle fibers were predominantly localized to the ribosome surface (Figure 3E), consistent with their ability to interact with internal ribosome entry site (IRES) elements in different mRNA populations, thereby coordinating selective translation. 40, 41 Furthermore, many proteins with altered stoichiometry in skeletal muscle fibers were located near functional regions such as the mRNA exit tunnel (Figure 3F), which selectively regulate translational elongation and arrest of specific peptides. 42 In summary, our data suggest that the stoichiometry of skeletal muscle ribosomal proteins exhibits heterogeneity, resulting in differences between skeletal muscle fibers.
We next set out to identify fast- and slow-twitch fiber signatures and explore the mechanisms of their transcriptional regulation. Comparing the fast- and slow-twitch fiber clusters defined by UMAP in the two datasets (Figures 1G–H and 4A–B), transcriptomic and proteomic analyses identified 1366 and 804 differentially abundant features, respectively (Figures 4A–B, Supplementary Datasets 9–12). We observed the expected differences in signatures related to sarcomeres (e.g., tropomyosin and troponin), excitation-contraction coupling (SERCA isoforms), and energy metabolism (e.g., ALDOA and CKB). Moreover, transcripts and proteins regulating protein ubiquitination were differentially expressed in fast- and slow-twitch fibers (e.g., USP54, SH3RF2, USP28, and USP48) (Figures 4A–B). Moreover, the microbial protein gene RP11-451G4.2 (DWORF), which has previously been shown to be differentially expressed across lamb muscle fiber types43 and enhances SERCA activity in cardiac muscle44, was significantly upregulated in slow skeletal muscle fibers (Figure 4A). Similarly, at the individual fiber level, significant differences were observed in known signatures such as metabolism-related lactate dehydrogenase isoforms (LDHA and LDHB, Figure 4C and Supplementary Figure 8A)45,46 as well as previously unknown fiber-type-specific signatures (such as IRX3, USP54, USP28, and DPYSL3) (Figure 4C). There was a significant overlap of differentially expressed features between the transcriptomic and proteomic datasets (Supplementary Figure 8B), as well as a fold change correlation driven mainly by the more pronounced differential expression of sarcomere features (Supplementary Figure 8C). Notably, some signatures (e.g. USP28, USP48, GOLGA4, AKAP13) showed strong post-transcriptional regulation only at the proteomic level and had slow/fast twitch fiber type-specific expression profiles (Supplementary Figure 8C).
A and B volcano plots comparing slow and fast clusters identified by the uniform manifold approximation and projection (UMAP) plots in Figures 1G–H. Colored dots represent transcripts or proteins that are significantly different at FDR < 0.05, and darker dots represent transcripts or proteins that are significantly different at log change > 1. Two-way statistical analysis was performed using the DESeq2 Wald test with Benjamini–Hochberg adjusted p values (transcriptomics) or the Limma linear model method with empirical Bayesian analysis followed by Benjamini–Hochberg adjustment for multiple comparisons (proteomics). C Signature plots of selected differentially expressed genes or proteins between slow and fast fibers. D Enrichment analysis of significantly differentially expressed transcripts and proteins. Overlapping values are enriched in both datasets, transcriptome values are enriched only in the transcriptome, and proteome values are enriched only in the proteome. Statistical analysis was performed using the clusterProfiler package with Benjamini-Hochberg adjusted p-values. E. Fiber type-specific transcription factors identified by SCENIC based on SCENIC-derived regulator specificity scores and differential mRNA expression between fiber types. F. Profiling of selected transcription factors differentially expressed between slow and fast fibers.
We then performed an overrepresentation analysis of differentially represented genes and proteins (Figure 4D, Supplemental Data Set 13). Pathway enrichment for features that differed between the two datasets revealed expected differences, such as fatty acid β-oxidation and ketone metabolism processes (slow fibers), myofilament/muscle contraction (fast and slow fibers, respectively), and carbohydrate catabolic processes (fast fibers). Serine/threonine protein phosphatase activity was also elevated in fast fibers, driven by features such as the regulatory and catalytic phosphatase subunits (PPP3CB, PPP1R3D, and PPP1R3A), which are known to regulate glycogen metabolism (47) (Supplemental Figures 8D–E). Other pathways enriched in fast fibers included processing (P-) bodies (YTHDF3, TRIM21, LSM2) in the proteome (Supplementary Fig. 8F), potentially involved in post-transcriptional regulation (48), and transcription factor activity (SREBF1, RXRG, RORA) in the transcriptome (Supplementary Fig. 8G). Slow fibers were enriched in oxidoreductase activity (BDH1, DCXR, TXN2) (Supplementary Fig. 8H), amide binding (CPTP, PFDN2, CRYAB) (Supplementary Fig. 8I), extracellular matrix (CTSD, ADAMTSL4, LAMC1) (Supplementary Fig. 8J), and receptor-ligand activity (FNDC5, SPX, NENF) (Supplementary Fig. 8K).
To gain further insight into the transcriptional regulation underlying slow/fast muscle fiber type characteristics, we performed transcription factor enrichment analysis using SCENIC49 (Supplementary Data Set 14). Many transcription factors were significantly enriched between fast and slow muscle fibers (Figure 4E). This included transcription factors such as MAFA, which has previously been linked to fast muscle fiber development,50 as well as several transcription factors not previously associated with muscle fiber type-specific gene programs. Among these, PITX1, EGR1, and MYF6 were the most enriched transcription factors in fast muscle fibers (Figure 4E). In contrast, ZSCAN30 and EPAS1 (also known as HIF2A) were the most enriched transcription factors in slow muscle fibers (Figure 4E). Consistent with this, MAFA was expressed at higher levels in the UMAP region corresponding to fast muscle fibers, while EPAS1 had the opposite expression pattern (Figure 4F).
In addition to known protein-coding genes, there are numerous non-coding RNA biotypes that may be involved in the regulation of human development and disease. 51, 52 In transcriptome datasets, several non-coding RNAs exhibit fiber type specificity (Figure 5A and Supplementary Dataset 15), including LINC01405, which is highly specific for slow fibers and is reported to be decreased in muscle from patients with mitochondrial myopathy. 53 In contrast, RP11-255P5.3, corresponding to the lnc-ERCC5-5 gene (https://lncipedia.org/db/transcript/lnc-ERCC5-5:2) 54, exhibits fast fiber type specificity. Both LINC01405 (https://tinyurl.com/x5k9wj3h) and RP11-255P5.3 (https://tinyurl.com/29jmzder) exhibit skeletal muscle specificity (Supplementary Figures 9A–B) and have no known contractile genes within their 1 Mb genomic neighborhood, suggesting that they play a specialized role in regulating fiber types rather than regulating neighboring contractile genes. The slow/fast fiber type-specific expression profiles of LINC01405 and RP11-255P5.3, respectively, were confirmed using RNAscope (Figures 5B–C).
A. Non-coding RNA transcripts are significantly regulated in slow- and fast-twitch muscle fibers. B. Representative RNAscope images showing the slow- and fast-twitch fiber type specificity of LINC01405 and RP11-255P5.3, respectively. Scale bar = 50 μm. C. Quantification of myofiber type-specific non-coding RNA expression as determined by RNAscope (n = 3 biopsies from independent individuals, comparing fast and slow muscle fibers within each individual). Statistical analysis was performed using a two-tailed Student’s t-test. Box plots show the median and the first and third quartiles, with whiskers pointing to the minimum and maximum values. D. De novo microbial protein identification workflow (created with BioRender.com). E. The microbial protein LINC01405_ORF408:17441:17358 is specifically expressed in slow skeletal muscle fibers (n=5 biopsies from independent participants, comparing fast and slow muscle fibers in each participant). Statistical analysis was performed using the Limm linear model method combined with an empirical Bayesian approach, followed by the Benjamini-Hochberg method for multiple comparisons with p-value adjustment. Box plots show the median, first and third quartiles, with whiskers pointing to the maximum/minimum values.
Recently, studies have shown that many putative non-coding transcripts encode transcribed microbial proteins, some of which regulate muscle function. 44, 55 To identify microbial proteins with potential fiber type specificity, we searched our 1000 fiber proteome dataset using a custom FASTA file containing the sequences of non-coding transcripts (n = 305) found in the 1000 fiber transcriptome dataset (Figure 5D). We identified 197 microbial proteins from 22 different transcripts, 71 of which were differentially regulated between slow and fast skeletal muscle fibers (Supplementary Figure 9C and Supplementary Data Set 16). For LINC01405, three microbial protein products were identified, one of which showed similar slow fiber specificity to its transcript (Figure 5E and Supplementary Figure 9D). Thus, we identified LINC01405 as a gene encoding a microbial protein specific for slow skeletal muscle fibers.
We developed a comprehensive workflow for large-scale proteomic characterization of individual muscle fibers and identified regulators of fiber heterogeneity in healthy states. We applied this workflow to understand how nemaline myopathies affect skeletal muscle fiber heterogeneity. Nemaline myopathies are inherited muscle diseases that cause muscle weakness and, in affected children, present with a range of complications including respiratory distress, scoliosis, and limited limb mobility. 19,20 Typically, in nemaline myopathies, pathogenic variants in genes such as actin alpha 1 (ACTA1) result in a predominance of slow-twitch fiber myofiber composition, although this effect is heterogeneous. One notable exception is troponin T1 nemaline myopathy (TNNT1), which has a predominance of fast fibers. Thus, a better understanding of the heterogeneity underlying the skeletal muscle fiber dysregulation observed in nemaline myopathies may help to unravel the complex relationship between these diseases and myofiber type.
Compared with healthy controls (n=3 per group), myofibers isolated from nemaline myopathy patients with mutations in the ACTA1 and TNNT1 genes showed marked myofiber atrophy or dystrophy (Figure 6A, Supplementary Table 3). This presented significant technical challenges for proteomic analysis due to the limited amount of available material. Despite this, we were able to detect 2485 proteins in 272 skeletal myofibers. After filtering for at least 1000 quantified proteins per fiber, 250 fibers were subjected to subsequent bioinformatics analysis. After filtering, an average of 1573 ± 359 proteins per fiber were quantified (Supplementary Figure 10A, Supplementary Data Sets 17–18). Notably, despite the significant reduction in fiber size, the proteome depth of nemaline myopathy patient samples was only modestly reduced. Moreover, processing these data using our own FASTA files (including non-coding transcripts) allowed us to identify five microbial proteins in skeletal myofibers from nemaline myopathy patients (Supplementary Data Set 19). The dynamic range of the proteome was significantly broader, and total proteins in the control group correlated well with the results of a previous 1000-fiber proteome analysis (Supplementary Fig. 10B–C).
A. Microscopic images showing fiber atrophy or dystrophy and the predominance of different fiber types based on MYH in ACTA1 and TNNT1 nemaline myopathies (NM). Scale bar = 100 μm. To ensure reproducibility of staining in ACTA1 and TNNT1 patients, three patient biopsies were stained two to three times (four sections per case) before selecting representative images. B. Fiber type proportions in participants based on MYH. C. Principal component analysis (PCA) plot of skeletal muscle fibers in patients with nemaline myopathies and controls. D. Skeletal muscle fibers from patients with nemaline myopathies and controls projected onto a PCA plot determined from the 1000 fibers analyzed in Figure 2. For example, volcano plots comparing differences between participants with ACTA1 and TNNT1 nemaline myopathies and controls, and between participants with ACTA1 and TNNT1 nemaline myopathies. Colored circles denote proteins that were significantly different at π < 0.05, and dark dots denote proteins that were significantly different at FDR < 0.05. Statistical analysis was performed using the Limma linear model method and empirical Bayesian methods, followed by p-value adjustment for multiple comparisons using the Benjamini-Hochberg method. H. Enrichment analysis of significantly differentially expressed proteins across the entire proteome and in type 1 and 2A fibers. Statistical analysis was performed using the clusterProfiler package and Benjamini-Hochberg adjusted p-values. I, J. Principal component analysis (PCA) plots colored by extracellular matrix and mitochondrial gene ontology (GO) terms.
Because nemaline myopathies can influence the proportion of MYH-expressing myofiber types in skeletal muscle,19,20 we first examined the MYH-expressing myofiber types in patients with nemaline myopathies and controls. We determined myofiber type using an unbiased method previously described for the 1000 myofiber assay (Supplementary Figs. 10D–E) and again failed to identify pure 2X myofibers (Figure 6B). We observed a heterogeneous effect of nemaline myopathies on myofiber type, as two patients with ACTA1 mutations had an increased proportion of type 1 myofibers, whereas two patients with TNNT1 nemaline myopathy had a decreased proportion of type 1 myofibers (Figure 6B). Indeed, the expression of MYH2 and fast troponin isoforms (TNNC2, TNNI2, and TNNT3) was decreased in ACTA1-nemaline myopathies, whereas MYH7 expression was decreased in TNNT1-nemaline myopathies (Supplementary Figure 11A). This is consistent with previous reports of heterogeneous myofiber type switching in nemaline myopathies.19,20 We confirmed these results by immunohistochemistry and found that patients with ACTA1-nemaline myopathy had a predominance of type 1 myofibers, whereas patients with TNNT1-nemaline myopathy had the opposite pattern (Figure 6A).
At the single-fiber proteome level, skeletal muscle fibers from ACTA1 and TNNT1 nemaline myopathy patients clustered with the majority of control fibers, with TNNT1 nemaline myopathy fibers generally being the most severely affected (Figure 6C). This was particularly evident when plotting principal component analysis (PCA) plots of pseudo-inflated fibers for each patient, with TNNT1 nemaline myopathy patients 2 and 3 appearing the most distant from control samples (Supplementary Figure 11B, Supplementary Data Set 20). To better understand how fibers from myopathy patients compare to healthy fibers, we used detailed information obtained from proteomic analysis of 1,000 fibers from healthy adult participants. We projected fibers from the myopathy dataset (ACTA1 and TNNT1 nemaline myopathy patients and controls) onto the PCA plot obtained from the 1000-fiber proteomic analysis (Figure 6D). The distribution of MYH fiber types along PC2 in control fibers was similar to the fiber distribution obtained from the 1000-fiber proteomic analysis. However, most fibers in nemaline myopathy patients shifted down PC2, overlapping with healthy fast-twitch fibers, regardless of their native MYH fiber type. Thus, although patients with ACTA1 nemaline myopathy showed a shift towards type 1 fibers when quantified using MYH-based methods, both ACTA1 nemaline myopathy and TNNT1 nemaline myopathy shifted the skeletal muscle fiber proteome towards fast-twitch fibers.
We then directly compared each patient group with healthy controls and identified 256 and 552 differentially expressed proteins in ACTA1 and TNNT1 nemaline myopathies, respectively (Figure 6E–G and Supplementary Figure 11C, Supplementary Data Set 21). Gene enrichment analysis revealed a coordinated decrease in mitochondrial proteins (Figure 6H–I, Supplementary Data Set 22). Surprisingly, despite the differential predominance of fiber types in ACTA1 and TNNT1 nemaline myopathies, this decrease was completely independent of the MYH-based fiber type (Figure 6H and Supplementary Figures 11D–I, Supplementary Data Set 23). Three microbial proteins were also regulated in ACTA1 or TNNT1 nemaline myopathies. Two of these microproteins, ENSG00000215483_TR14_ORF67 (also known as LINC00598 or Lnc-FOXO1) and ENSG00000229425_TR25_ORF40 (lnc-NRIP1-2), showed differential abundance only in type 1 myofibers. ENSG00000215483_TR14_ORF67 has previously been reported to play a role in cell cycle regulation. 56 On the other hand, ENSG00000232046_TR1_ORF437 (corresponding to LINC01798) was increased in both type 1 and type 2A myofibers in ACTA1-nemaline myopathy compared with healthy controls (Supplementary Figure 12A, Supplementary Data Set 24). In contrast, ribosomal proteins were largely unaffected by nemaline myopathy, although RPS17 was downregulated in ACTA1 nemaline myopathy (Fig. 6E).
Enrichment analysis also revealed upregulation of immune system processes in ACTA1 and TNNT1 nemaline myopathies, while cell adhesion was also increased in TNNT1 nemaline myopathy (Figure 6H). The enrichment of these extracellular factors was reflected by the extracellular matrix proteins shifting PCA in PC1 and PC2 in a negative direction (i.e., toward the most affected fibers) (Figure 6J). Both patient groups showed increased expression of extracellular proteins involved in immune responses and sarcolemmal repair mechanisms, such as annexins (ANXA1, ANXA2, ANXA5)57,58 and their interacting protein S100A1159 (Supplementary Figures 12B–C). This process has been previously reported to be enhanced in muscular dystrophies60 but, to our knowledge, has not previously been associated with nemaline myopathies. Normal function of this molecular machinery is required for sarcolemmal repair following injury and for the fusion of newly formed myocytes with myofibers58,61. Thus, the increased activity of this process in both patient groups suggests a reparative response to injury caused by myofiber instability.
The effects of each nemaline myopathy were well correlated (r = 0.736) and showed reasonable overlap (Supplementary Figures 11A–B), indicating that ACTA1 and TNNT1 nemaline myopathy have similar effects on the proteome. However, some proteins were regulated only in ACTA1 or TNNT1 nemaline myopathy (Supplementary Figures 11A and C). The profibrotic protein MFAP4 was one of the most upregulated proteins in TNNT1 nemaline myopathy but remained unchanged in ACTA1 nemaline myopathy. SKIC8, a component of the PAF1C complex responsible for regulating HOX gene transcription, was downregulated in TNNT1 nemaline myopathy but not affected in ACTA1 nemaline myopathy (Supplementary Figure 11A). Direct comparison of ACTA1 and TNNT1 nemaline myopathy revealed greater reductions in mitochondrial proteins and increases in immune system proteins in TNNT1 nemaline myopathy (Figure 6G–H and Supplementary Figures 11C and 11H–I). These data are consistent with the greater atrophy/dystrophy observed in TNNT1 nemaline myopathy compared to TNNT1 nemaline myopathy (Figure 6A), suggesting that TNNT1 nemaline myopathy represents a more severe form of the disease.
To assess whether the observed effects of nemaline myopathy persist at the whole muscle level, we performed a bulk proteomic analysis of muscle biopsies from the same cohort of TNNT1 nemaline myopathy patients and compared them with controls (n=3 per group) (Supplementary Fig. 13A, Supplementary Data Set 25). As expected, controls were closely related in principal component analysis, whereas TNNT1 nemaline myopathy patients exhibited higher intersample variability similar to that seen in single fiber analysis (Supplementary Fig. 13B). Bulk analysis reproduced the differentially expressed proteins (Supplementary Fig. 13C, Supplementary Data Set 26) and biological processes (Supplementary Fig. 13D, Supplementary Data Set 27) highlighted by comparing individual fibers, but lost the ability to distinguish between different fiber types and failed to account for heterogeneous disease effects across fibers.
Taken together, these data demonstrate that single-myofiber proteomics can elucidate clinical biological features that are not detectable by targeted methods such as immunoblotting. Moreover, these data highlight the limitations of using actin fiber typing (MYH) alone to describe phenotypic adaptation. Indeed, although fiber type switching differs between actin and troponin nemaline myopathies, both nemaline myopathies decouple MYH fiber typing from skeletal muscle fiber metabolism toward a faster and less oxidative muscle proteome.
Cellular heterogeneity is critical for tissues to meet their diverse demands. In skeletal muscle, this is often described as fiber types characterized by different degrees of force production and fatigability. However, it is clear that this explains only a small part of the skeletal muscle fiber variability, which is much more variable, complex and multifaceted than previously thought. Technological advances have now shed light on the factors regulating skeletal muscle fibers. Indeed, our data suggest that type 2X fibers may not be a distinct skeletal muscle fiber subtype. Moreover, we identified metabolic proteins, ribosomal proteins and cell-associated proteins as major determinants of skeletal muscle fiber heterogeneity. By applying our proteomic workflow to patient samples with nematode myopathy, we further demonstrated that MYH-based fiber typing does not fully reflect skeletal muscle heterogeneity, especially when the system is perturbed. Indeed, regardless of the MYH-based fiber type, nematode myopathy results in a shift towards faster and less oxidative fibers.
Skeletal muscle fibers have been classified since the 19th century. Recent omics analyses have allowed us to begin to understand the expression profiles of different MYH fiber types and their responses to different stimuli. As described here, omics approaches also have the advantage of greater sensitivity for quantifying fiber type markers than traditional antibody-based methods, without relying on the quantification of a single (or a few) markers to define a skeletal muscle fiber type. We used complementary transcriptomic and proteomic workflows and integrated the results to examine the transcriptional and post-transcriptional regulation of fiber heterogeneity in human skeletal muscle fibers. This workflow resulted in the failure to identify pure 2X-type fibers at the protein level in the vastus lateralis of our cohort of healthy young men. This is consistent with previous single fiber studies that found <1% pure 2X fibers in healthy vastus lateralis, although this should be confirmed in other muscles in the future. The discrepancy between the detection of nearly pure 2X fibers at the mRNA level and only mixed 2A/2X fibers at the protein level is puzzling. MYH isoform mRNA expression is not circadian,67 suggesting that we are unlikely to have “missed” the MYH2 onset signal in seemingly pure 2X fibers at the RNA level. One possible explanation, although purely hypothetical, could be differences in protein and/or mRNA stability between MYH isoforms. Indeed, no fast fiber is 100% pure for any MYH isoform, and it is unclear whether MYH1 mRNA expression levels in the 70–90% range would result in equal MYH1 and MYH2 abundance at the protein level. However, when considering the entire transcriptome or proteome, cluster analysis can confidently identify only two distinct clusters representing slow and fast skeletal muscle fibers, regardless of their precise MYH composition. This is consistent with analyses using single-nucleus transcriptomic approaches, which typically identify only two distinct myonuclear clusters. 68, 69, 70 Furthermore, although previous proteomic studies have identified type 2X fibers, these fibers do not cluster separately from the rest of the fast fibers and show only a small number of differentially abundant proteins compared to other fiber types based on MYH. 14 These results suggest that we should return to the early 20th century view of muscle fiber classification, which divided human skeletal muscle fibers not into three distinct classes based on MYH, but into two clusters based on their metabolic and contractile properties. 63
More importantly, myofiber heterogeneity should be considered along multiple dimensions. Previous “omics” studies have pointed in this direction, suggesting that skeletal muscle fibers do not form discrete clusters but are arranged along a continuum. 11, 13, 14, 64, 71 Here, we show that, in addition to differences in contractile and metabolic properties of skeletal muscle, myofibers can be differentiated by features related to cell-cell interactions and translation mechanisms. Indeed, we found ribosome heterogeneity in skeletal muscle fibers that contributes to heterogeneity independent of slow and fast fiber types. The underlying cause of this substantial myofiber heterogeneity, independent of slow and fast fiber type, remains unclear, but it may point to specialized spatial organization within muscle fascicles that optimally respond to specific forces and loads,72 specialized cellular or organ-specific communication with other cell types in the muscle microenvironment73,74,75 or differences in ribosome activity within individual myofibers. Indeed, ribosomal heteroplasmy, either through paralogous substitution of RPL3 and RPL3L or at the level of 2′O-methylation of rRNA, has been shown to be associated with skeletal muscle hypertrophy76,77. Multi-omic and spatial applications combined with functional characterization of individual myofibers will further advance our understanding of muscle biology at the multi-omic level78.
By analyzing the proteomes of single myofibers from patients with nemaline myopathies, we also demonstrated the utility, effectiveness, and applicability of single myofiber proteomics to elucidate the clinical pathophysiology of skeletal muscle. Moreover, by comparing our workflow to global proteomic analysis, we were able to demonstrate that single myofiber proteomics yields the same depth of information as global tissue proteomics and extends this depth by accounting for interfiber heterogeneity and myofiber type. In addition to the expected (albeit variable) differences in fiber type ratio observed in ACTA1 and TNNT1 nemaline myopathies compared to healthy controls,19 we also observed oxidative and extracellular remodeling independent of MYH-mediated fiber type switching. Fibrosis has previously been reported in TNNT1 nemaline myopathies.19 However, our analysis builds on this finding by also revealing increased levels of extracellular secreted stress-related proteins, such as annexins, involved in sarcolemmal repair mechanisms, in myofibers from patients with ACTA1 and TNNT1 nemaline myopathies.57,58,59 In conclusion, increased annexin levels in myofibers from patients with nemaline myopathy may represent a cellular response to repair severely atrophic myofibers.
Although this study represents the largest single-fiber whole-muscle-omics analysis of humans to date, it is not without limitations. We isolated skeletal muscle fibers from a relatively small and homogeneous sample of participants and a single muscle (the vastus lateralis). Therefore, it is impossible to exclude the existence of specific fiber populations across muscle types and at extremes of muscle physiology. For example, we cannot exclude the possibility of a subset of ultrafast fibers (e.g., pure 2X fibers) emerging in highly trained sprinters and/or strength athletes79 or during periods of muscular inactivity66,80. Furthermore, the limited sample size of participants prevented us from investigating sex differences in fiber heterogeneity, as fiber type ratios are known to differ between men and women. Furthermore, we were unable to conduct transcriptomic and proteomic analyses on the same muscle fibers or samples from the same participants. As we and others continue to optimize single-cell and single-myofiber analyses using omics analysis to achieve ultra-low sample input (as demonstrated here in the analysis of fibers from patients with mitochondrial myopathy), the opportunity to combine multi-omics (and functional) approaches within single muscle fibers becomes apparent.
Overall, our data identify and explain transcriptional and post-transcriptional drivers of skeletal muscle heterogeneity. Specifically, we present data that challenge a long-standing dogma in skeletal muscle physiology associated with the classical MYH-based definition of fiber types. We hope to renew the debate and ultimately rethink our understanding of skeletal muscle fiber classification and heterogeneity.
Fourteen Caucasian participants (12 men and 2 women) voluntarily agreed to take part in this study. The study was approved by the Ethics Committee of Ghent University Hospital (BC-10237), complied with the 2013 Helsinki Declaration, and was registered at ClinicalTrials.gov (NCT05131555). General characteristics of the participants are presented in Supplementary Table 1. After obtaining verbal and written informed consent, participants underwent a medical examination before final inclusion in the study. Participants were young (22–42 years), healthy (no medical conditions, no smoking history), and moderately physically active. Maximal oxygen uptake was determined using a step ergometer for assessing physical fitness as described previously. 81
Muscle biopsy samples were collected at rest and in the fasting state three times, 14 days apart. Because these samples were collected as part of a larger study, participants consumed placebo (lactose), an H1-receptor antagonist (540 mg fexofenadine), or an H2-receptor antagonist (40 mg famotidine) 40 minutes before the biopsy. We have previously demonstrated that these histamine receptor antagonists do not affect resting skeletal muscle fitness81, and no state-related clustering was observed in our quality control plots (Supplementary Figures 3 and 6). A standardized diet (41.4 kcal/kg body weight, 5.1 g/kg body weight carbohydrate, 1.4 g/kg body weight protein, and 1.6 g/kg body weight fat) was maintained for 48 hours prior to each experimental day, and a standardized breakfast (1.5 g/kg body weight carbohydrate) was consumed in the morning of the experimental day. Under local anesthesia (0.5 ml 1% lidocaine without epinephrine), muscle biopsies were obtained from the vastus lateralis muscle using percutaneous Bergström aspiration.82 Muscle samples were immediately embedded in RNAlater and stored at 4°C until manual fiber dissection (up to 3 days).
Freshly isolated myofiber bundles were transferred to fresh RNAlater medium in a culture dish. Individual myofibers were then manually dissected using a stereomicroscope and fine tweezers. Twenty-five fibers were dissected from each biopsy, paying special attention to selecting fibers from different areas of the biopsy. After dissection, each fiber was gently immersed in 3 μl of lysis buffer (SingleShot Cell Lysis Kit, Bio-Rad) containing proteinase K and DNase enzymes to remove unwanted proteins and DNA. Cell lysis and protein/DNA removal were then initiated by brief vortexing, spinning the liquid in a microcentrifuge, and incubation at room temperature (10 min). The lysate was then incubated in a thermal cycler (T100, Bio-Rad) at 37°C for 5 min, 75°C for 5 min, and then immediately stored at -80°C until further processing.
Illumina-compatible polyadenylated RNA libraries were prepared from 2 µl of myofiber lysate using the QuantSeq-Pool 3′ mRNA-Seq Library Prep Kit (Lexogen). Detailed methods can be found in the manufacturer’s manual. The process begins with first-strand cDNA synthesis by reverse transcription, during which unique molecular identifiers (UMIs) and sample-specific i1 barcodes are introduced to ensure pooling of samples and reduce technical variability during downstream processing. cDNA from 96 myofibers is then pooled and purified with magnetic beads, after which RNA is removed and second-strand synthesis is performed using random primers. The library is purified with magnetic beads, pool-specific i5/i7 tags are added, and PCR amplified. A final purification step produces Illumina-compatible libraries. The quality of each library pool was assessed using the High Sensitivity Small Fragment DNA Analysis Kit (Agilent Technologies, DNF-477-0500).
Based on Qubit quantification, pools were further pooled at equimolar concentrations (2 nM). The resulting pool was then sequenced on a NovaSeq 6000 instrument in standard mode using the NovaSeq S2 Reagent Kit (1 × 100 nucleotides) with 2 nM loading (4% PhiX).
Our pipeline is based on Lexogen’s QuantSeq Pool data analysis pipeline (https://github.com/Lexogen-Tools/quantseqpool_analysis). Data were first demultiplexed with bcl2fastq2 (v2.20.0) based on the i7/i5 index. Read 2 was then demultiplexed with idemux (v0.1.6) based on the i1 sample barcode and UMI sequences were extracted with umi_tools (v1.0.1). Reads were then trimmed with cutadapt (v3.4) in multiple rounds to remove short reads (<20 in length) or reads consisting solely of adapter sequences. Reads were then aligned to the human genome using STAR (v2.6.0c) and BAM files were indexed with SAMtools (v1.11). Duplicate reads were removed using umi_tools (v1.0.1). Finally, alignment counting was performed using featureCounts in Subread (v2.0.3). Quality control was performed using FastQC (v0.11.9) at several intermediate stages of the pipeline.
All further bioinformatics processing and visualization were performed in R (v4.2.3), primarily using the Seurat (v4.4.0) workflow. 83 Therefore, individual UMI values and metadata matrices were transformed into Seurat objects. Genes expressed in less than 30% of all fibers were removed. Low-quality samples were removed based on a minimum threshold of 1000 UMI values and 1000 detected genes. Ultimately, 925 fibers passed all quality control filtering steps. UMI values were normalized using the Seurat SCTransform v2 method, 84 including all 7418 detected features, and between-participant differences were regressed out. All relevant metadata can be found in Supplementary Dataset 28.
Post time: Sep-10-2025