Bioenergetic defects in muscle fibers of RYR1 mutant knock-in mice associated with malignant hyperthermia

Mutations in the skeletal muscle ryanodine receptor gene ( RYR1 ) can cause susceptibility to malignant hyperthermia (MH), a potentially lethal genetic condition triggered by volatile anesthetics. MH is associated with hypermetabolism, which has directed research interest into oxidative phosphorylation and muscle bioenergetics. The most common cause of MH in the United Kingdom is the c.7300G > A RYR1 variant, which is present in ~ 16% of MH families. Our study focuses on the MH susceptible G2435R-RYR1 knock-in mouse model, which is the murine equivalent of the human c.7300G > A genotype. Using a combination of transcriptomics, protein expression, and functional analysis, we investigated adult muscle fiber bioenergetics in this mouse model. RNA-Seq data showed reduced expression of genes associated with mitochondria and fatty acid oxidation in RYR1 mutants when compared with WT controls. Mitochondrial function was assessed by measuring oxygen consumption rates in permeabilized muscle fibers. Comparisons between WT and homozygous G2435R-RYR1 mitochondria showed a significant increase in complex I – facilitated oxidative phosphorylation in mutant muscle.

Mutations in the skeletal muscle ryanodine receptor gene (RYR1) can cause susceptibility to malignant hyperthermia (MH), a potentially lethal genetic condition triggered by volatile anesthetics. MH is associated with hypermetabolism, which has directed research interest into oxidative phosphorylation and muscle bioenergetics. The most common cause of MH in the United Kingdom is the c.7300G>A RYR1 variant, which is present in~16% of MH families. Our study focuses on the MH susceptible G2435R-RYR1 knock-in mouse model, which is the murine equivalent of the human c.7300G>A genotype. Using a combination of transcriptomics, protein expression, and functional analysis, we investigated adult muscle fiber bioenergetics in this mouse model. RNA-Seq data showed reduced expression of genes associated with mitochondria and fatty acid oxidation in RYR1 mutants when compared with WT controls. Mitochondrial function was assessed by measuring oxygen consumption rates in permeabilized muscle fibers. Comparisons between WT and homozygous G2435R-RYR1 mitochondria showed a significant increase in complex I-facilitated oxidative phosphorylation in mutant muscle. Furthermore, we observed a genedose-specific increase in reactive oxygen species production in G2435R-RYR1 muscle fibers. Collectively, these findings provide evidence of metabolic defects in G2435R-RYR1 knock-in mouse muscle under basal conditions. Differences in metabolic profile could be the result of differential gene expression in metabolic pathways, in conjunction with mitochondrial damage accumulated from chronic exposure to increased oxidative stress.
The ryanodine receptor isoform 1 (RyR1) is a calcium release channel expressed predominantly in the sarcoplasmic reticulum (SR) of skeletal muscle. Mutations in the RYR1 gene result in dysregulation of intracellular Ca 21 homeostasis and are the primary genetic cause for malignant hyperthermia (MH) susceptibility, a dominantly inherited disorder that presents with potentially fatal hypermetabolic reactions after exposure to volatile anesthetic agents or succinylcholine (1). Clinical diagnosis of MH susceptibility uses a combination of genetic screening and pharmacological challenge testing (2). At present, the European Malignant Hyperthermia Group recognizes 48 functionally characterized RYR1 variants for use in prospective DNA diagnosis of MH (www.emhg.org). Within this list, the RYR1 c.7300G.A variant is the most prevalent within the United Kingdom, found in almost 16% of MH families (3). This missense variant results in a p.G2434R amino acid change which, despite being the most common variant, generally presents with a weaker phenotype compared with other variants (4). Given its prevalence, a new mouse RyR1 knock-in model (G2435R-RYR1), equivalent to the human c.7300G.A genotype, was recently developed to further our understanding of MH (5).
Individuals who are susceptible to MH may have no phenotype without exposure to anesthesia or can present with a range of muscle-related features such as myopathy, exertional heat illness, rhabdomyolysis, exercise intolerance, weakness, and cramping (6,7). Whereas the primary feature of the hypermetabolic reaction during anesthesia is a massive increase in [Ca 21 ] i , mitochondrial dysfunction, which has been found in other MH mouse models (Y522S-RYR1, R163C-RYR1, and T4826I-RYR1) (8)(9)(10), may contribute to the nonanesthesia phenotypes.
Structural anomalies such as mitochondrial swelling and deformation, which appear to increase in severity with age, have been observed in homozygous (HOM) T4826I-RYR1 and heterozygous (HET) Y522S-RYR1 mice (8,11,12). Evidence of reduced mitochondrial content has been seen in all HET RYR1 knock-in MH mouse models investigated (9)(10)(11). The mitochondrial injury observed in HET Y522S-RYR1 and HET R163C-RYR1 muscle has been attributed to increased oxidative stress because of increased production of reactive oxygen species (ROS) (8,9,13). All of these murine models have also been shown to have hypermetabolism after exposure to an anesthetic challenge and elevated resting calcium ([Ca 21 ] rest ) in their skeletal muscle fibers (12,14,15).
However, there are limitations in the methods previously used to study mitochondrial function. Studies using isolated mitochondria are limited by the absence of their normal cellular environment and lack of interaction with other organelles. This can be detrimental in the study of conditions such as MH because it is intimately dependent on SR function and calcium regulation. Likewise, mitochondrial studies in cultured myotubes, which are immature myofiber precursors, have limitations because of differences in their calcium signaling and organelle structure compared with adult muscle fibers (16,17 the G2435R-RYR1 MH-susceptible mouse model using a combination of transcriptomics, protein expression, and muscle fiber functional analysis. Similar to previously characterized RYR1 knock-in mouse models, HET and HOM G2435R-RYR1 mice have elevated [Ca 21 ] rest and fulminant MH reactions in response to anesthetic exposure (5). We performed respirometry and real-time oxidative stress analysis on adult muscle fibers, allowing measurements to be made within the context of native muscle structure and signaling properties. We hypothesized that G2435R-RYR1 mouse muscles will exhibit differential gene expression governing metabolic pathways that contribute toward mitochondrial dysfunction. We used HET T4826I-RYR1 knock-in mice for comparison because they are a well-documented model with previously demonstrated alterations in their metabolic profile (10,12).

Differential gene expression
RNAseq was performed on poly(A)-selected mRNA samples extracted from WT, G2435R-RYR1 HET, G2435R-RYR1 HOM, and T4826I-RYR1 HET mice. A sample-sample distance heatmap (Fig. 1A) was used to visualize overall gene expression between samples, and through the dendrogram we can see clustering of WT and G2435R-RYR1 HET samples showing a close relationship between these two genotypes. In contrast, the T4826I-RYR1 HET samples formed their own distinct cluster alongside a pair of G2435R-RYR1 HOM samples, which suggests that these two genotypes have more divergent transcriptome profiles than G2435R-RYR1 HETs compared with WT. Principal component analysis (Fig. 1B) displays highest degrees of variability between the T4826I-RYR1 HET and G2435R-RYR1 HOM samples as shown by the distant clustering on the PC1 axis (31% variance). Variability on the PC2 axis is greatest within the G2435R-RYR1 HOM genotype because one sample does not seem to cluster (26% variance).
Transcriptomes from each RYR1 knock-in genotype were compared with WT to identify differentially expressed genes. Differences with adjusted p values , 0.05 were considered significant following the protocol outlined in the DESeq2 R package. In the WT versus G2435R-RYR1 HOM and WT versus T4826I-RYR1 HET comparisons, a cohort of 393 (98 up-regulated and 295 down-regulated) and 206 (82 up-regulated and 124 down-regulated) differentially expressed genes were found, respectively. In contrast, the WT versus G2435R-RYR1 HET knock-in comparison showed only one gene differentially expressed, Blvrb, which was found to be up-regulated in all RYR1 knock-in genotypes.

Pathway and gene ontology analysis
Pathway and gene ontology analyses were conducted on WT versus G2435R-RYR1 HOM and WT versus T4826I-RYR1 HET gene lists. Pathway analysis of down-regulated genes in the WT versus G2435R-RYR1 HOM comparison showed a set of seven significantly enriched pathway terms, many of which are involved in cellular energy production ( Table 1). Within this set, five pathway terms, "fatty acid b-oxidation WikiPathways1269," "mitochondrial long-chain (LC) fatty acid b-oxidation WP401," "fatty acid biosynthesis WP336," "tricarboxylic acid (TCA) cycle WP434," and "fatty acid oxidation WP2318," were also found down-regulated in the WT versus T4826I-RYR1 HET comparison, but only the last term was statistically significant. The gene ontology analysis of down-regulated genes showed that the top 10 enriched biological processes from the WT versus G2435R-RYR1 HOM comparison were all involved in either fatty acid oxidation or fatty acid transport into mitochondria ( Table 2). The gene ontology analyses of down-regulated genes in the WT versus T4826I-RYR1 HET gene list shared several of the same terms surrounding fatty acid metabolism but were not statistically significant (adjusted p values . 0.05) ( Table 2). The only shared observation between the WT versus G2435R-RYR1 HOM and the WT versus T4826I-RYR1 HET comparisons was an enrichment of genes in the cellular component gene ontology terms "mitochondrion" and "mitochondrial matrix" (Table S1 and Table S2). The WT versus G2435R-RYR1 HOM comparison featured the highest number of differentially expressed genes in the mitochondrion term. Among these, genes of interest include Ndufs1, Ndufs2, Sdha, and Uqcrc2, which encode for mitochondrial complexes I, II, and III of the electron transport chain (Table S1). In contrast, analysis of up-regulated genes from each comparison showed Table 1 Pathway analysis using Enrichr. The top 10 downregulated pathways are shown for the WT versus G2435R-RYR1 HOM and WT versus T4826I-RYR1 HET comparisons, in combined score order (log of p-value from Fisher's exact test multiplied by z-score of deviation from expected rank). This pathway analysis was produced using significantly downregulated genes (adjusted p-values < 0.05) generated from DESeq2, as input for Enrichr. Statistically significant pathway terms are annotated with an asterisk (adjusted p-values < 0.05)  Table 2 Gene ontology enrichment analyses using Enrichr. The top 10 most enriched biological processes (GO_BP) are shown, generated using downregulated gene lists from both the WT versus G2435R-RYR1 HOM and WT versus T4826I-RYR1 HET comparisons. Statistically significant ontology terms are annotated with an asterisk (adjusted p-value < 0.05) Additional results from these analyses can be found in the supporting information (Tables S1-S4) Regulation of acyl-CoA biosynthetic process (GO:0050812) Cardiolipin metabolic process (GO:0032048) 1.00E 1 00 Lclat1, Taz Acetyl-CoA metabolic process (GO:0006084) 1.00E 1 00 Acss1, Nudt7 Regulation of acetyl-CoA biosynthetic process from pyruvate (GO:0010510) Fatty acid b-oxidation using acyl-CoA dehydrogenase (GO:0033539) 1.00E 1 00 Acadm, Etfa no significant enrichment of pathways and gene ontology terms, aside from a few cellular component gene ontology terms found in the WT versus G2435R-RYR1 HOM comparison (Table S3 and Table S4).
Oxygen flux was then normalized to a common reference state, CI1CII (ETS) , to obtain flux control ratios (FCR) (Fig. 2B). These ratios provide a measure of function as a proportion of ETS capacity, which is independent of mitochondrial content. In contrast to raw data, a comparison of normalized data showed that FCR differed significantly among RYR1 genotypes in the LEAK, CI (OXPHOS), and CI1CII (OXPHOS) states (all p , 0.001). Pairwise FCR comparisons within the LEAK state did not differ from raw data after normalization. However, CI (OX-PHOS) FCR showed significantly higher values in G2435R-RYR1 HOM compared with WT (p = 0.014), suggesting a greater reliance on complex I function because it forms a higher proportion of the ETS capacity. This effect appears to carry over to the CI1CII (OXPHOS) combined state, which is also significantly higher in G2435R-RYR1 HOM (p = 0.004). On the other hand, T4826I-RYR1 HET knock-ins showed significantly lower CI (OX-PHOS) and CI1CII (OXPHOS) FCR compared with G2435R-RYR1 HET and HOM knock-ins, which implies that different RYR1 variants can alter skeletal muscle mitochondrial function in different ways.

Oxidative stress and protein expression
ROS generation in fibers with exposure to hydrogen peroxide (H 2 O 2 ) was assessed using the CellROX Green reagent. The increase in the CellROX Green fluorescence signal intensity generated over 30 min after treatment with 1 mM H 2 O 2 was significantly different among genotypes (Fig. 3, A and B). The twoway mixed analysis of variance (ANOVA) showed a significant interaction between genotype and time on ROS generation, F (10, 480) = 13.05, p , 0.001. When compared with WT mice, G2435R-RYR1 HOM mice showed significant increase in ROS at each time point, and significant differences were seen in G2435R-RYR1 HET mice at 20 min and beyond (Fig. 3A). Area under the curve analysis showed that the area for WT fibers (10,6133 6 10,923) was significantly less than for either G2435R-RYR1 HET (15,4722 6 17,844, p = 0.028) or G2435R-RYR1 HOM (21,2481 6 11,617, p = 0.013) fibers. Comparisons between HET and HOM G2435R-RYR1 knock-ins also showed a significant difference (p = 0.008) related to gene dose (Fig. 3B).
Western blot analyses quantifying citrate synthase expression, a marker for mitochondrial content (Fig. 3C), showed no significant differences between genotypes at the protein level (WT versus HET p = 0.997, HET versus HOM p = 0.878, WT versus HOM p = 0.911). Based on their important roles in the ROS pathway we measured the expression of both NADPH-oxidase 2 (NOX2) and superoxide dismutase 1 (SOD1). NOX2 catalyzes the reduction of molecular oxygen into superoxide ions and SOD1 is responsible for their breakdown (18). NOX2 expression levels were significantly higher in G2435R-RYR1 HOM versus WT muscle (0.74 6 0.07 versus 0.44 6 0.07, p = 0.047) (Fig. 3D). In contrast, the opposite was seen with SOD1, which was significantly decreased in G2435R-RYR1 HOM compared with WT muscle (0.68 6 0.08 versus 0.33 6 0.08, p = 0.017, Fig. 3E). Interestingly, despite the difference in the slopes of the increase in CellROX Green intensity, neither NOX2 nor SOD1 protein expression in G2435R-RYR1 HET muscle was significantly different from WT (NOX2 WT versus HET p = 0.08, SOD1 WT versus HET p = 0.25).

Discussion
Transcriptome profiles between WT and RYR1 knock-in mice revealed differential gene expression at rest. Gene expression profiles between WT and G2435R-RYR1 HET mice show great similarity, whereas G2435R-RYR1 HOM and T4826I-RYR1 HET comparisons featured several hundred differentially expressed genes when compared with WT. This disparity demonstrates that RYR1 mutations can impact skeletal muscle gene expression in a variant-specific manner.
Pathway and gene ontology analysis found that a large proportion of down-regulated genes in G2435R-RYR1 HOM and T4826I-RYR1 HET knock-in muscles were localized in mitochondria, many of which were responsible for regulating fatty acid metabolism. Free fatty acids (FFA) are an important source of energy in muscle. They are metabolized by first undergoing b-oxidation to generate acetyl-CoA, which is a substrate for the Krebs cycle in the mitochondrial matrix (19). Impaired fatty acid oxidation and transport results in the accumulation of FFA in skeletal muscle, which causes lipotoxicity and is associated with insulin resistance, oxidative stress, and mitochondrial dysfunction (19)(20)(21). This appears to be a mutant RYR1-related phenotype because several of the genes related to this process (acyl-CoA dehydrogenase medium chain (Acadm), 2,4dienoyl-CoA reductase 1 (Decr1), and electron transfer flavoprotein subunit a (Etfa)) were shared between the WT versus G2435R-RYR1 HOM and WT versus T4826I-RYR1 HET comparisons. Furthermore, it was previously speculated that there is reduced fatty acid oxidation in the R163C-RYR1 knock-in mouse model because there is reduced expression of enzymes that regulate this process in R163C mitochondria (9,22).
The human orthologues of differentially expressed genes found in G2435R-RYR1 HOM (Acadm, acyl-CoA dehydrogenase very long chain (Acadvl), hydroxyacyl-CoA dehydrogenase trifunctional multienzyme complex subunit a (Hadha), Hadhb, carnitine palmitoyltransferase 2 (Cpt2), Etfa, and electron transfer flavoprotein dehydrogenase (Etfdh)) and T4826I-RYR1 HET mouse muscle (Acadm and Etfa) are associated with fatty acid oxidation disorders (23). Connections between FFA and MH pathology were first described in MH susceptible pigs after significantly higher levels of FFA and phospholipase activity were found in mitochondria (24)(25)(26). Several years later, elevated FFA production in MH-susceptible human muscle was linked to increased sensitivity to halothane (27)(28)(29). It was hypothesized that FFA and lipid peroxidation contributed to calcium dysregulation by inducing calcium release from the SR, Figure 3. G2435R-RYR1 mice show increased oxidative stress. A, This graph shows fluorescence signal changes with 5 mM CellROX Green fluorescent dye, which reflect ROS levels in FDB muscle fibers treated with 1 mM of H 2 O 2 for 30 min, normalized, and compared with WT mice (n = 41 fibers from four mice for WT, n = 28 fiber from four mice for G2435R-RYR1 HET, and n = 30 fibers from three mice for G2435R-RYR1 HOM in at least five independent experiments). Significant interactions between variables were assessed using a two-way mixed ANOVA before differences at each time point were analyzed using a one-way ANOVA. Significant post hoc comparisons to WT are annotated using asterisks *p , 0.05, **, p , 0.01, ***p , 0.001. B, comparison of area under curve from CellROX Green assay between G2435R-RYR1 HET, HOM, and WT fibers. C, D, and E, Western blot analyses showing Citrate Synthase, NOX2, and SOD1 protein bands in mice gastrocnemius muscle, respectively, with the bar charts representing band densitometry intensities normalized to GAPDH signals and averaged WT signals (n = 3 mice for each genotype). Data represent mean 6 S.E. These data were analyzed using a one-way ANOVA, and significant post hoc comparisons between genotypes are annotated with asterisks *p , 0.05, **p , 0.01, ***p , 0.001. with potential to modulate the MH phenotype (6,24,25,30). The results seen in our RNAseq analysis suggest that this line of research warrants further investigation. Although biliverdin reductase B (Blvrb) was the only transcript up-regulated in all three RYR1 knock-ins tested, it does not participate in any of the enriched pathways and gene ontology terms identified and is therefore unlikely to represent a common mechanistic link across the genotypes studied. Blvrb encodes a protein associated with heme metabolism and redox regulation (31, 32) and its up-regulation is consistent with a response to increased oxidative stress in all genotypes.
Other areas of interest regarding specific genes include the down-regulation of NADH:ubiquinone oxidoreductase core subunit S1 (Ndufs1) and Ndufs2 (encoding complex I), succinate dehydrogenase complex flavoprotein subunit A (Sdha) (encoding complex II), and ubiquinol-cytochrome C reductase core protein 2 (Uqcrc2) (encoding complex III), which were found to be differentially expressed in the WT versus G2435R-RYR1 HOM comparison. Reduced gene expression of mitochondrial complexes suggests multicomplex impairment of the electron transport chain because of an imbalance in turnover. Investigation into mitochondrial function using high resolution respirometry showed significant functional differences in G2435R-RYR1 HOM and T4826I-RYR1 HET genotypes. The results from two different markers for mitochondrial content were not consistent in our study, because citrate synthase protein expression showed no differences between genotype whereas the CIV (MAX) output suggested evidence of lower mitochondrial content in G2435R-RYR1 HOM muscle (33). However, we did not see this change in G2435R-RYR1 HET and T4826I-RYR1 HET muscles despite evidence in previous studies that the latter had lower mitochondrial content (10).
The LEAK FCR provide an indication of ETS uncoupling, reflecting the proportion of ETS capacity related to nonphosphorylating respiration (proton leakage). Comparison of LEAK FCR shows that mitochondria from both G2435R-RYR1 HET and HOM genotypes have normal ETS coupling efficiency compared with WT, whereas T4826I-RYR1 HETs have a lower LEAK FCR indicating higher ETS coupling efficiency. This was surprising because MH-susceptible muscle is often associated with mitochondrial uncoupling at rest (9,34,35). The elevated CI (OXPHOS) FCR observed in G2435R-RYR1 HOM mice suggests that they are more reliant on complex I function because it takes up a higher proportion of the ETS capacity. This is thought to be due to an increase in electron flow through complex I, which may be an adaptation to compensate for lower mitochondrial content and to maintain the energy demands of the muscle.
Compared with both G2435R-RYR1 genotypes, the electron flow through complex I in phosphorylating states was lower in T4826I-RYR1 HETs. WT versus T4826I-RYR1 HET comparisons were not significantly different, and this may have been because of compensation by the increased ETS coupling efficiency observed. In a previous study of T4826I-RYR1 mitochondrial function in myotubes using Seahorse ® technology, a reduction in basal respiration and ETS capacity was reported, but this was not observed in our findings (10). These contrasting results between studies is perhaps due to methodological differences and our use of adult muscle fibers as opposed to cultured myotubes. To assess the proportion of electron flow contributing to ATP production, we analyzed maximum OXPHOS capacity corrected for leak respiration (netOXPHOS), and results indicate a higher ATP production rate in G2435R-RYR1 HOMs (Fig. S1) Consistent with studies of other MH-susceptible RYR1 knock-in mice, we have shown that both HET and HOM G2435R-RYR1 muscle have elevated levels of oxidative stress (8,9,12,36). Two major oxidative stress proteins, NOX2 and SOD1, were quantified, and their differences in expression suggest that there is an elevated baseline level of superoxide ions in the G2435R-RYR1 muscles that is the result of both increased superoxide production and a reduced ability to metabolize it. In addition, the increased FCR in complex I of the G2435R-RYR1 HOM mice may lead to increased superoxide production because it is one of the main sources of superoxide in the mitochondria (37,38). Superoxide will likely damage cells and impair function, which may explain the deficient and swollen mitochondria described in other RYR1 knock-in MH mouse models (5,8). Our study also shows that G2435R-RYR1 HOM and HET mice have a significant and gene-dose dependent elevation of CellROX Green-detected ROS production in response to H 2 O 2 -induced oxidative stress compared with WT mice. The increase in ROS level has been attributed to an elevation in superoxide and/or hydroxyl radicals, because CellROX Green is found to be reactive to both species but not to H 2 O 2 (39)(40)(41). A decrease in antioxidant capacity in combination with an increase in ROS production could be responsible for the higher levels of oxidative stress in G2435R-RYR1 mice (37,42).
The lack of differences seen in gene expression and mitochondrial function between WT and G2435R-RYR1 HET mice highlights the complexity of the MH phenotype. Individuals with the equivalent c.7300G.A human MH variant are generally associated with weaker clinical phenotypes, and our observations with G2435R-RYR1 HET mice reflect this (4). The metabolic phenotype in MH is complicated, but most of the evidence suggests bioenergetic impairment in susceptible individuals. Here we have demonstrated evidence of mitochondrial dysfunction and increased oxidative stress in the G2435R-RYR1 knock-in mouse model. These outcomes have likely been perpetuated by elevated [Ca 21 ] rest and possible deficiencies in fatty acid metabolism, both of which are heavily involved in regulating muscle bioenergetics. Given its complex nature, future research exploring fatty acid metabolism and glycolysis would be ideal for a more comprehensive view of how RYR1 mutations affect metabolism in skeletal muscle.

Experimental procedures
Animals Animals were housed in pathogen-free conditions with free access to food, water, and 12-h light-and-dark cycles. All experiments were undertaken with United Kingdom Home Office approval. For transcriptomics, 8-week-old male mice were used, and a mixture of male and female mice between 8-14 weeks old were used for functional studies.

Preparation of samples for RNA-Seq
WT mice for sequencing were created within the G2435R-RYR1 colony by breeding HET mice together. A total of 12 mice of four genotypes, WT, G2435R-RYR1 HET, G2435R-RYR1 HOM, and T4826I-RYR1 HET knock-ins (age 8 weeks, n = 3 for each genotype), were killed by cervical dislocation. Hind limbs were removed and kept in oxygenated Krebs ringer solution for 30 min at 37°C before soleus muscles were dissected and transferred into RNAlater (Invitrogen) at 280°C until further processing.
The stored soleus muscles were homogenized in 1 ml of TRI-zol©, and RNA was extracted using a chloroform/isopropanol extraction method, followed by an RNA clean-up step using the RNAeasy microkit (Qiagen) with on-column DNase treatment, according to manufacturer's instructions. The Agilent 4200 Tapestation was then used to assess RNA integrity. RNA samples with integrity scores of RINe . 7.0 were considered acceptable for sequencing.
Truseq ® stranded mRNA library preparation Poly(A)-selected mRNA libraries were created from the isolated RNA by the Next Generation Sequencing Facility at St. James's Hospital, Leeds. RNA samples were first quantified using the Qubit TM RNA BR Assay Kit (Invitrogen) before being used to create libraries using the Truseq Stranded mRNA library preparation kit (Illumina), according to manufacturer's guidelines. All cDNA libraries were quality-checked using the Agilent D1000 screentape and quantified using the Quant-iT TM PicoGreen ® dsDNA assay.

Illumina HiSeq ® next-generation sequencing
Twelve mRNA-enriched libraries were indexed and compiled into a single pool at equimolar concentrations. This library pool was then distributed across two HiSeq 3000 150-bp paired-end lanes, achieving an average of 25 million reads per sample. A forward and reverse FASTQ file was generated for each sample on each lane.

Differential gene expression analysis
Raw sequencing data (FASTQ files) were initially processed using the Medical Advanced Research Computer, a high-performance computer cluster at the University of Leeds. The FASTQ files for each technical replicate across HiSeq lanes were combined and trimmed using Cutadapt to remove adapters and bad-quality base calls (,10) (43). Quality assessment before and after trimming was performed using FastQC (RRID: SCR_014583) (44). All samples were aligned using the STAR aligner and quantified using featureCounts (45,46). Processed sequencing data were then imported into RStudio and analyzed using the DESeq2 package for differential gene expression (47,48).

Pathway and gene ontology
Lists of differentially expressed genes, generated from DESeq2, were functionally characterized to identify enriched pathway terms. This was achieved using the online enrichment analysis tool Enrichr (49,50). Gene lists were compared with 176 pathway terms in the "WikiPathways_2019_Mouse" geneset library with a gene coverage of 4558 genes.

High-resolution respirometry
Oxygen consumption over time was measured using Oroboros respiratory analyzers (Oroboros Instruments, Innsbruk, Austria) at 2-s intervals with polarographic oxygen sensors and expressed as mass-specific oxygen flux (pmol/s*mg). The analyzer was calibrated daily in air-saturated solution before experimentation. Assays were initiated by injecting oxygen into each chamber to raise the oxygen concentration to .400 nmol/ml before starting a substrate-uncoupler-inhibitor titration (SUIT) protocol. Re-oxygenation of the chambers was performed to maintain oxygen concentration between 200-500 nmol/ml to prevent limitation due to oxygen diffusion (51). Each assay was performed at 37°C, with chamber stirrers set at 750 rpm.

CIV (MAX)
Measurement of mass-specific complex IV oxygen flux was used as an alternative proxy marker for estimations of mitochondrial content. Ascorbate (2 mM) and N,N,N9,N9-tetramethyl-p-phenylenediamine (0.5 mM) were applied to each sample at the end of the standard SUIT protocol (33). The measurement of CIV (MAX) is taken at the peak of the corresponding trace, and the residual chemical background is finally removed from this value by applying sodium azide (10 mM) to the sample.

Adult single fiber preparation
Single fibers were isolated from flexor digitorum brevis (FDB) muscles using collagenase digestion (53). FDB muscles were incubated at 37°C, 5% CO 2 in culture medium containing 0.2% (w/v) collagenase type I (Abnova) and 0.4% BSA (Sigma-Aldrich) for 2.5 h. Following incubation, the fibers were separated by gentle trituration in 3 ml of collagenase-free culture DMEM (Gibco) containing 100 units/ml of penicillin and 100 mg/ml of streptomycin (Gibco) and 10% heat-inactivated FBS (Gibco). The fibers were then allowed 30 min to settle on 96well plates coated with Cultrex basement membrane matrix (Trevigen).

CellROX green ROS assay
CellROX Green dye 5 mM (Thermo Fisher Scientific) dissolved in mammalian Ringer solution (133 mM NaCl, 5 mM KCl, 1 mM MgSO 4 , 10 mM HEPES, 5.5 mM glucose, and 2 mM CaCl 2 ) was applied to the dissociated fibers for 30 min at 37°C, 5% CO 2 . The fibers were washed twice with Ringer solution and then transferred to the stage of a Nikon TE2000 inverted microscope. H 2 O 2 1 mM dissolved in Ringer solution was perfused onto the fibers at ;100 ml/s for 30 s. CellROX Green was illuminated at 485 nm with the X-Cite 120 metal halide light source (EXFO, Ontario, Canada) and fluorescence emission was measured at 520 nm using a 103, 0.3-numerical aperture objective. The data were collected with an intensified 10-bit digital intensified charge-coupled device at 1 fps (Stanford Pho-tonics) from 3-6 individual fibers and analyzed using Piper image analysis software (Stanford Photonics).

Data handling and statistical analysis
Gene expression data were analyzed with the DESeq2 package, which performs significance testing using a Wald test. This generates lists of differentially expressed genes with p values adjusted for multiple testing using the Benjamini-Hochberg adjustment (47,54). Pathway and gene ontology analysis of gene lists was done using Enrichr, which calculates enrichment scores using statistics based on the Fisher's exact test (49,50).
Respirometry data are presented as oxygen flux/muscle mass (pmol/s*mg) or are normalized as flux control ratios, which excludes the confounding effects of both mitochondrial quality and quantity. Oxygen flux readings were exported from specialized software DatLab5 (Oroboros Instruments) into Microsoft Excel before analysis using Statistical Product and Service Solutions (SPSS) statistics analysis software (version 25.0). Respirometry data were not normally distributed when assessed by the Shapiro Wilk test (p , 0.05) and were subsequently analyzed using a Kruskal-Wallis H test with a post hoc Dunn's test for pairwise comparisons. Data were summarized as box plots using GraphPad Prism 8.
CellROX Green ROS assay data are presented as mean fluorescence changes from baseline, normalized to WT fluorescence (F 1 -F 0 ). A two-way mixed ANOVA was used to investigate interaction effects between time and genotype on ROS production using SPSS statistics analysis software (version 25.0). One-way ANOVAs were used to compare the ROS differences between genotypes at each time point, and pairwise comparisons were made using Tukey post hoc tests. Area under the curve calculations between WT, G2435R-RYR1 HET, and HOM mice were also compared using a one-way ANOVA with Tukey post hoc tests. Western blotting band densities were obtained using ImageJ software, and expression levels of the protein(s) of interest were normalized to the intensity of the GAPDH signal in the same lane. The normalized band densities for citrate synthase, SOD1, and NOX2 were compared between WT, G2435R-RYR1 HET, and HOM mice using a one-way ANOVA with Tukey post hoc tests using SPSS statistics analysis software (version 25.0).

Data availability
All data are contained within this article or in the supporting information provided. All sequencing data files have been uploaded to the NCBI Sequence Read Archive accession: PRJNA594852, release date 2020-08-20. Conflict of interest-The authors declare that they have no conflicts of interest with the contents of this article.