Gene Expression Profiling of Skeletal Muscles

Next-generation sequencing provides an opportunity for an in-depth biocomputational analysis to identify gene expression patterns between soleus and tibialis anterior, two well-characterized skeletal muscles, and analyze their gene expression profiling. RNA read counts were analyzed for differential gene expression using the R package edgeR. Differentially expressed genes were filtered using a false discovery rate of less than 0.05 c, a fold-change value of more than twenty, and an association with overrepresented pathways based on the Reactome pathway over-representation analysis tool. Most of the differentially expressed genes associated with soleus are coded for components of lipid metabolism and unique contractile elements. Differentially expressed genes associated with tibialis anterior encoded mostly for glucose and glycogen metabolic pathway regulatory enzymes and calcium-sensitive contractile components. These gene expression distinctions partly explain the genetic basis for skeletal muscle specialization, and they may help to explain skeletal muscle susceptibility to disease and drugs and further refine tissue engineering approaches.


Introduction
Skeletal muscles (SM) provide voluntary body movement and locomotion, posture and body position, energy production, fatty acid (FA) β-oxidation, carbohydrate metabolism, and soft tissue support. Skeletal muscles are composed of myofibers, which are formed by fused myoblasts. Myofibers have defining metabolic and contractile properties.
SM contribute to glucose, lipid, and protein metabolism. SM of healthy and well-fed mammals rely on carbohydrate or glucose metabolism as their primary sources of energy from aerobic respiration, anaerobic respiration, or both. SM can also use FAs as a fuel source via β-oxidation when mammals consume a high-fat diet or are in a starving state. As lipid resources become depleted, SM proteins can be broken down for energy use. Oxidative myofibers utilize mostly aerobic respiration. Glycolytic myofibers use primarily anaerobic respiration. Myofibers that utilize both are called oxidative-glycolytic.
Based on contractile properties, myofibers are defined by which myosin heavy chain is present. Myosin heavy chains are a significant component of the thick filament in a sarcomere. Myosin heavy chains (MyHCs) differentiate myofibers into types I (MYH7), II(A) (MYH2), II(B) (MYH4), and II(D or X) (MYH1) [1]. The four contractile properties relate to the three metabolic distinctions. Type I myofibers are oxidative [2]. Type II(A) myofibers are oxidative-glycolytic [2]. Types II(B) and II(D) myofibers are glycolytic [2]. Type II(D) myofibers combine with type II(A) or (B) to create hybrids called II(AD) and II(DB).
The hindlimb SM Soleus (So) and Tibialis Anterior (Ta) were used to determine the gene expression profile of oxidative and glycolytic myofibers by RNA-Seq in one-month-All animal experiments were performed following the Institutional and National Health and Medical Research Council guidelines. The experimental protocol was approved by the Institutional Animal Care and Use Committee at Oregon State University. Biopsies from Pitx2 FL/Z mice were utilized as previously described for this study [5].

RNA-Sequencing (RNA-Seq) Data Analysis
Five So and two Ta myofiber biopsies were collected at P30. RNA was extracted, sequenced, and analyzed as previously described [5]. The raw RNA-Seq data files for these samples are publicly available through the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) under the reference numbers SRP127367 (So) and SRP145066 (Ta). The reads were aligned using the mouse genome (mm10) and TopHat2. The genomic feature counts were then obtained using HTSeq [5], resulting in a total of 24,421 elements identified. The read counts were normalized to reduce transcript length bias [6]. Longer transcripts appear more prevalent and overshadow smaller transcripts, if not corrected before performing differential expression analysis. Quantitative changes in gene expression levels between the So and Ta groups were discovered using the quasilikelihood F test from the R package edgeR [7], which resulted in 6123 differentially expressed (DE) genes (p < 0.05 c, False Discovery Rate (FDR)-corrected) between So and Ta. The "c" in 0.05 c refers to the p-value being FDR-corrected in all instances. A smaller set of DE genes based on a log-fold change of ±1.0 or greater was used, resulting in 2481 DE genes that passed our fold change cutoff.

RNA-Seq Data Quality Check
So and Ta were processed and subject to RNA-Seq as described above. Principal component analysis of standardized samples showed a clear separation of So and Ta based on RNA expression (Figure 1a). We use red font to denote increased expression in So relative to Ta and blue for decreased expression in So relative to Ta. A volcano plot (negative log of the gene's FDR-corrected p-value vs. log 2 fold change) of all 24,421 annotated genes showed that the subset of 2481 genes included many highly significant DE genes with increased expression differences in So and Ta (Figure 1b). A heatmap of sample-standardized RNA counts of the final subset of 2481 genes shows clear clustering of myofibers from So and Ta based on gene expression (Figure 1c), as we expect from the biplot in Figure 1a.
Molecular pathways enriched for DE genes from the So and Ta samples were identified using the Reactome database [8]. The 2481 DE genes were matched to molecular pathways, and those pathways were statistically analyzed for over-or under-representation using the Reactome database and its corresponding ReactomePA R package [8,9]. The number of genes analyzed by the Reactome R package decreased to 2198 after removing genes with no Entrez ID. Hypergeometric testing, provided by the ReactomePA R package, identified 19 Reactome pathways as overrepresented with an adjusted p-value of less than 0.05 c. Three pathways had fewer than ten associated genes and were removed from further analysis. The DE genes related to the 16 Reactome pathways were ranked from highest to lowest log 2 (FC) values ( Figure S1). The ranking revealed that the 16 Reactome pathways could be consolidated into broader categories, such as contraction, ion and amino acid transport, and glucose and glycogen metabolism ( Figure S1). We identified 119 DE genes with FC ≥ 20, of which 14 were on the Reactome-curated list (Figure 1d). Examination of known gene functions of the remaining 105 high FC genes revealed that 36 could be manually categorized into lipid metabolism, glycogen metabolism, glucose metabolism, or contraction. standardized RNA counts of the final subset of 2481 genes shows clear clustering of myofibers from So and Ta based on gene expression (Figure 1c), as we expect from the biplot in Figure 1a. Molecular pathways enriched for DE genes from the So and Ta samples were identified using the Reactome database [8]. The 2481 DE genes were matched to molecular pathways, and those pathways were statistically analyzed for over-or under-representation using the Reactome database and its corresponding ReactomePA R package [8,9]. The number of genes analyzed by the Reactome R package decreased to 2198 after removing genes with no Entrez ID. Hypergeometric testing, provided by the ReactomePA R package, identified 19 Reactome pathways as overrepresented with an adjusted p-value of less than 0.05 c. Three pathways had fewer than ten associated genes and were removed from further analysis. The DE genes related to the 16 Reactome pathways were ranked from highest to lowest log2(FC) values ( Figure S1). The ranking revealed that the 16 Reactome pathways could be consolidated into broader categories, such as contraction, ion and amino acid transport, and glucose and glycogen metabolism ( Figure S1). We identified 119 DE genes with FC ≥ 20, of which 14 were on the Reactome-curated list (Figure 1d). Examination of known gene functions of the remaining 105 high FC genes revealed that 36 could be manually categorized into lipid metabolism, glycogen metabolism, glucose metabolism, or contraction. Genes with an FDR-corrected p-value (p_c) less than 0.05 c are gray dots, and those with a p-value higher than 0.05 c are black. A select set of highly significant DE genes with large fold changes for So (log FC > 3, −log(p_c) > 50) are shown in red and for Ta (log FC < −3, −log(p_c) > 75) in blue. (c) Heatmap of z-scored log-transformed read counts. Positive z-scores are in orange, and negative z-scores are in blue. Each sample is a column, and each row represents a statistically significant DE gene Genes with an FDR-corrected p-value (p_c) less than 0.05 c are gray dots, and those with a p-value higher than 0.05 c are black. A select set of highly significant DE genes with large fold changes for So (log FC > 3, −log(p_c) > 50) are shown in red and for Ta (log FC < −3, −log(p_c) > 75) in blue. (c) Heatmap of z-scored log-transformed read counts. Positive z-scores are in orange, and negative z-scores are in blue. Each sample is a column, and each row represents a statistically significant DE gene from the 2481 set. (d) Workflow used to reduce the number of genes for further analysis. A subset of the 2481 DE genes with extremely high Genes 2021, 12, 1718 4 of 19 differential expression was selected; 119 DE genes had a log-fold expression greater than twenty (purple). Separately, 97 of the 2481 genes mapped to sixteen statistically overrepresented molecular pathways using the Reactome database R package (p < 0.05 c) and had a log-fold expression change >3 (orange).

Lipid Metabolism
SM primarily metabolize glucose. Oxidative muscles can also use FAs along with glucose and glycogen as fuel sources. Free FA molecules are transported from the liver via cholesterol transport, which shuttles cholesterol, triglycerides, and protein-bound FAs. Three cholesterol-related genes, Pon1, Tspo2, and Apoa2, were expressed in So by 68, 51, and 27-fold, respectively ( Figure 2, Table 1). Pon1 protects against lipid oxidation for highdensity lipoprotein (HDL) cholesterol and low-density lipoprotein (LDL) cholesterol [10]. Tspo2 traffics free cholesterol in erythroid cells [11]. Apoa2 regulates steroid concentrations by modulating cholesterol transport [12], the precursor of steroid hormones.
Once free FA molecules pass through the plasma membrane, they enter the mitochondrial matrix to be broken down through β-oxidation. Specific protein channels handle different FA aliphatic tail lengths ranging from greater than 22 to fewer than 5 carbons. Acsm5 was expressed 30-fold higher in So ( Figure 2, Table 1). Acsm5 encodes for an enzyme that catalyzes the activation of FAs with aliphatic tails of 6 to 12 carbons by CoA to produce acyl-CoA, the first step in FA metabolism [13]. Dio1, expressed 69-fold higher in So ( Figure 2, Table 1), is involved in FA oxidation and oxidative phosphorylation uncoupling primarily in the liver, kidney, and thyroid [14], in addition to altering the thyroid hormone balance of triiodothyronine (T 3 ) and thyroxine (T 4 ) [15]. Thyroid hormone signaling regulates numerous genes involved in skeletal muscle homeostasis, function, and metabolism [16].
Several of the highest FC DE genes in So were involved in lipid, sphingolipid, and derived lipid biosynthesis. FA molecules can be elongated and modified into complex and derived lipid molecules. Tecrl, involved in FA elongation in polyunsaturated FA biosynthesis [17] and myoblast differentiation upon TNF activation [18], was expressed 55-fold higher in So ( Figure 2, Table 1). The increased expression of Tecrl could indicate that the So muscle is still developing in P30 mice.
One of the highest FC genes in this study was Cyp4f39 with a 269-fold higher expression in Ta, which encodes for a FA ω-hydroxylase involved in acyl-ceramide synthesis [19] ( Figure 2, Table 1). Cyp4f39 is involved in cell surface protection, cell recognition, signaling, membrane transportation, increased rigidity, hydrolyzation of very long-chain FAs, and steroid hormone synthesis [19]. Ceramides have been linked to insulin resistance in type II diabetes [20]. Positive or negative fold change represents increased gene expression in So or Ta, respectively.

Figure 2.
Lipid metabolism associated genes in So and Ta myofibers. Substrates and products of the enzymatic pathways for cholesterol transport, β-oxidation, FA elongation, and ceramide de-novo synthesis are shown based on their cellular location. The gene transcripts involved an enzymatic process and had increased expression in So or Ta, labeled in red or blue font, respectively. Font size correlates to relative FC.

Glycogen Metabolism
Glycogen metabolism combines two inverse processes: glycogenesis and glycogenolysis. During glycogenesis, glucose or glucose 6-phosphate molecules are converted into glycogen, using 1 adenosine triphosphate (ATP) and 1 ATP equivalent to add one glucose molecule to glycogen (Figure 3a). Only the glycolytic fibers store glycogen as a source of energy. Glycogen is tapped during times of nutritional insufficiency or rapid demand for Lipid metabolism associated genes in So and Ta myofibers. Substrates and products of the enzymatic pathways for cholesterol transport, β-oxidation, FA elongation, and ceramide de-novo synthesis are shown based on their cellular location. The gene transcripts involved an enzymatic process and had increased expression in So or Ta, labeled in red or blue font, respectively. Font size correlates to relative FC.

Glycogen Metabolism
Glycogen metabolism combines two inverse processes: glycogenesis and glycogenolysis. During glycogenesis, glucose or glucose 6-phosphate molecules are converted into glycogen, using 1 adenosine triphosphate (ATP) and 1 ATP equivalent to add one glucose molecule to glycogen (Figure 3a). Only the glycolytic fibers store glycogen as a source of energy. Glycogen is tapped during times of nutritional insufficiency or rapid demand for muscle activity. This study cannot determine the direction of glycogen metabolism alone, only inferring potentially emphasized steps in either So or Ta. The glycogenesis-involved enzymes phosphoglucomutase (Pgm2, Pgm2l1) and glycogen synthase (Gys2) had increased expression in Ta ( Figure 3a; Table 2). Phosphoglucomutase reversibly converts glucose 6-phosphate to glucose 1-phosphate. Increased expression of phosphoglucomutase could contribute to either glycogen metabolic processes.
Glycogen synthase is a critical enzyme in glycogenesis that converts glucose-1-phosphate into glycogen with the assistance of pyro-phosphorylase and a branching enzyme [25]. Increased expression of Gys2 in Ta may imply that glycogenesis is the priority for this muscle ( Figure 3a; Table 2). Gys2 activity is regulated by dephosphorylation and phosphorylation by protein phosphatase 1 (PP1) and cyclic-adenosine monophosphate (AMP)-dependent protein kinase. PP1 is an enzyme complex with regulatory and catalytic subunits, of which two regulatory subunits, Ppp1r1c and Ppp1r3g, were expressed higher in So ( Figure 3b; Table 2). Ppp1r3g regulates glycogen accumulation [26] and, along with the predicted inhibitory regulatory subunit Ppp1r1c, could be controlling glycogenesis in the So myofibers. The increased expression of the PP1 regulatory subunits in So could indicate a tighter control mechanism that prevents glycogen build-up. Epm2a, a factor that regulates PP1, was increased in Ta ( Figure 3b; Table 2). Epm2a encodes for laforin, which interacts with PP1 and regulates the dephosphorylation of glycogen to promote glycogen branching [27]. Laforin is also involved in preventing cytotoxicity by protein ubiquitination [28], and it may control glycogen synthesis by ubiquitinating PP1, which would prevent the activation of glycogen synthase.
Glycogenesis must be inhibited to trigger glycogenolysis or glycogen catabolism. Glycogenolysis neither requires nor produces energy when breaking glycogen into glucose molecules. The process is triggered by the accumulation of gluconeogenic precursors, such as lactate or alanine. Glycogen synthase is inactivated to inhibit glycogenesis. The inactivation triggers the activation of phosphorylase kinase, which phosphorylates glycogen phosphorylase. Phosphorylase kinase is an enzyme complex of four subunits [29], three of which (Phka1, Phkg1, Phkb) also had increased expression in Ta ( Figure 3c; Table 2). Phosphorylase kinase could be regulating the glycogenolysis in Ta without inhibiting glycogen synthase by ubiquitination. Glycogen phosphorylase pairs with a debranching enzyme to break off glucose 1-phosphate molecules from glycogen; a gene encoding for each (Pygm and Agl, respectively) was expressed higher in Ta (Figure 3a; Table 2). Phka1 phosphorylase kinase α 1 −8 activates glycogen phosphorylase [29] Phkb phosphorylase kinase β −7.5 activates glycogen phosphorylase [29]  Positive or negative FC represents increased gene expression in So or Ta, respectively.

Glucose Metabolism
Glucose metabolism consolidates glycolysis and gluconeogenesis enzymatic processes. The former process breaks glucose into an intermediate pyruvate, and the latter is the reverse process. Glycolysis converts glucose into pyruvate in 10 steps, consuming one glucose, two NAD+, two Pi, four ADP, and two ATP molecules to generate two pyruvate, two NADH, two H+, four ADP, and four ATP molecules. The net products are two pyruvate, two ATP, two NADH, and two H+ molecules. Gluconeogenesis converts two pyruvate, four ATP, two GTP, two NADH, two H+, and two H2O molecules into fructose 6-phosphate or glucose 6-phosphate. Overall, gluconeogenesis consumes six ATP, and glycolysis generates a net two ATP. Therefore, the whole glycolysis/gluconeogenesis cycle costs four ATP.
Three critical enzymes regulate the complete glycolysis process, and they are hexokinase, phosphofructokinase-1 (PFK1), and pyruvate kinase. Genes encoding for PFK1 (Pfkm) and pyruvate kinase (Pkm) had increased expression in Ta (Figure 4a; blue font, Table 2). PFK1 converts fructose 6-phosphate to F1,6-BP and is activated by high AMP concentrations and fructose 2,6-bisphosphate (F2,6-BP) [31]. Pyruvate kinase converts phosphoenolpyruvate to pyruvate when activated by high F1,6-BP concentration and is deactivated by an elevated ATP concentration. Four other isoenzymes involved in glycolysis and gluconeogenesis, aldolase (Aldoa, Aldoc), triosephosphate isomerase (Tpi1), phosphoglycerate kinase (Pgk1), and enolase (Eno3 in Ta and Eno2 in So), were increased in Ta (Figure 4a, Table 2). Represented are the addition (+Ca 2+ ) and the removal of calcium ions (−Ca 2+ ). Substrates and products of the enzymatic steps in glycogen metabolism include the interconversion of inactive to active states of regulatory enzymes. Enzymes involved at specific stages are encircled and adjacent to the corresponding reaction arrow. Gene symbols located underneath enzyme names in red or blue font were increased in So or Ta, respectively. Font size correlates to relative FC.

Glucose Metabolism
Glucose metabolism consolidates glycolysis and gluconeogenesis enzymatic processes. The former process breaks glucose into an intermediate pyruvate, and the latter is the reverse process. Glycolysis converts glucose into pyruvate in 10 steps, consuming one glucose, two NAD+, two Pi, four ADP, and two ATP molecules to generate two pyruvate, two NADH, two H+, four ADP, and four ATP molecules. The net products are two pyruvate, two ATP, two NADH, and two H+ molecules. Gluconeogenesis converts two pyruvate, four ATP, two GTP, two NADH, two H+, and two H2O molecules into fructose 6phosphate or glucose 6-phosphate. Overall, gluconeogenesis consumes six ATP, and glycolysis generates a net two ATP. Therefore, the whole glycolysis/gluconeogenesis cycle costs four ATP.
Three critical enzymes regulate the complete glycolysis process, and they are hexokinase, phosphofructokinase-1 (PFK1), and pyruvate kinase. Genes encoding for PFK1 (Pfkm) and pyruvate kinase (Pkm) had increased expression in Ta (Figure 4a; blue font, Table 2). PFK1 converts fructose 6-phosphate to F1,6-BP and is activated by high AMP concentrations and fructose 2,6-bisphosphate (F2,6-BP) [31]. Pyruvate kinase converts phosphoenolpyruvate to pyruvate when activated by high F1,6-BP concentration and is deactivated by an elevated ATP concentration. Four other isoenzymes involved in glycolysis and gluconeogenesis, aldolase (Aldoa, Aldoc), triosephosphate isomerase (Tpi1), phosphoglycerate kinase (Pgk1), and enolase (Eno3 in Ta and Eno2 in So), were increased in Ta (Figure 4a, Table 2). Represented are the addition (+Ca 2+ ) and the removal of calcium ions (−Ca 2+ ). Substrates and products of the enzymatic steps in glycogen metabolism include the interconversion of inactive to active states of regulatory enzymes. Enzymes involved at specific stages are encircled and adjacent to the corresponding reaction arrow. Gene symbols located underneath enzyme names in red or blue font were increased in So or Ta, respectively. Font size correlates to relative FC.
When blood glucose concentrations are low, glucose can be generated through gluconeogenesis, an 11-step process for converting pyruvate to glucose. The enzymes pyruvate carboxylase, PEP carboxykinase, and fructose 1,6-bisphosphatase (FBPase1) are critical regulatory points in gluconeogenesis. FBPase1 (Fbp2) and PEP carboxykinase (Pck1) were increased in Ta and So, respectively (Figure 4a, Table 2). FBPase1 converts F1,6-BP to fructose 6-phosphate in the presence of a high concentration of AMP and F2,6-BP, while increased citrate levels deactivate it. PEP carboxykinase turns oxaloacetate into phosphoenolpyruvate and is enabled by a high ADP concentration.
All myofibers utilize glycolysis, whether or not they use aerobic and/or anaerobic respiration. Glycolysis occurs under both aerobic and anaerobic conditions because it does not require oxygen. Most of the DE genes associated with glucose metabolism had prominent expression in Ta compared to So. As a muscle that uses anaerobic metabolism to gain fuel, Ta may better utilize the anaerobic pathway due to the increased expression levels of glycolytic and regulatory enzymes.
(Pck1) were increased in Ta and So, respectively (Figure 4a, Table 2). FBPase1 converts F1,6-BP to fructose 6-phosphate in the presence of a high concentration of AMP and F2,6-BP, while increased citrate levels deactivate it. PEP carboxykinase turns oxaloacetate into phosphoenolpyruvate and is enabled by a high ADP concentration.
All myofibers utilize glycolysis, whether or not they use aerobic and/or anaerobic respiration. Glycolysis occurs under both aerobic and anaerobic conditions because it does not require oxygen. Most of the DE genes associated with glucose metabolism had prominent expression in Ta compared to So. As a muscle that uses anaerobic metabolism to gain fuel, Ta may better utilize the anaerobic pathway due to the increased expression levels of glycolytic and regulatory enzymes.

Contraction
The Reactome pathway overrepresentation analysis identified contraction as another significant difference between So and Ta. The sliding filament theory explains the molecular basis of muscle contraction, where myosin and actin filaments interact in a particular manner to produce contractile force. Most contraction-related DE genes exposed were associated with the thin filament, thick filament, and Z-line of the sarcomere (Figure 5a, Table 3).
The thin filament is composed of actin subunits. Actc1 had increased expression in Ta (Figure 5b, Table 3). Tropomyosin wraps outside the actin helix to stabilize the strand and rotate the thin filament to expose the binding sites. Two isoforms of tropomyosin were increased in So (Tpm3) [32] or Ta (Tpm1) [33] (Figure 5b, Table 3). The troponin complex is a regulatory element of the actin filament and is associated with calcium binding, inhibitory regulation, and tropomyosin binding. Tropomyosin-binding element encoding genes Tnnt1 [34] and Tnnt2 [35] were increased in So, while Ta had increased expression of Tnnt3 [36] (Figure 5b, Table 3). The calcium and inhibitory troponin units, Tnnc1, Tnni1, and Tnni3, had an increased expression in So (Figure 5b, Table 3). The actin filaments are capped by tropomodulin Tmod1, which was expressed higher in Ta (Figure 5b, Table 3).
The thick filament consists of several strands of heavy myosin chains, which are two myosin heavy chains wrapped around one another with myosin heads at the ends that attach and flex to move along the thin filament. Myosin heavy chain isoforms unique to type I (Myh7) [1] and type II(B) (Myh4) [37] myofibers were increased in So and Ta, respectively (Figure 5c, Table 3).
Myosin light chains and myosin-binding proteins regulate myosin head movement. Myosin light chains are classified into essential (ELCs) and regulatory (RLCs). ELCs stabilize the myosin head, while RLCs stiffen the myosin neck domain. ELCs encoded by Myl6b, Myl3, and Myl4 [2] and the RLCs encoded by Myl2, Myl10 [38] were increased in So (Figure 5c, Table 3) while only Mylpf [39] was increased in Ta. Mybpc2 [40] and Mybph that encode for myosin-binding proteins that assist in movement along the thin filament were increased in Ta (Figure 5c, Table 3).
The thin and thick filaments are held in an antiparallel configuration by a structure called the Z-line at the sarcomere ends (Figure 5d). The titin cap (Tcap) [41] is involved in mechano-electrical links between Z-lines and T-tubules and was expressed higher in So (Figure 5d, Table 3). Actinin (Actn3) [42] is involved in stabilizing the antiparallel formation and was increased in Ta (Figure 5d, Table 3). These structural differences could be altering the electrical activity and, in turn, the contractile properties of the different myofibers in the So and Ta muscles.

Discussion
Two broad-scale patterns emerged along with the three categories of overrepresented molecular pathways when examining the gene expression profiles of So and Ta muscles ( Figure S1). First, DE genes involved in lipid, glucose, and glycogen metabolism pathways associated more with one muscle than the other (Tables 1 and 2). Second, each muscle expressed unique genes encoding for similarly functioning isoenzymes, particularly for genes involved in the contraction and ion transport categories (Tables 3 and S1). We put these differences into context below, starting with the pathways that differentiate the muscle types. DE genes with the highest FC difference were not prevalent among genes associated with the overrepresented pathways. Therefore, DE genes with a greater than 20-fold expression difference were included for further investigation. Similarly seen in a previous study [54], we found that approximately 10% of DE transcripts were mainly related to fatty acid metabolism, structural components, and neuromuscular junction assembly.

Lipid Metabolism
FA lipids are an alternative energy source during fasting, starvation, and endurance exercise in oxidative muscles. Both So and Ta muscles contain oxidative myofibers, and both should express genes associated with FA metabolism (Figure 2). However, DE genes explicitly associated with FA transport and catabolism were increased in expression in So, which contains a higher percentage of oxidative fibers in its overall myofiber composition. In contrast, the 269-fold increased expression of Cyp4f39 in Ta suggests that lipid catabolism to generate ceramides is more active in Ta. Ceramide accumulation has been linked to insulin resistance in type II diabetes [55]. The increased Cyp4f39 expression may allow Ta to utilize ceramides to shift from being insulin-sensitive and from facilitating efficient glucose uptake [20] to becoming insulin-resistant. Alternatively, increased expression of Cyp4f39 and presumed increased ceramide levels in Ta might be associated with differential apoptosis, cell cycling, or autophagy.

Glycogen and Glycosaminoglycan Metabolism
Glycogen acts as an energy reserve for myofibers by storing excess glucose molecules that are then utilized when blood glucose levels are low. Glycolytic or oxidative-glycolytic myofibers retain glycogen because they rely heavily on glucose as energy, and Ta has a high percentage of those myofibers; therefore, the increased expressions of the glycogen metabolic enzyme genes Gys2, Pygm, Phka1, Phkg1, and Phkb in Ta support the known metabolic composition of the muscle relative to So (Figure 3a; Table 2).
Hydrolysis of muscle glycogen to glucose occurs in lysosomes that engulf glycogen granules. The lysosome-associated DE genes, Slc35d3, Atp6v0d2, and Galns, exhibited increased expression in So, suggesting that lysosomal-related organelles in So and Ta myofibers perform slightly different functions (Table S1). Slc35d3 is associated with the biosynthesis of platelet-dense granules [56]. Platelet-dense granules contain high concentrations of calcium, adenine nucleotides, pyrophosphate, and polyphosphate molecules that enhance autophagy in lysosome-related organelles. Atp6v0d2 is a part of vacuolar ATPases involved in proton translocation into vacuoles, lysosomes, or the Golgi apparatus to lower the pH [57]. The presence of Atp6v0d2 suggested that lysosome-related organelles in So need to reduce their pH levels routinely. If lysosome-related organelles in So myofibers have an enhanced autophagy, the pH levels within those organelles would fluctuate as the cell engulfs and metabolizes more molecules from the extracellular space and may require proton pumps.
Galns is involved in glycosaminoglycan biosynthesis. Glycosaminoglycans are cell surface proteins with branches of sugars projecting into the extracellular matrix to support cell identity, adhesion, and growth. Galns encodes for a protein that breaks keratan sulfate off glycosaminoglycans [58] in lysosomes. Stab2, B3gat1, and Cspg5 are also involved in glycosaminoglycan synthesis, indicating that So and Ta synthesize distinct glycosaminoglycans associated with cell adhesion and proliferation (Tables S1 and S2). Stab2 is involved in the endocytosis of metabolic waste products, including circulating hyaluronic acid (HA) that promotes cell proliferation once the plasma concentration decreases. B3gat1 is involved in generating HNK-1 carbohydrate (CD57) cell surface epitope associated with cell adhesion [59]. Cspg5 is involved in chondroitin sulfate synthesis, another molecule added to proteoglycans and connected to cell adhesion, growth, migration, and receptor binding in the central nervous system (Table S1). The high Cspg5 expression in So suggests that neurons associated with the So muscle are marked differently.

Glucose Metabolism
Genes associated with glucose metabolism were increased in Ta more than So despite both muscles utilizing glycolysis. Yet, Eno2 and Pck1 were increased in expression in So compared to the 11 isoenzymes in Ta, perhaps due to the higher percentage of type II(B) glycolytic myofibers in Ta ( Figure 4, Table 2). Glycolysis-related enzymes encoded by Pfkm, Pkm Pfkfb3, Pfkfb1, and Pfkfb4 were increased in Ta (Figure 4b, Table 2). In contrast, two gluconeogenesis regulatory enzymes that regulate the reversion of pyruvate to glucose in skeletal muscle encoded by Fbp2 and Pck1 were expressed higher in Ta and So, respectively. Gluconeogenesis in Ta might not be as tightly controlled across several steps as glycolysis seems to be based on the 3-to-1 regulatory enzyme ratio. Gluconeogenesis may be controlled by the PFK2/FBPase2 complex instead of relying on the energy-costing enzyme PEP carboxylase seemingly utilized by So. Pfkfb4 had a higher DE than Pfkfb3. Pfkfb4 opposes Pfkfb3 as it redirects glucose to the pentose phosphate pathway to promote the detoxification of reactive oxygen species and lipid and nucleotide biosynthesis. Fbp2, which was also increased in expression in Ta, converts F1,6-BP back into fructose 6-phosphate. These combined findings may indicate that Ta initiates gluconeogenesis at the PFK2/FBPase2 step of the glucose metabolism process to reduce energy loss.

Contraction
Several components that contribute to sarcomere formation were identified as DE between So and Ta, including MyHCs, myosin light chains, troponin subunits, and NMJassociated proteins (Figure 5a-e). The presence of type I (Myh7) and type II(B) (Myh4)associated MyHC isoforms among DE genes provided proof of the concept for this analysis. The fold change in differential expression of Myh7 and Myh4, 34.7-and 57.7-fold, respectively (Table 3), reflected the appropriate percentage of types I and II(B) myofibers in each muscle. So and Ta both have type II(A) and II(D) myofibers, so differential expression of their specific MyHCs was predictably absent from our analysis. Our data also support a previous analysis correlating Myh4 and Myh7 with other genes associated with myofiber types I and II(B) [54], respectively (Tables 3, S1 and S2). The myosin ELCs and RLCs mainly correlated to the appropriate muscle, except for Myl4 and Myl6 (Figure 5c, Table 3). The troponin complex subunits encoded by Tnnt1, Tnni1, and Tnnc1 were unique and highly expressed in So (Figure 5b, Table 3), suggesting that the troponin complex in So might be modified relative to Ta to keep the actin filament in an open position longer, facilitating longer durations of contraction.
Factors affecting NMJ specification and calcium handling have been theorized to be among the non-myofibrillar, muscle-specific systems that allow for increased maximal isometric stress after P28 [60]. Synaptic modeling affects the motility or twitching of the skeletal muscle. Muscle identity is determined before innervation, yet innervation maintains muscle identity. During regeneration, only innervated So muscles can upregulate slow isoform mRNAs [61]. Yet, muscles do not completely change into another phenotype when neural cues are altered [62]. Intrinsic signaling from nerves and other sources, and environmental cues, dictate adult muscle phenotype [61]. Genes are associated with neurotransmission and synaptic vesicles [43] at the presynaptic side of NMJs, in addition to clathrin-mediated endocytosis [47,49], and cell adhesion [48] of the postsynaptic side was highly expressed in So. The release and retrieval of the synaptic vesicles and their contents may play an important role in the slow-twitch mechanism. On the contrary, Ta had increased expression of genes associated with the presynaptic side of NMJs only, which modulate neuronal excitability [44], broaden the action potential [45], and slow down the degeneration of NMJs [46]. These genes may explain how Ta and its associated neurons handle the stimulation necessary for rapid contraction. The DE genes related to the ion transport pathways may be associated with the postsynaptic side of NMJs in Ta myofibers, such as the voltage-dependent calcium channel subunits (Table S1). However, all of this information and conjuncture is based on a small number of differentially expressed genes.

Conclusions
Two hindlimb muscles, So and Ta, which perform distinct locomotion functions, are characterized by unique metabolic and contractile properties. This focused investigation identified several genes uniquely tied to either So or Ta when these two muscle groups are compared. Approximately 10% of the mouse transcriptome was differentially expressed, including genes involved in lipid, glucose, and glycogen metabolism. Very highly differentially expressed genes such as Cyp4f39 highlight previously undescribed potential differences in fatty acids between the SM that may underlie susceptibility to pathologies. The high DE of Pfkm and Gys2 point to known differences in glucose and glycogen metabolism between the SM but add molecular details to the picture. Among contractionrelated genes, RLCs, ELCs, and troponin complex subunits were higher expressed in So. DE of genes involved in NMJs suggests that electrochemical signaling is different between the two muscle groups, supporting the contractile differences.
Overall, this detailed study into the gene expression differences between these two muscle groups adds to the abundance of data on molecular differences among SM. Some similar results were found in single nuclei studies, which also identified unique expressions among individual skeletal muscle groups [54,63]. This study also highlights the continuing need to maintain and update gene and pathway ontology databases, especially adding more information on gene expression in different cell types.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/genes12111718/s1, Figure S1: Evaluation of Reactome overrepresented pathway analysis, Table S1. Differentially Expressed Genes in So and Ta involved in Ion Transport and Glycosaminoglycan biosynthesis pathways, Table S2. Differentially Expressed Genes in So and Ta organized by Immune Response, Signaling, and Cellular function pathways.