Bridging the muscle genome to phenome across multiple biological scales

Skeletal muscle is a tissue with highly hierarchical organization across multiple scales of biological organization. Muscle tissue also has an enigmatic capacity to meet highly diverse functional and physiological demands, ranging from slow, submaximal and prolonged postural contractions to fast and powerful maximal contractions. At the finer organizational scales of molecules, genes and proteins, the modulatory capacity of muscle is thought to correspond with the ensemble of motor proteins. The isoform makeup of these proteins can be highly heterogeneous and this heterogeneity shapes muscle contractile function at the coarser scales of sarcomeres, fibers, tissues and the organ. At these coarser scales, muscle physiological function is also influenced by sarcomere architecture, and lengths during rest and operation. Hence, to define and relate the molecular- and protein-scale underpinnings of adaptations and reconcile them across muscles with different and specialized functions, an interdisciplinary multi-scale approach is needed. With few exceptions (Horak et al., 2016; Narici et al., 2016), such multi-scale analyses are lacking.

The precisely defined sarcomeric structure, combined with muscle's deeply hierarchical architecture, influences the functional scope of muscle fibers such that the force production of skeletal muscle is highly length – and velocity – dependent (Lieber and Ward, 2011). Muscle contraction and force generation result from interactions between the thick (myosin-containing) and thin (actin-containing) filaments sliding past each other within sarcomeres (Huxley and Niedergerke, 1954; Huxley and Hanson, 1954). Thick and thin filament composition also influences muscle contractility. In particular, the myosin heavy chain (MHC) isoform composition within sarcomeres has been shown to define the speed of cross-bridge cycling and contractile efficiency, and also to influence the energy metabolism of different fiber types (Burke et al., 1971; Hettige et al., 2020; Pette and Staron, 2000; Schiaffino and Reggiani, 2011; Wells, 1965). This isoform dependency of converting chemical energy into mechanical energy by myosin has been used to establish a paradigm of characterizing muscle fibers as primarily fast twitch (MHC IIa, MHC IIx, MHC IIb containing) and slow twitch (MHC I containing) (Pellegrino et al., 2003). Varying populations of fast or slow twitch fibers can be identified in different muscles, from embryonic development to adult stages (Agbulut et al., 2003; Pette and Staron, 2000; Stål et al., 1994; Staron and Johnson, 1993), and thus the relative proportions of fast and slow muscle fibers provide an alternative avenue for muscles to vary contractile capabilities required to meet functional demands.

Several techniques have been used to study specific structural and functional attributes of muscle types, but such studies have mostly focused on single scales of biological organization rather than using a multi-scale approach. At the whole-muscle scale, muscle ergometry studies have established differences in physiological properties between fast- and slow-twitch fibers (Close, 1964; Hill, 1938; Huxley, 1957) on the basis of force–length (FL), and force–velocity–power (FVP) relationships (Burkholder et al., 1994; Lieber and Fridén, 2000). At the protein scale, SDS-PAGE gel electrophoresis, immunohistochemistry and in vitro motility assays have been used to reveal the relative populations, ATPase rate and shortening velocities of slow- or fast-twitch fiber types (Höök et al., 1999; Pette and Staron, 2000; Schiaffino and Reggiani, 1994). Recent advances in RNA sequencing and other bioinformatics approaches have permitted exploration of the transcriptomic profile of muscles, allowing the identification of genetic and expression profiles of individual muscle types, but often lacking interrogation of functional relevance (Hettige et al., 2020; Terry et al., 2018). Therefore, multi-scale approaches that integrate knowledge about muscle contractile properties and protein isoform composition are necessary to link the gene, protein and whole-muscle scales, to facilitate the integration of active and passive sarcomeric components, and to define the aggregate behavior of heterogeneous mixtures of myosins.

Here, we took an interdisciplinary approach to identify and characterize the key differences underlying the distinct functional properties of mouse soleus and tibialis anterior (TA) at the genome, protein and physiology scales. These muscles have distinct functional characteristics: mouse soleus is a fairly slow ankle plantar flexor (Asmussen and Maréchal, 1989; Holt and Askew, 2012) characterized by slow myosin isoform expression and relatively slow contraction velocities. It also has a fairly parallel-fibered architecture (Burkholder et al., 1994) with reported pennation angles ranging from 5 to 11.4 deg (Brooks and Faulkner, 1988; Charles et al., 2016; Lal et al., 2021) and a slow-twitch, oxidative fiber composition (Asmussen and Maréchal, 1989; Burkholder et al., 1994), with a low energy cost (Westerblad et al., 2010). Mouse TA is a predominately fast-twitch ankle dorsiflexor (Burkholder and Lieber, 1998) characterized by fast myosin isoform expression, relatively fast contraction velocities and a bipennate architecture (Burkholder et al., 1994; Charles et al., 2016; Moo et al., 2016), with reported pennation angles ranging from 11.7 to 24 deg and greater pennation angles being determined using non-destructive diffusion tensor imaging (Burkholder et al., 1994; Charles et al., 2016; Heemskerk et al., 2005; Lovering et al., 2013). Correspondingly, TA muscles are composed of mainly fast-twitch, glycolytic fibers (Burkholder et al., 1994) with greater energy cost (Boss et al., 2018). A comparison of these structurally and functionally different muscles is useful for determining how differences in gene expression patterns (measured using Tag-Seq) influence protein expression, whilst not necessarily directly correlating to protein isoform distribution (measured via SDS-PAGE) and unloaded translation speed for actin (measured via in vitro motility assay) and whole muscle (measured via ergometry).

We hypothesized that while myosin isoform expression strongly influences contraction force and speed (Pellegrino et al., 2003), contraction velocity is further refined at the muscle scale by additional, potentially contrasting and/or reinforcing, aspects of muscle fiber architecture (Bottinelli et al., 1991; Close, 1965; Lieber and Ward, 2011; Winters et al., 2011) that ultimately fine tune force–velocity curvature and the optimal force (F_opt) and velocity (V_opt) for peak power (Lutz et al., 2001). Because of the established correlation between specific myosin and titin isoform expression (Eckels et al., 2018; Freiburg et al., 2000; Linke et al., 1996; Prado et al., 2005; Trombitás et al., 1998), we further tested the hypothesis that active and passive muscle properties are tuned in combination to establish muscle functional phenotype. Our multi-scale approach indeed demonstrates that muscle shortening velocity is not a direct extrapolation from molecular actin sliding velocity, highlighting the importance of considering the influences of architectural, active and passive muscle properties when assessing the causes of whole-muscle function. Thus, contrary to general assumptions, overall muscle performance is imperfectly reflected by gene-, protein- and fiber-scale characteristics, and is instead shaped by intricate fine-tuning across scales into the specialized function of individual muscles.

All experiments were performed in accordance with approved IACUC protocols at both University of Massachusetts, Lowell and Harvard University. Wild-type mice (mean±s.d. body mass 32.1±7.4 g; total n=20; male n=9, female n=11) were obtained from an out-bred colony from no known specific strain and individually housed with ad libitum food and water access. Detailed procedures conducted at each organizational scale are described below.

Collected muscle samples were stored in RNAlater™ stabilization solution (Thermo Fisher Scientific) at −80°C until RNA extraction. Total RNA was extracted from five TA and soleus muscles (male n=3 and female n=2) using a Qiagen fibrous tissue total RNA extraction mini kit. The concentration and quality of RNA were measured using the Qubit RNA Broad-Range assay and Agilent 2100 Bioanalyzer RNA 6000 Nano assay. A 10 ng sample of total RNA from the extracted samples was used for Tag-Seq library preparation. The RNA integrity values (RIN) for these samples ranged from 5.4 to 7.8, with a mean (±s.d.) for all the samples of 7.06±0.96, which meets the recommended criteria in the Takara SMART-Seq Stranded protocol for Illumina sequencing platforms. The cDNA libraries were quantified using a Qubit dsDNA HS kit (Thermo Fisher Scientific) and library sizes were quantified using an Agilent bioanalyzer 2100 High Sensitivity dsDNA assay. An extra quantification step was performed using qPCR with KAPA Library Quantification qPCR Kit for Illumina sequencing platforms. The cDNA libraries with concentration >3 ng µl⁻¹ and average fragment size of 300–450 bp were sequenced with an Illumina NextSeq 500 sequencer using the sequencing guidelines suggested by Illumina.

Quality of the short-read fragments generated from sequencing was assessed using the FastQC quality control tool for high throughput sequence data. Fastq files with a per-base quality score below 20 were preprocessed using the sliding window quality filtering (window size 4) in Trimmomatic v0.32. Preprocessed short reads were aligned to Mus musculus GRCm38.p4 reference genome annotation using the alignment pipeline available on GitHub and archived in Eli-Meyers (http://eli-meyer.github.io/TagSeq_utilities/index.html?)

Count files generated from the read alignments were obtained and used for downstream data analysis in R v3.6.1. Prior to the analysis, the genes with a normalized read count <3 were removed to increase the sensitivity of differential gene expression. The differential gene expression analysis was performed between soleus and TA tissues using R library DESeq2. Significantly differentially expressed genes were identified at a 95% significance level using Benjamini–Hochberg corrected P-values.

Variance stabilizing transformation was applied to the read counts to improve the homogeneity of variance and to normalize the counts, followed by removing batch effects and noise in the dataset prior to data visualization.

Principal component analysis (PCA) was performed on selected differentially expressed genes. Differentially expressed genes were ranked based on the loading values of the first principal component and genes with the highest contribution to the variance were selected using the cutoff (0.07) defined using the broken stick model (Jackson, 1993). Selected genes were visualized in an expression heatmap and the genes were grouped into four clusters (I–IV) based on their expression profiles using the Ward.D2 clustering method.

PSI values for each exon coding for the titin PEVK region were determined using the approach developed by Schafer et al. (2015). Mus musculus GRCm38.100 genome annotation was used to extract titin exon coordinates and utilized as the reference exon annotation. Titin expression data for soleus (N=6) and TA (N=6) were obtained from the GEO database (accession number: GSE139204; 10 month old wild-type mice). Gathered RNA sequencing data were aligned to Mus musculus GRCm38.100 genome using the STAR alignment tool and used for the subsequent PSI calculations. Calculated PSI values range from 0 to 1, with 1 being constitutive expression of the respective exon, values <1 indicating splicing events and 0 denoting complete exclusion of the exon from the expressed transcript pool.

Myosin for the in vitro motility assay was extracted from chicken pectoralis, mouse cardiac, TA and soleus muscles as previously described (Bárány, 1967; Malmqvist et al., 2004), but with modifications. Skeletal muscle tissue (∼100–150 mg) was finely chopped and homogenized on ice after rinsing in 1× PBS buffer. Extraction buffer containing 300 mmol l⁻¹ KCl, 150 mmol l⁻¹ imidazole pH 6.8, 10 mmol l⁻¹ Na₄P₂O₇ pH 6.8, 1 mmol l⁻¹ MgCl₂, 2 mmol l⁻¹ DTT and 1:100 ratio of protease inhibitor cocktail (MedChemExpress) was added to the muscle mince and stirred for 20 min at 4°C. After centrifugation at 200,000 g for 1.5 h at 4°C, the supernatant was diluted 20-fold in 2 mmol l⁻¹ DTT. The protein precipitated overnight at 4°C and was then centrifuged at 30,000 g for 30 min at 4°C. The concentrated pellet was dissolved in 600 mmol l⁻¹ KCl, 2 mmol l⁻¹ MgCl₂, 5 mmol l⁻¹ DTT and 10 mmol l⁻¹ Tris-HCl pH 7.6. Myosin was stored in 50% glycerol at −20°C and myosin concentration was determined using spectrophotometer readings at 280 nm, corrected at 320 nm. Porcine cardiac actin was purified from acetone powder as previously described (Spudich and Watt, 1971).

To determine MHC isoform composition, purified chicken pectoralis, mouse cardiac, soleus and TA samples from above were prepared for SDS-PAGE by dilution with 5× sample buffer [1.25 mol l⁻¹ Tris pH 6.8, 48% (w/v) urea, 50% (w/v) sucrose, 5% (w/v) SDS, 0.5% (w/v) Bromophenol Blue] to a final concentration of 0.13 mg ml⁻¹ for chicken pectoralis, 0.16 mg ml⁻¹ for cardiac muscle, 1.13 mg ml⁻¹ for soleus and 0.88 mg ml⁻¹ for TA. Myosin isoforms in these samples were separated as previously described (Blough et al., 1996) using a Mini-PROTEAN Tetra Cell (Bio-Rad Laboratories, Hercules, CA, USA), with modifications. The stacking gel was 4% acrylamide:bis-acrylamide (with 30% glycerol) and the separating gel contained 8% acrylamide:bis-acrylamide [with 35% (v/v) glycerol]. Both upper and lower gels were degassed prior to adding ammonium persulfate and TEMED, and the upper gel chamber was leak-proofed with a gel of ∼1% agarose in running buffer containing 0.025 mol l⁻¹ Tris, 0.192 mol l⁻¹ glycine and 0.1% SDS. Gels were run for ∼24–26 h at 4°C, with additional 1% betamercaptoethanol in the upper buffer chamber. Gels were stained with Coomassie Blue dye and imaged using a Bio-Rad ChemiDoc XRS Imaging System (Bio-Rad Laboratories, CA, USA). The densitometric analysis for determining band intensity in the gel images was done using ImageJ as described previously (Sundar et al., 2020). Values were averaged across three biological replicates and tests of statistical significance between myosin isoform pairs in the soleus and TA were performed with GraphPad Prism 9 software, as an unpaired t-test where P-values less than 0.05 were considered significant.

An unregulated in vitro motility assay was performed as per Orzechowski et al. (2014) with modifications. Non-functional myosin was pelleted by diluting myosin with 1.1 μmol l⁻¹ actin and 1 mmol l⁻¹ ATP in myosin buffer containing 300 mmol l⁻¹ KCl, 25 mmol l⁻¹ imidazole, 1 mmol l⁻¹ EGTA, 4 mmol l⁻¹ MgCl₂ and 10 mmol l⁻¹ DTT, pH 7.4. After centrifugation at 100,000 g for 30 min at 4°C, the myosin concentration in the supernatant was determined by Bradford assay (Bio-Rad Laboratories) and diluted to 100 μg ml⁻¹ in myosin buffer. For slide preparation, a nitrocellulose-coated cover-slip was placed on a standard glass slide with double-sided tape, thereby forming a channel between the glass slide and the nitrocellulose surface. One volume (25 μl) of diluted myosin was loaded into the chamber and incubated for 1 min, followed by a wash with 1 mg ml⁻¹ bovine serum albumin (BSA) diluted in the myosin buffer to avoid any non-specific binding. Two volumes of actin buffer containing 55 mmol l⁻¹ KCl, 5 mmol l⁻¹ EGTA, 4 mmol l⁻¹ MgCl₂ and 10 mmol l⁻¹ DTT, pH 7.4 was flowed through the chamber before blocking any non-functional myosin binding with 1 μmol l⁻¹ actin diluted in actin buffer for 2 min. This was followed by two sequential washes of actin buffer containing 1 mmol l⁻¹ ATP, and then four washes with actin buffer, to remove any unbound actin. Actin labeled with TRITC was then added and allowed interaction with the myosin surface for 1 min. One volume of motility buffer containing 10 mmol l⁻¹ calcium (pCa 4.0), 2 mmol l⁻¹ dextrose, 160 units of glucose oxidase, 2 μmol l⁻¹ catalase and several concentrations of ATP (0–1 mmol l⁻¹) was added to initiate motility of actin filaments.

Data from three technical replicates were analyzed using the MTrackJ plugin in ImageJ (https://imagescience.org/meijering/software/mtrackj/). An average of 12–15 filaments were analyzed from each of the five movies recorded for each experiment to obtain the sliding velocity of actin under each condition. Data were then averaged to obtain mean actin sliding velocity with a standard deviation (s.d.) measurement and plotted in GraphPad Prism 8.0.2 where the Michaelis–Menten equation was fitted as Y=(V_max×x)/(K_m+x) (1); where V_max is the maximal velocity of actin translation across myosin and K_m is the substrate concentration at half-maximal velocity. Statistical significance of V_max values between the soleus and TA was determined by a two-tailed unpaired t-test with Welsh's correction and P<0.05 was considered significant.

Muscle dissection and isolation was performed as previously detailed (Askew et al., 1997; Askew and Marsh, 1997; Hakim et al., 2013). Briefly, mice were deeply anesthetized by administering isoflurane via a nose cone (3% induction, 0.8% maintenance). Skin and surrounding connective tissues were carefully removed from one hindlimb. Room temperature (∼23°C) Ringer's solution of pH 7.4 (composition in mmol l⁻¹: 144 NaCl, 6 KCl, 2 CaCl₂, 1 MgSO₄, 1 NaH₂PO₄, 10 glucose, 10 Hepes; Daut and Elzinga, 1989) was frequently applied to exposed tissue via syringe. The distal tendons of the lateral and medial gastrocnemius were carefully dissected away from their attachment to the calcar and reflected. Similarly, the distal tendon of the plantaris was cut and the muscle was carefully removed to reveal the soleus in the posterior compartment. The proximal and distal tendon of the soleus were then tied with braided 4/0 silk suture (Ethicon) at each myotendinous junction. Resting muscle length (L_mus) was measured from knot to knot with the ankle held at 90 deg (Rassier et al., 1999). Both tendons were then cut and the muscle was immediately placed in a bath of circulated oxygenated Ringer's solution (pH 7.4 and 25°C, monitored via a thermometer). The suture at the distal end of the soleus was tied to a stopper and the proximal end to the lever arm of a 1 N ergometer (Aurora Scientific Inc., model 300C-LR, Aurora, ON, Canada) (Fig. 1A). After mounting, the preparation was set to L_mus. In vitro testing of the soleus commenced while the experimental limb was further dissected to isolate the TA and secure a custom nerve-cuff around the sciatic nerve prior to mounting of the TA (Fig. 1B). The distal tendon of the TA was tied to the lever arm of a 5 N ergometer (Aurora Scientific Inc., model 305C) using braided 3/0 silk suture (Ethicon). The patella tendon was secured to a metal rod on the stereotaxic stage that held the motor and the subject, to minimize in-series compliance of the preparation. The TA was covered with Ringer’s solution-soaked gauze and kept at a near-constant temperature of 25°C using an adjustable infrared lamp and a thermocouple. Both preparations were set to L_mus and allowed to rest for 10 min following rigging.

To elicit contractions, electrical stimuli were transmitted via parallel platinum paddle electrodes (soleus) or a custom-made bipolar nerve cuff (TA), respectively, using a Grass stimulator (S48, Grass Technologies, West Warwick, RI, USA). Stimulation parameters were optimized by muscle to the following literature values: stimulation duration (soleus and TA: 300 ms), pulse width (soleus: 1.0 ms; TA: 0.2 ms), stimulation frequency (soleus and TA: 150 Hz) (Askew and Marsh, 1997; Hakim et al., 2013; Holt and Askew, 2012; James et al., 1995; Moran et al., 2005). For each experiment, stimulation was optimized to elicit supramaximal contractions (soleus and TA: 1–4 V). A 3 min rest period was held between contractions. Changes of force, length and stimulation voltage were recorded at 4 kHz via a NIDAQ A/D converter into IGOR pro (v.8.0; Wavemetrics, Lake Oswego, OR, USA).

FL curves (Fig. 1E) were constructed for twitch and tetanic conditions to determine muscle optimal length (L_o). Starting from a length shorter than slack length (L_slack), muscles were subjected to isometric contractions at increasing lengths. Fit lines were constructed for active and passive force using third-order polynomial fits in IGOR pro and maximum active force (P_o) and L_o were determined from these data fits. After normalizing active FL curves using tetanic L_o, we measured the tetanic and twitch plateau lengths. Tetanic L_o was used to determine starting length for remaining tests and P_o was used to monitor preparation viability at regular intervals; experiments in which P_o decreased by more than 10% over the course of the experiment were excluded from further analysis. The width of the ascending limb of the active FL curve (henceforth called the plateau width) was measured as the range of strain from 90% to 100% P_o on the ascending FL curve limb. Slack length was measured for each preparation where passive force departed from zero. As an indication of relative differences in passive stiffness of the two muscles, we measured the force borne by the preparation at 110% L_o and the passive FL curve slope between L_o and 110% L_o.

Passive muscles were stretched to slightly above L_o and subjected to isotonic contractions against loads that were successively decreased between trials (Fig. 1F). Velocity was calculated as the derivative of length and power was calculated as force×velocity. Using a custom implementation of the Hill equation (P+a)×(v+b)=b×(P_o+a) in IGOR pro (Tijs et al., 2021), we then obtained Hill-fits for FV and PV data, and normalized FVP data to muscle mass-specific force (N g⁻¹) and power (W kg⁻¹), optimal length (L_o s⁻¹) and maximum shortening velocity (%V_max). F_opt and V_opt were measured as the force and velocity, respectively, at which peak power was generated by the preparation. Force–velocity curvature was measured as a/P_o (Wilkie, 1949). After testing, muscles were carefully excised, free tendon was removed, and the belly was systematically blotted on Kim-wipes to remove moisture and weighed to the nearest 0.001 g.

General linear model (GLM) designs (Systat 11) were used to test for statistical differences between muscles (mean±s.d.) in each of the following variables: P_o, L_o, plateau shift, plateau width, slack length, passive stiffness, peak power, V_max, F_opt and V_opt. For all tests, significance level (α) was set at 0.01. We disaggregated sex and individual in our GLM design and simplified the model when these factors were not statistically significant. Pooled samples were then analytically re-assessed between muscle groups.

We studied gene expression profiles from Tag-Seq to determine differences in transcriptomic signatures between the soleus and TA. Analysis of differential gene expression revealed 1061 genes with significant differences in expression between soleus and TA samples. PCA revealed a separation of slow- and fast-twitch tissue along PC1, accounting for 73% of the variation in the dataset. Meanwhile, there was no differentiation of muscles along PC2 (Fig. 2A). Given the lack of distinct sex-specific patterns, sex was not considered in subsequent gene analyses (Fig. 2A). Loading values along PC1 were used to rank genes and identify the strongest transcriptomic features that distinguished soleus and TA. The broken stick approach (Jackson, 1993) revealed 28 genes (descriptions provided in Table S1) above the cutoff value of 0.07 for the gene loading values of PC1 (Fig. 2B). An expression heatmap of genes above the selected cutoff showed four distinct gene clusters, which were designated as gene clusters I–IV (Fig. 2C). As distinct expression pattern differences were observed within the clusters, and between the muscle types (Fig. 2C), each cluster was analyzed to determine potential gene ontologies that correlated with the associated genes.

Genes that code for MHC isoforms IIb (Myh4) and I (Myh7), which are typically expressed in fast-twitch and slow-twitch muscles, respectively, were identified in clusters II and III. These two MHC isoforms also had distinct relative patterns of regulation in the soleus and TA (Fig. 2C). Gene ontologies classifying genes in each of the other clusters are provided in Table S2. Cluster II has a common gene ontology in relation to muscle contraction, assembly and organization, with significant transcriptomic differences observed for soleus and TA in myosin binding protein C2 (Mybpc2) and myosin heavy chain 4 (Myh4) genes. In cluster III, Myl3, Myh7, Tnnc1 and Tpm3 describe the protein ensemble making up the muscle sarcomeres. The genes in the clusters I and IV do not share any common gene ontologies.

Because of the distinct regulatory differences of Myh4 and Myh7 genes in the soleus and TA, the expression levels of the genes coding for MHC isoform expression, specifically Myh1, Myh2, Myh4 and Myh7, were compared as a measure of fragments per kilobase of transcript per million mapped reads (FPKM) (Fig. 3A). Gene expression profiles of soleus muscle showed high levels of expression of the Myh7 gene and relatively lower levels of the faster myosin isoforms. High levels of Myh4 and Myh1 were identified in the TA, with low expression of Myh2 and Myh7 genes. Myh1 and Myh2 were expressed at lower levels in both muscle types compared with Myh4 and Myh7.

To specifically identify differences in MHC isoform distribution at the protein scale, SDS-PAGE was used. Our results revealed clear differences in myosin isoform distribution where the TA exhibited only fast isoforms while the soleus contained a mixture of fast and slow isoforms (Fig. 3C). Based on densitometric analysis of the gel bands and western blotting (not shown), the specific myosin isoforms were found to be statistically significant with average distributions of slow I (Myh7; 51%; P<0.05) and fast IIa (Myh2; 49%; P<0.05) in the soleus, and predominantly fast IIb (Myh4; 69%; P<0.05) and fast IIx (Myh1; 31%) in the TA (Fig. 3B). Despite low levels of Myh2 gene expression in the soleus, the protein expression of fast IIa myosin was much higher, indicating translation-level amplification of Myh2 protein production.

Mean shortening velocities for myosin isolated from soleus and TA, across increasing ATP concentrations, were measured using an in vitro motility assay (Fig. 4A). Consistent with myosin expression profiles (Figs 2 and 3), myosin from soleus had a significantly lower actin sliding velocity (V_actin) across all ATP concentrations when compared with myosin from TA, with a 3.3-fold decrease in velocity at a saturating ATP concentration of 1 mmol l⁻¹ (t=13.64; d.f.=3.645; P<0.005) (Fig. 4B). The Michaelis–Menten constant (K_m), which defines the ATP concentration required for half-maximal activation, was calculated to be 42±1.7 and 6.4±1.9 μmol l⁻¹ for TA and soleus myosin, respectively. This result is consistent with differences in the protein expression profile (Fig. 3), with TA consisting of only fast MHC isoforms whereas soleus consists of an almost even mix of fast and slow MHC isoforms.

Muscle ergometry studies on soleus and TA samples were used to account for tissue-scale factors that determine whole-muscle performance, and to delineate contributions of factors shaping active versus passive muscle functional attributes. Data from male and female samples were pooled for analyses (Figs 5 and 6) as GLM revealed no significant sex differences for all response variables. Significant differences between the soleus and TA were, however, observed in the active (Fig. 5A) and passive (Fig. 6A) FL and FVP (Fig. 7A) relationships.

Muscle mass-specific active tetanic force, measured at optimal length (L_o) was approximately 1.3-fold lower for soleus (18.0±4.4 N g⁻¹) than for TA (23.4±2.4 N g⁻¹) (Fig. 5A,D), a difference that was statistically significant (ANOVA; F_1,21=15.3; P<0.001). The FL plateau width, measured at 90% P_o on the ascending limb for the soleus (21.8±3.1% L_o) and TA (20±2.7% L_o), was not statistically significantly different (F_1,21=2.15; P=0.16) under tetanic activation (Fig. 5B,C,F). When compared with the twitch plateau width for the soleus (12.3±2.0% L_o) and TA (15.5±1.9% L_o), the tetanic plateau width was statistically significantly greater (soleus F_1,14=25.49; P<0.0001; TA F_1,12=9.14; P<0.01) and the left-shift of L_o with increased activation (twitch to tetanic) was 8.2±0.2% for both muscles and statistically significant (F_1,21=17.2; P<0.001).

Passive FL relationships are shown in Fig. 6A. We found passive force at 110% L_o to be statistically significantly different (F_1,19=61.86; P<0.0001), and lower for the soleus (0.012±0.005 N g⁻¹) than for the TA (0.091±0.035) (Fig. 6B). There was no statistically significant difference in the slack length (F_1,20=0.41; P<0.53) of soleus (−15±6.3% L_o) and TA (−16.7±5.6% L_o) (Fig. 6C). However, the passive FL curve slope was statistically significantly different (F_1,20=56.3; P<0.0001), suggesting that the passive soleus (β=0.05±0.029) is less stiff than the TA (β=0.8±0.36) across a physiologically relevant range of operating lengths (±10% L_o) (Fig. 6D). Because titin is the primary filament responsible for passive stiffness in the muscle sarcomere, we calculated the PSI for the soleus and TA titin genes using RNA-Seq analysis of data retrieved from the GEO database (GSE139204) to determine whether differences in passive muscle FL properties were related to titin expression differences. Alternative splicing-induced differences in the predicted size of expressed titin isoforms between the muscles revealed four regions of higher exon inclusion for the soleus compared with the TA (Fig. 6E), consistent with previous findings that demonstrate that slower muscles typically express longer titin isoforms (Prado et al., 2005).

The relationship between muscle mass-specific force and fiber length-specific shortening velocity provides key information about muscle performance. The FV relationship revealed a significantly (F_1,18=59.5; P<0.001) lower peak shortening velocity (V_max) for the soleus (4.6±0.5 L_o s⁻¹) compared with the TA (6.7±0.4 L_o s⁻¹) (Fig. 7A). The curvature of the FV relationship, as defined by the Hill constant (a/P_o), was statistically significantly (F_1,17=22.3; P<0.001) greater for the soleus (0.30±0.09) than for the TA (0.49±0.06) (Fig. 7A). Finally, peak muscle mass-specific power was statistically significantly (F_1,23=91.3; P<0.001) lower for the soleus (60.0±25.6 W kg⁻¹) than for the TA (202.2±48.3 W kg⁻¹), and optimal velocity (V_opt) and force (F_opt) – the velocity and force, respectively, where peak power is produced – were also statistically significantly (V_optF_1,23=39.2; P<0.001; F_optF_1,23=31.8; P<0.001) lower for soleus (1.1±0.3 L_o s⁻¹ and 5.6±1.6 N g⁻¹) than for TA (1.9±0.4 L_o s⁻¹ and 8.9±1.2 N g⁻¹), as shown by the left-shifted power–velocity peak for the soleus, when compared with the TA (Fig. 7B).

We took a multi-scale interdisciplinary approach to parse out the fundamental underpinnings of the protein ensemble involved in force and power generation across skeletal muscle lengths and speeds, specifically with respect to the thick filament motor protein myosin, which has a well-recognized and major role in determining muscle contractile performance. We studied mouse soleus and TA muscles, which are distinct in their fiber-type composition, architecture, energy metabolism, size and physiological function, to determine the differences between these muscles in MHC and titin isoform expression and distribution, in vitro translation of actin, and whole-muscle force–length–velocity–power properties. Our results highlight that myosin isoform-dependent differences in V_actin at the molecular scale does not completely account for differences in whole-muscle shortening velocities (V_max) between the TA and soleus. In addition, our multi-scale approach has provided a deeper understanding of the fine-tuning of transcriptional, translational, active and passive influences on whole-muscle force generation and power output between different muscles, culminating in specialized muscle function.

Our determination of the differences in gene expression profiles between the slow-twitch soleus and fast-twitch TA highlight the rich diversity that defines muscle fiber characteristics. PCA on soleus and TA transcriptomes revealed genes encoding myosin isoforms as the key determinant differentiating two muscle types (Fig. 2). Clusters II and III, identified through Ward.D2 cluster analysis, specifically showed opposite expression patterns between the soleus and TA, highlighting their significant transcriptome-based influence on each muscle type. On the one hand, over-expression of fast-type isoforms Mybpc2 and Myh4 in Cluster II highlights the influence of these genes in defining the fast-twitch phenotype of the TA, compared with the under-expression in the soleus. Significant variability in Cluster III genes includes opposite expression patterns in Myl3, Myh7, Tnnc1 and Tpm3 among the two muscle types (Fig. 2C). The genes Myl3 and Myh7 are typically associated with shortening velocity in slow cardiac tissue (Sitbon et al., 2021) and have high expression levels in the soleus, as compared with the TA. In addition, high levels of slow-twitch muscle-associated Tnnc1 and Tpm3 genes in Cluster III, and Tnni1 and Tnnt1 in Cluster IV for the soleus, also highlight the precise fine-tuning of isoform expression involved in altering muscle physiology. This is because Tpm3, Tnnc1, Tnni1 and Tnnt1 are critical in defining actin cytoskeletal assembly and Ca²⁺-based thin filament regulation (el-Saleh et al., 1986).

Quantification of the overall gene expression levels of the four MHC isoforms (Fig. 3A) showed that the TA is predominantly expressing Myh4, with lower levels of Myh1 and Myh2. However, the soleus predominantly expresses Myh7 along with lower levels of Myh2 and Myh1. At the protein scale, both the TA and soleus express varied myosin isoforms (Fig. 3C). There are only fast isoforms in the TA (∼2/3 fast MHC IIb and 1/3 fast MHC IIx), while the soleus has a close to 50:50 mix of fast MHC IIa and slow MHC I isoforms (Fig. 3B). Myosin protein levels measured here are consistent with previous studies showing that muscles with larger physiological cross-sectional area (PCSA) also tend to contain more fast fibers (Burkholder et al., 1994; Pellegrino et al., 2003). However, analysis of myosin isoform expression in mouse muscle during postnatal development has also demonstrated expression of low levels of fast MHC IIa (Agbulut et al., 2003). While we did not detect any MHC IIa via SDS-PAGE, we did detect low levels of fast MHC IIa via western blotting (not shown) and Tag-Seq (Fig. 3A). There is a large degree of variability regarding the correlation between expression levels of RNA and protein (Gry et al., 2009). Consistent with this previous observation, gene and protein expression levels were largely well aligned for the MHC isoforms expressed in the TA, despite no presence of fast MHC IIa at the protein scale. However, the gene expression levels for the soleus were ∼5-fold higher for Myh7 compared with Myh2, while the protein expression levels were close to equal. These comparative data highlight that transcription may not directly relate to protein translation levels in the muscle, necessitating measurement of both RNA and protein levels.

Expression of MHC isoforms has repeatedly been shown to influence muscle contractile properties (Bottinelli and Reggiani, 2000; Pellegrino et al., 2003; Reiser et al., 2009), and is also a major determinant of shortening velocity of single muscle fibers. Here, the ATP dependence of V_actin was determined using the in vitro motility assay (Fig. 4A). The 6.6-fold higher K_m observed for TA when compared with soleus represents a significantly lower affinity for ATP by the fast MHC isoforms in TA, along with a faster ADP release rate, thereby requiring a higher ATP concentration to reach half-maximal activation but with faster sliding velocities (Canepari et al., 2012). Consistent with the exclusive expression of fast myosin isoforms (Fig. 3), the V_actin at saturating ATP (1.0 mmol l⁻¹) for TA (3.91±0.04 µm s⁻¹) was in agreement with the V_actin measured for pure fast (MHC IIb) myosin (Pellegrino et al., 2003). In contrast, V_actin for soleus was 1.14±0.05 µm s⁻¹, close to that of pure slow I myosin (Pellegrino et al., 2003), despite this muscle having an almost 50:50 mixture of fast IIa and slow I MHC isoforms. These results indicate that differences in aggregate behavior of isoform distributions in specific muscles can be highly non-linear, with the slow isoforms having a dominating influence on shortening velocities compared with the fast isoforms, as shown previously (Egan et al., 2017; Harris et al., 1994).

The aggregate differences observed at the in vitro scale translate to the whole-muscle scale under the influence of active as well as passive muscle characteristics. At the actomyosin scale, slow myosin isoforms are capable of producing greater force than fast isoforms as a result of the prolonged attachment time of the slower isoforms (Cuda et al., 1997; Harris et al., 1994). This is consistent with the actin filament velocity we obtained for the mix of fast and slow isoforms found in the soleus, being closer to what would be seen for pure slow myosin (Pellegrino et al., 2003). Despite the apparent higher force production for the soleus myosin ensemble, the muscle mass-specific force measured for the soleus at the intact muscle scale was lower than for the TA. While our force measurements are comparable to those taken previously for these muscles (Close, 1965; Sacco and Jones, 1992), it cannot be ruled out that our different methods for testing the soleus (in vitro) and TA (in situ) affected our force measurements, despite the fact that we maintained the same temperature across preparations. The observed cross-scale discrepancy may also result from the more pennate structure of TA, which would result in more myofibers, and hence myosins, per cross-sectional area as compared with a less pennate muscle such as the soleus. These results highlight the need to consider muscle architecture when extrapolating from molecular- to tissue-scale experiments, as performed here. A related limitation of our study, therefore, is that we did not convert force measurements to stress (F/PCSA) because of the highly variable PCSA measurements for these small muscles in the literature (see Introduction).

Muscle force production is subject to change with respect to external factors. Most muscles rotate skeletal joints as they contract and therefore operate across varying muscle and sarcomere lengths. The link between joint range of motion and sarcomere length change may vary with activation level and can be relaxed by muscle fiber architecture and compliance. As muscle activation was reduced from tetanic to twitch, we observed a right-shift in optimal length by 8%, comparable to the 7% for mouse soleus and EDL (Askew and Marsh, 1997), and 5% (maximal to submaximal only) for rat MG (Tijs et al., 2021).

To gain a simplified understanding of and compare the range of operating lengths where the soleus and the more pennate TA produce near-maximal force, we measured the FL plateau width (Fig. 5B–D) (Winters et al., 2011). Despite fiber architectural differences, both muscles generated near-maximal force (≥90% P_o) across a broad range of tetanic operating lengths (±10% L_o). The plateau width at twitch activation was 8% narrower for both muscles, consistent with a previously reported plateau width reduction from maximal to submaximal activation in the rat MG (Tijs et al., 2021). Greater excursion abilities are commonly attributed to the in-series sarcomeres of relatively parallel-fibered muscles such as the soleus (Lieber and Ward, 2011; Winters et al., 2011) and are consistent with the relatively low stiffness of the soleus compared with the TA (Fig. 6D), so we predicted a greater plateau width for the soleus versus the TA. However, at the organ scale, the soleus attaches to the bones it moves with significantly less tendon than the TA, where added series elastic compliance may further decouple the length change of muscle fibers and sarcomeres from the excursion of skeletal joints (Konow et al., 2012). Because of muscle architectural differences, we acknowledge that lengthening both muscles by a percentage of whole-muscle length is unlikely to be representative of the changes occurring at the sarcomere level, and this difference may have caused some of the FL differences observed. Overall, however, our results underscore that structural differences between the soleus and TA transcend organizational scales to shape their respective functions in postural control (soleus) and limb-swing initiation (TA).

Understanding the differences in passive muscle properties between the soleus and TA also requires insights from multiple organizational scales. Passive muscle properties are driven by a combination of muscle architecture, titin isoform expression and extracellular matrix (ECM) organization and can be tuned by varying titin isoform length and/or the prominence of the ECM. Our results show that soleus stiffness is 16-fold lower (Fig. 6A,D) than that of the TA. This is consistent with the longer titin PEVK region that is expressed in soleus muscles (Fig. 6E). The PEVK region is the predominant region of alternative splicing in titin and is composed of two sequence motifs: the 28 amino acid PPAK motif and the variable length poly-E motif. The primary building block is the PPAK motif and the increased exons in the soleus muscle, as evidenced by the four regions of higher exon inclusion seen in Fig. 6E, result in a longer expressed titin isoform in the soleus compared with the TA. Therefore, consistent with observed correlations between PEVK length and passive tension in previous studies (Li et al., 2012; Prado et al., 2005), the increased length of titin in the soleus muscle likely underlies the muscle's low passive stiffness (Fig. 6).

Muscle architecture also contributes to this difference as muscles such as the TA with more pennate architecture typically have more prominent ECM, and therefore greater passive stiffness (Granzier and Irving, 1995; Purslow, 2020). Interestingly, our Tag-Seq data did not show significant differences in the transcript levels for the primary collagen genes between the soleus and TA (92±50 and 182±99, respectively, values expressed in FPKM). However, the increased ECM component in pennate muscles may be regulated at the translational or post-translational level.

Comparison of the force–velocity relationships (Fig. 7A) also highlights the influence of the myosin isoform and architecture differences on shortening velocities. For example, peak unloaded shortening velocity (V_max) was greater for the TA than for the soleus (6.7±0.4 and 4.6±0.5 L_o s⁻¹, respectively). Our V_max measurements for the soleus fell within the range measured previously (Asmussen and Maréchal, 1989; Holt and Askew, 2012) but we measured lower average V_max for the TA than a previous study (Burkholder and Lieber, 1998). However, consistent with our whole-muscle V_max measurements, myosin purified from TA muscle also exhibited a higher V_actin when compared with myosin purified from soleus muscle. Interestingly, the 3.4-fold difference between TA and soleus V_actin (V_actin TA/V_actin soleus) was dampened to a ∼1.6-fold difference in V_max at the whole-muscle scale (Fig. 7A), a depression that cannot directly be explained by measurement variability. Likely explanations for this difference include incomplete recruitment of all fibers during muscle contraction (Holt et al., 2014) and the presence of internal loads at the whole-muscle scale (Ross et al., 2020) that are not present with the protein isolates used for the in vitro motility assay. For instance, the TA is more pennate, and contains more ECM, than the soleus. Hence, the TA would likely have a larger internal (inertial) load, which may dampen V_max. It is also interesting to speculate that fiber pennation may be reducing the myosin-based velocity vector for TA (as a function of the cosine of the pennation angle), thus further reducing V_max, compared with the soleus.

The curvature of the Hill fit (Hill, 1938) applied to the muscle's force–velocity relationship (a/P_o) was greater for the soleus than for the TA (Fig. 7D), indicating that shortening velocity is more load sensitive for the soleus than TA. Previous studies have shown that load slows ADP release in striated muscle myosins, which would ultimately decrease cross-bridge detachment rate and thus fiber shortening velocity under load (Greenberg et al., 2014; Newhard et al., 2019; Weiss et al., 2001). Because slow myosins are more sensitive to load and mouse soleus has a 50:50 mix of slow MHC I and fast MHC IIb myosin, it is likely that the slow myosin isoform expression contributes to the increased force–velocity curvature.

The application of load to the ADP binding pocket of myosin is also influenced by differences in myosin light chain phosphorylation, which directly affects the load dependence of ADP release (Greenberg et al., 2009) as well as the loaded shortening velocities for MHC I (Greenberg et al., 2010; Karabina et al., 2015). Given the importance of the myosin light chains in supporting and stiffening the elongated alpha-helical neck region of myosin, it is interesting to speculate that the higher expression of Myl3 (Fig. 2C) may play a role in the increased sensitivity to load in soleus myosin.

By contrast, the TA showed a more consistent increase in velocity over all loads applied and therefore increased power compared with the soleus. The right-ward shift in peak power towards faster shortening velocities for the TA than for the soleus (Fig. 7B) illustrates the functional effects on increased fast myosin expression (Fig. 3C), consistent with the propensity for TA to prioritize power over force production (Hill, 1950). This propensity benefits TA during metabolically costly activities such as toe-off during running, compared with the lower energy consuming and slower contracting activities of the soleus, the predominant role of which is postural maintenance.

Skeletal muscle is often likened to a physiological ‘black box’ because of its multi-scale organization and function. Prevalent views of muscle as either diverse or constrained in its structure–function have prevailed in disciplines operating at different scales for decades. We propose that these divergent views are predominantly shaped by a lack of cross-scale research interaction. Twenty-first century convergence approaches can now let us crack the code of several biological black boxes, including muscle (e.g. Cavigelli et al., 2021). The utility of multi-scale approaches, as we have applied them here, likely extends to other areas in biology with complex interactions, ranging from microbial networks to ecosystem interactions. For muscle in particular, we anticipate leveraging the carefully characterized muscle types we have presented here and our multi-scale toolkit to determine, for example, how and why muscles with different fiber-type composition and architecture show different patterns of recovery from stretch injury and respond in different ways to changes in dietary nutrients.

Our integrative approach to studying muscles with disparate function identified the influence of myosin isoform distributions in determining muscle functional performance, along with additional, functionally relevant traits across scales of biological organization. We provide a comprehensive view of variations from gene expression to the protein scale, with respect to myosin isoform distributions, in addition to differences in gene expression signatures of associated proteins including titin, ECM and actin regulatory and cytoskeletal proteins. Apart from major differences in myosin isoform distribution at the protein scale, our analysis of the contributions of specific muscle architecture, as it relates to structural properties of protein–protein interactions, provide a deeper understanding of the cumulative effect of structural and functional traits in shaping overall muscle function. This understanding also highlights the need for interdisciplinary approaches. A deeper understanding of the multi-scale traits that shape muscle form and function will require an interdisciplinary toolkit but will also advance clinical and therapeutic efforts to harness the regenerative properties of skeletal muscle and extend its functionality beyond exercise injury and myopathies, and into older age of our longer-living population.

We thank Eric Hardiman, Krina Shah, Magdolin Shenouda and Devin Jenness for assistance with muscle experiments, and Jack Lepine for sequencing assistance. Thanks to students in the 2020 Biology of Muscle course at University of Massachusetts, Lowell for tracking actin filaments. Andrew Biewener and the Concord Field Station (Harvard University) were generous with facilities and resources, and Pedro Ramirez's animal care was invaluable. Two anonymous reviewers provided comments to improve our manuscript. This is publication #5 for the UMass Movement Center (UMOVE).