RESEARCH ARTICLE -Omic Approaches to Understanding Muscle Biology Proteome and transcriptome proﬁling of equine myoﬁbrillar myopathy identiﬁes diminished peroxiredoxin 6 and altered cysteine metabolic pathways

Proteome and transcriptome proﬁling of equine myoﬁbrillar myopathy identiﬁes diminished peroxiredoxin 6 and altered cysteine metabolic pathways. Physiol Genomics 50: 1036–1050, 2018. myoﬁbrillar myopathy (MFM) causes exertional muscle pain and is characterized by myoﬁbrillar disarray and ectopic desmin aggregates of unknown origin. To investigate the pathophysiology of MFM, we compared resting and 3 h postexercise transcriptomes of gluteal muscle and the resting skeletal muscle proteome of MFM and control Arabian horses with RNA sequencing and isobaric tags for relative and absolute quantitation analyses. Three hours after exercise, 191 genes were identiﬁed as differentially expressed (DE) in MFM vs. control muscle with (cid:2) 1 log 2 fold change (FC) in genes involved in sulfur compound/cysteine metabolism such as cystathionine-beta-synthase ( CBS , 2 4.51), a cysteine and neutral amino acid membrane transporter ( SLC7A10 , 2 1.80 MFM), and a cationic transporter (SLC24A1, 2 1.11 MFM). In MFM vs. control at rest, 284 genes were DE with (cid:2) 1 log 2 FC in pathways for structure morphogenesis, ﬁber organization, tissue development, and cell differentiation including (cid:2) 1 log 2 FC in cardiac alpha actin ( ACTC1 1 2.5 MFM), cytoskeletal desmoplakin ( DSP 1 2.4 MFM), and basement membrane usherin ( USH2A 2 2.9 MFM). Proteome analysis revealed signiﬁcantly lower antioxidant peroxiredoxin 6 content (PRDX6, 2 4.14 log 2 FC MFM), higher fatty acid transport enzyme carnitine palmitoyl transferase (CPT1B, 1 3.49 MFM), and lower sarcomere protein tropomyosin (TPM2, 2 3.24 MFM) in MFM vs. control muscle at rest. We propose that in MFM horses, altered cysteine metabolism and a deﬁciency of cysteine-containing antioxidants combined with a high capacity to oxidize fatty acids and generate ROS during aerobic exercise causes chronic oxidation and aggregation of key proteins such as desmin. The objectives of the present study were to compare the skeletal muscle transcriptome and proteome of healthy and MFM affected Arabian horses at rest, to compare the alterations in gene expression that occurred with aerobic exercise and to compare the transcriptome of MFM and control horses 3 h after exercise. Procedures were approved by the Oregon State University Institutional Animal Care and Use Committee #4480, and horse owners provided informed consent.


BACKGROUND
Arabian horses are adept at endurance races, routinely competing in distances of 50 -100 miles. Exertional rhabdomyolysis (ER), literally the dissolution of skeletal muscle with exercise, is a common occurrence in endurance horses, affect-ing 4 -18% of competitors (68). While overexertion and environmental factors can cause ER, a genetic predisposition to ER is a common occurrence in some breeds (61). To date, however, genetic causes of ER have not been identified in the Arabian breed. We recently identified a muscle disease with a suspected genetic predisposition in Arabian horses called myofibrillar myopathy (MFM) that causes muscle stiffness and intermittent ER (45,63).
Across species, MFM comprises a group of muscle diseases defined by morphological characteristics of Z-disc disarray, disruption of myofiber alignment, and ectopic accumulation of proteins that contain a variety of cytoskeletal and Z-disc material, depending upon the causative genetic mutation (56). The extent of ectopic protein aggregation, degree of change in mitochondrial staining, and severity of muscle weakness in equine MFM is frequently less than that described for human patients with MFM (53,56,63). Approximately half of all MFM cases in human medicine are caused by mutations in genes encoding sarcomeric and extrasarcomeric proteins, including desmin, filamin C, plectin, ZASP, myotilin, ␣B-crystallin, and BAG3, while the remaining diseases are due to still unresolved gene defects (20,34) Although common in humans, Arabian horses with MFM do not appear to have impaired cardiac function as they have no clinical indicators of heart failure. Even when over 20 yr of age, MFM horses complete endurance rides without abnormal elevations in heart rates that are routinely evaluated at compulsory veterinary evaluations during and after endurance races (45). The etiopathology of MFM in horses and the protein(s) or gene(s) responsible remain unknown.
RNA-Seq and proteomic analyses have been of value in defining pathophysiological mechanisms for MFM and muscular dystrophies in human patients (24,27). Three studies have evaluated the equine exercise transcriptome in healthy Thoroughbred or Arabian racehorses performing nearmaximal to maximal exercise (7,43,52). To date, no studies have investigated the gluteal muscle equine transcriptome or proteome after aerobic exercise. We hypothesized that transcriptome and proteome analyses of skeletal muscle before and several hours after aerobic exercise would provide critical insights into the pathophysiology of MFM in horses.

Aim
The objectives of the present study were to compare the skeletal muscle transcriptome and proteome of healthy and MFM affected Arabian horses at rest, to compare the alterations in gene expression that occurred with aerobic exercise and to compare the transcriptome of MFM and control horses 3 h after exercise. Procedures were approved by the Oregon State University Institutional Animal Care and Use Committee #4480, and horse owners provided informed consent.

Horses
Horses for the study were recruited by contacting a large group of active endurance riders in the Pacific Northwestern US and selecting four Arabian and one Arabian-cross horse (age 15.8 Ϯ 7.1 yr, 3 castrated males, 2 females) that had clear evidence of MFM based on clinical history of an exertional myopathy and histopathologic evaluation of muscle samples including desmin immunohistochemistry. Six control Arabian horses (13.0 Ϯ 6.2 yr, 4 castrated males, 2 females) that were located on the same property and had similar feeding and exercise routines and no apparent histopathology in muscle biopsies were used as a matched control. All MFM horses had abnormal desmin-positive aggregates in gluteal muscle fibers (data reported previously) (63). ER was previously documented in 5 MFM horses, with the last episode at least 6 mo before taking part in the study. Control horses had participated in endurance activities for a minimum of 3 yr with no prior evidence of ER and had no evidence of muscle pathology in gluteal muscle biopsies. The muscle biopsy specimens utilized in the current experiment were aliquots of those obtained for previously published studies (45,63).

Standardized Exercise Test
An aerobic exercise test was used that all horses could complete regardless of their current fitness. The details of the test and the metabolic responses of the MFM and control horses included in the present study have been published previously (45). In brief, horses were stall rested for 24 -48 h, and then a rider with a global positioning system-enabled watch guided the horse through a 47 min exercise test consisting of precise intervals of walk, trot, and canter over 6.4 -7.2 km (4.0 -4.5 miles). One owner of one control horse elected not to perform the exercise test, leaving five MFM and five control horses to complete the test with no differences in heart rate response, blood lactate, blood glucose, or plasma electrolyte concentrations identified between MFM and control horses (45). No horses displayed clinical signs of muscle stiffness during or after the exercise test. Plasma creatine kinase activities were within normal limits before and 3 h following exercise in all but one MFM horse that had a mild ϳ2-fold elevation in creatine kinase activity (1116 U/l; normal range 145-633 U/l) following exercise (45).

Muscle Samples
A 6 mm Bergstrom needle biopsy of the middle gluteal muscle was obtained at a depth of~6 cm at a standardized site 17 cm along a line running from the dorsal tuber coxa to the tail head before and 3 h after exercise ceased on alternate sides of the gluteal muscle (63). Three hours postexercise was chosen because in a previous equine study, little differential expression (DE) of genes was noted immediately after exercise, whereas significant changes occurred 4 h after exercise (43). Samples for transcriptome analysis were immediately frozen in liquid nitrogen and stored at Ϫ80°C until analysis.

Muscle Histochemistry
Samples for muscle histopathology obtained at rest were affixed to cork squares with optimal cutting temperature media and frozen in isopentane suspended in liquid nitrogen. Muscle fiber types were determined with myosin ATPase staining on 7 m thick cryostat sections with 5 min preincubation at pH 4.4 (13). Muscle fiber type composition for type 1, 2A, and 2X fibers was determined by counting a minimum of 300 muscle fibers per horse. A one-way ANOVA and Tukey post hoc test was used to compare fiber type proportions between MFM and controls. Hematoxylin and eosin (HE) staining was performed on 7 m thick cryostat sections.

Immunohistochemistry
Immunohistochemistry was used to investigate desmin aggregation, regeneration, and myosin isoforms that had altered gene expression as well as peptide expression in the proteomic profile. Staining was performed for desmin (1:100 anti-desmin-mouse monoclonal D33; Agilent Technologies, Santa Clara CA) and MHCd (1:50 mouse monoclonal anti-MHCd; Leica Biosystems, Buffalo Grove, IL) as previously described (63,64). For MYH13 staining, flash-frozen muscle samples 7 m thick were placed on charged slides and air-dried overnight at 25°C. Endogenous enzyme blocking was performed in 0.3% hydrogen peroxide/Tris-buffered saline (TBS). After rinsing in tap water, pH adjustment was made with TBS, pH 7.4 (Scytek Laboratories, Logan, UT). Subsequently, standard micropolymer complex staining steps were performed at room temperature on the IntelliPath Flex Autostainer. All staining steps were followed by rinses in TBS Autowash buffer (Biocare Medical, Concord, CA). After blocking for nonspecific protein with Background Punisher for 10 min (Biocare), sections were incubated in rabbit polyclonal anti-MYH13 (PA5-70713; Thermo Fisher Scientific, Rockford, IL) diluted 1:50 in normal antibody diluent (Scytek). After rinsing, Promark Rabbit on Farma horseradish peroxidase micropolymer (Biocare) was applied for 30 min. Immunoreactivity was detected by using Romulin AEC for 5 min (Biocare). Slides were counterstained with CATHE hematoxylin.

Transcriptomics
RNA extraction and sequencing. Total muscle RNA was isolated with TRIzol/chloroform extraction and a biopulverizer (BioSpec Products, Fartlesville, OK) as previously described (22). DNase treatments were performed in columns (Direct-zol RNA MiniPrep Plus; Zymo, Irvine, CA) with TURBO DNase (Thermofisher, Wilmington, DE) according to manufacturer's instructions. Quantification and quality of RNA, along with degree of rRNA contamination, were assessed with the Pico chip on the Agilent Bioanalyzer 2100 (Santa Clara, CA), and samples with an RNA integrity number Ն7 were used for library preparation and quantification. Library construction with a strand-specific polyA ϩ capture protocol (TruSeq Stranded RNA, Illumina) was performed and sequencing completed using the Illumina HiSeq 2000 genome analyzer (100 bp PE) at a targeted 35 million reads/sample.
Mapping and assembling. All paired end RNA-Seq reads were initially assessed for quality with FASTQC (1). Samples that passed through the quality threshold of 30 (Q Ն 30) were aligned with the STAR aligner (16) to Ensembl horse CDS (EquCab2.86). A total of 26,993 sequences were used to create the index file and map RNA-Seq reads. The final number of reads aligned ranged from 31 to 38 M per sample, with uniquely mapped reads falling in the range of 87.92-94.12% (Supplemental Table S1). (The online version of this article contains supplemental material.) Count data for each sample were generated from the STAR-aligned BAM files using the internal flag in STAR. Differential gene expression (DE) was identified with EdgeR (51) generalized linear models (GLMs) available through R/Bioconductor (R Core Team, 2016), making assumptions to shrink genespecific variance toward an underlying common variance for all genes (41,51). The number of counts per million (CPM) were set to include a minimum of five horse samples. Pairwise differences among means and linear combinations of model parameters were used to evaluate the DE between control and MFM affected horses. The fixed effects model can be represented as: ϩ ␣i ϩ ␤ j ϩ ␣␤ij, where is the overall mean, ␣i represents effects due to any confounding effects,␤ j represents effects due to specific treatment, and ␣␤ij represents interaction among effects. Such models are also flexible to incorporate nested effects, and specifically for these data we represent the model as: ϩ ␣ treatment(horse), where is the overall mean, and ␣treatment(horse) is the treatment effect nested within horse to accommodate repeated measures over each individual horse. There were no confounding covariates included in our model. All gene lists derived from pairwise comparisons were false discovery rate (FDR) corrected for multiple testing at a threshold of P Ͻ 0.05.
Enrichment and pathway analysis of transcriptomics. Significant gene lists (P Ͻ 0.05) derived from pairwise comparison of transcriptomics data were analyzed with Cluego (2.3.2) (6), a Cytoscape (version 3.4.0) framework-based plug-in to visualize biological terms and large gene data sets and perform pathway analysis. The ENSEMBL gene IDs were mapped to the Gene Ontology (GO) database using all evidence without electronic annotation (IEA), which excludes computational annotations. Data were assessed at the biological level in GO. Default enrichment/depletion analysis was performed with a suitable background comprising~11,000 genes, and P values were adjusted for multiple testing by Bonferroni step-down, which is the default in Cluego. The final list of genes evaluated in this study was derived from those genes that had enriched GO terms. Further information regarding the function of these genes was obtained from gencards.org.

Transcriptome Profiling
After using the CPM feature to ensure that each transcript was present across at least five horses, we reduced the initial 26,993 transcripts identified in Ensembl to 11,330 transcripts, constituting the final transcripts used for identification of DE genes. Using the resulting list of genes from the final transcripts, we made pair-wise comparisons (pre-vs. postexercise, MFM vs. control at rest, MFM vs. control postexercise, MFM postexercise vs. control pre-exercise) with a GLM model. FDR-corrected significant genes (P Ͻ 0.05) were used for GO term analysis.

Quantitative Real-time Polymerase Chain Reaction
Genes with Ͼlog 2 fold higher or lower expression in MFM vs. control muscle were also evaluated by quantitative real-time polymerase chain reaction (qRT-PCR). Primer design is provided in Table 1. Thermocycling was conducted with EvaGreen dye (Biotium, Fremont, CA), ROX Reference Dye (Invitrogen, Life Technologies), and Hot Start Taq DNA Polymerase (New England BioLabs, Ipswich MA), using the QuantStudio 3 Real-Time PCR System (ThermoFisher Scientific, Rockford IL). PCR reactions were run in duplicate (20 l volume reactions). Each reaction contained 2 l of sample cDNA, 2 l of 2.5 mM dNTPs, 2 l of 10ϫ PCR buffer, 1 l of EvaGreen dye, 1.5 l of 1:10 ROX Reference Dye dilution, 0.125 l of Hot Start Taq DNA Polymerase, 2 l of 1.6 M forward primer, 2 l of 1,600 M reverse primer, and 7.4 l of sterile nuclease free distilled water. Reactions were run for 40 cycles under the following conditions: denaturation at 95°C for 10 min, annealing at 60°C for 1 min; melt curve stages at 95°C for 15 s, 60°C for 1 min, and 95°C for 15 s. Cycle thresholds (CT) were automatically calculated by the QuantStudio 3 Real-Time PCR System. For each gene of interest 100% geometric efficiency was established. Nontemplate controls run for each gene showed no amplification. Mann-Whitney test was used to determine statistical significance with P Ͻ 0.05.

Coding Sequence of Desmin
Potential mutations in the coding sequences of desmin were investigated as a putative cause of MFM by using initial raw fastq files realigned using the STAR two-pass alignment method (13), and the resulting BAM files were sorted; duplicates were marked (http:// broadinstitute.github.io/picard/), realigned, and recalibrated with GATK (44). The resulting BAM files were then used as an input to call variants with SAMtools (38). The resulting VCF file was filtered by BCFtools (46), and a final list of single nucleotide polymorphisms (SNPs) was compared with a case vs. control setup using SnpSift (10). The resulting SNP data were also annotated for variant function (synonymous or nonsynonymous).

Proteomics
Sample preparation, peptide fractionation, and mass spectrometry. We reconstituted 20 mg samples of frozen ground tissue from the same gluteal muscle biopsy obtained at rest from a subset of horses used for RNA-Seq For in-solution proteolytic digestion and isobaric tags for relative and absolute quantitation (iTRAQ) labeling, a 100 g aliquot of each sample was transferred to a new 1.5 ml microfuge tube and brought to the same volume with protein extraction buffer and 8 mM MMTS. All samples were diluted fourfold with ultrapure water and trypsin (Promega, Madison, WI), added in a 1:35 ratio of trypsin to total protein. Samples were incubated overnight for 16 h at 37°C after which they were frozen at Ϫ80°C for 30 min and dried in a vacuum centrifuge. Each sample was then cleaned with a 4 ml Extract Clean C18 SPE Primers sets used for quantitative RT-PCR sequencing of cardiac alpha actin (ACTC1), desmoplakin (DSP), Usherin (USH2A), delta notch-like protein 1 (DLK1), cystathione beta synthase (CBS), and housekeeping gene glyceraldehyde 3 phosphate dehydrogenase (GAPDH). *Whether exon junction was on the forward or reverse primer. cartridge from Grace-Davidson (Deerfield, IL); eluates were vacuum dried and suspended in dissolution buffer (0.5 M TEAB, pH 8.5) to a final 2 g/l concentration. We labeled 40 g for each sample with iTRAQ 8-plex reagent (AB Sciex, Foster City, CA) per manufacturer's protocol. After labeling, the 40 g aliquot of each sample was multiplexed together into one 1.5 ml tube and vacuum-dried. The multiplexed sample was cleaned with a 3 CC Oasis MCX SPE cartridge (Waters, Milford, MA), and the eluate was dried in vacuo.
The iTRAQ-labeled samples were resuspended in buffer A (10 mM ammonium formate pH 10 in 98:2 water-acetonitrile) and fractionated offline by high-pH C18 reversed-phase (RP) chromatography (69). A MAGIC 2002 HPLC (Michrom BioResources, Auburn, CA) was used with a C18 Gemini-NX column, 150 mm ϫ 2 mm internal diameter, 5 M particle, 110 Å pore size (Phenomenex, Torrence, CA). Buffer A was 10 mM ammonium formate, pH 10 in 98:2 water-acetonitrile, and buffer B was 10 mM ammonium formate, pH 10 in 10:90 water-acetonitrile. The flow rate was 200 l/min with a gradient from 5 to 35% buffer B over 60 min, followed by 35-60% over 5 min. Fractions were collected every 2 min, and UV absorbances were monitored at 215 and 280 nm. Peptide containing fractions were divided into two equal numbered groups: "early" and "late." The first early fraction was concatenated with the first late fraction, and so on. Concatenated samples were dried in vacuo, resuspended in loading solvent (98:2:0.01, water-acetonitrile-formic acid), and 1-1.5 g aliquots were run on a Velos Orbitrap mass spectrometer (Thermo Fisher Scientific, Waltham, MA) as described previously (39) with the following modifications: lock mass was not used, HCD activation time was 20 ms, dynamic exclusion duration was 15 s, and the minimum signal threshold for data-dependent trigger was 20,000 counts. The mass spectrometer RAW data were converted to mzXML using MS Convert software from the ProteoWizard Toolkit (PMID: 23051804) and to MGF files using TINT raw-to-mgf converter (inhouse tool).
Protein identification and statistical analysis. MGF-formatted files were used for identification of peaks by using several search engines to limit potential pitfalls of using a single algorithm/method and to perform more exhaustive survey of the spectra (57). Multiple search engines such as OMSSA (version 2. Peptides and proteins were inferred from the spectrum identification results with PeptideShaker (version 1.16.4) (65). Peptide spectrum matches (PSMs), peptides, and proteins were validated at a 1.0% FDR estimated from the decoy-hit distribution. The resulting file was imported into Reporter package (version 0.7.2) (http://compomics.github.io/projects/reporter.html), which was included as complementary tool within the SearchGUI platform. Peptides that could belong to different members of the same protein family were handled by using stringent filtering and the Occam's razor approach. The three rules used by this approach were 1) For two proteins to be clustered, the sum of the probabilities of their shared peptides must be at least 95%.
2) The proteins must share at least 50% of their evidence. This is determined by summing the probabilities of the shared peptides and comparing this value with the summed probabilities of all of the peptides for each individual protein. If the sum of the probabilities of the shared peptides is greater than or equal to half of the sum of the peptide probabilities for either of the individual proteins, a cluster is formed. 3) A protein may be included in an existing cluster if it meets the above criteria with a member protein of the cluster. For two proteins to be clustered in a group, the sum of the probabilities of their shared peptides must be at least 95%.
To normalize data to controls, ratios were created per protein by dividing mean MFM values by mean control values. The resulting data table was exported as a text file, log transformed, and analyzed to determine significance between treatments with Perseus (62). P values were adjusted for multiple testing by the permutation-based FDR method included within Perseus. Individual spectra for all proteins with Ͼlog 2 fold change were manually checked for each sample/ animal to ensure there were no missing data. The biological significance of the final list of confident proteins was assessed with the Panther protein classification tool (59).

Transcriptomic vs. Proteomic Analyses
To find intersecting IDs between transcriptomics and proteomic data sets, we first converted all confident proteins (673 protein IDs) expressed in MFM and control horses to gene names with the ENSEMBL BioMart resource. Significant gene IDs (P Ͻ 0.05) from MFM vs. Control transcriptomics at rest (284 genes) and MFM vs. Control transcriptomics after exercise (191 genes) were evaluated for genes that had overlap with Gene IDs identified from confident proteins in the proteomics data. Because long-term changes in gene expression are required to alter protein expression, proteomic analysis was performed only on at-rest and not on postexercise samples (24).

Muscle Fiber Type Composition and Histopathology
There was no significant difference in skeletal muscle fiber type composition between MFM and control horses (Fig. 1, A and B). Fiber type composition was type 1 (MFM 17 Ϯ 3%; control 17 Ϯ 3%), type 2A (MFM 56 Ϯ 7%; control 52 Ϯ 8%), and type 2X fibers (MFM 27 Ϯ 4%; control 31 Ϯ 8%). Neither acute myofiber necrosis nor macrophage infiltration was evident in HE stains. Mature myofibers with small internalized myonuclei were present in 5/5 MFM horses and in 2/6 control horses. Increased endomysial or perimysial connective tissue was not evident in all MFM horses; however, two of the most severely affected MFM horses appeared to have a slight increase in endomysial connective tissue in some regions (Fig. 1,  C and D).
Immunohistochemistry. Immunohistochemistry identified desmin-positive aggregates in scattered myofibers of MFM but not control horses (Fig. 2, A and B). Neither small basophilic fibers with centrally located nuclei nor fibers staining positively for developmental myosin MYHd typical of regeneration were observed in any horses (Fig. 2, C and D). As a validation step for proteomic expression, immunohistochemical staining was performed for MYH13. Positive staining for MYH13 was found in control samples of extraocular muscle fibers and was also found in the walls of small blood vessels in MFM and control samples (Fig. 2, E-G). Only one desmin-positive skeletal muscle fiber in one MFM horse had mild focal MYH13 staining (Fig. 2E). Expression of MYH13 in proteomic data could reflect overlap in amino acid sequence with vascular smooth muscle myosin.

Transcriptomics
Control horses: effect of aerobic exercise. The comparison of transcriptomes of control horses at rest vs. after aerobic exercise identified 330 DE genes (adj. P Ͻ 0.05) (Supplemental Table  S1). After GO enrichment analysis, 12 genes had significant GO terms in biological processes involving corticosteroid rec-    Table S1). These genes had different GO terms compared with rest vs. after exercise in control horses (Fig. 3B). After enrichment analysis, 27 genes had significant GO terms in biological processes that involved regulation of intracellular signal transduction (GO1902531) and response to stimulus (GO:0048583), small GTPase signaling (GO:0007264), and cell differentiation (GO:0030154) (Fig.  3B). When the additional stringent filter of Ͼ1.0 log 2 fold  Table 2).

MFM Horses Postexercise vs. Resting Control Horses
This analysis was performed to compare the extremes: DE expression after exercise in MFM horses to a healthy resting horse. Significant DE of 883 genes was found when comparing MFM muscle after exercise to resting control muscle. This encompassed 305 genes with 22 significant GO terms (Table  3), many involving inflammation/immune response and cell activation (Fig. 3E). Four sarcomeric genes had Ͼ1.0 log 2 higher expression in MFM (ACTC1, 12.75 MFM): ankrin repeat domain 1 (ANKRD1, 11.50 MFM), a myofibrillar stretch sensor and regulator of lipid metabolism; myosin heavy chain 13 (MYH13 11.63 MFM, P ϭ 1.17E-33); myosin heavy chain 3 (MYH3 11.27 MFM, P ϭ 1.19E-06); as well as cytoskeletal gene DSP (11.36 MFM). Other genes with Ͼ1.0 log 2 altered expression that overlapped   Table S1).

qRT-PCR Validation of DE Transcripts
The transcripts with Ͼ2 log 2 fold change in expression in resting RNA-Seq data (MFM vs. control) included ACTC1 (1log 2 FC 2.5), DSP (1log 2 FC 2.4), DLK1 (2log 2 FC 2.4), and USH2A (2log 2 fold 2.9). A similar pattern of increased or decreased expression was found for these genes by qRT-PCR (Table 4). A significant increase in expression was found for DLK1 (P ϭ 0.004), and USH2A approached statistical significance (P ϭ 0.052) with the small number of samples in the comparison (5 MFM vs. 6 control). Additionally, lower CBS expression in RNA-Seq data for MFM vs. control after exercise (2log 2 FC 4.5) was also found with qRT-PCR, and this difference approached statistical significance (P ϭ 0.056) (Table 4).

Desmin Coding Sequence from RNA-Seq
Thirty-eight variants existed in the coding sequence of desmin between MFM and control horses; however, not one variant was consistently present in each MFM horse and consistently absent in control horses. There was one synonymous variant (Asp89Asp) that was present in five MFM as well as five control horses. The other identified variants were intronic (n ϭ 37), of which nine were potential splice region variants that did not segregate exclusively in MFM horses.

Proteomics
Confident protein identification. Using the target/decoy strategy on MFM and Control samples, we found 28,232 proteins to match our proteome spectra. SearchGUI identified 11,876 of these proteins with coverages ranging from Ͼ0.19 to 100%. The default setup of FDR (1%) and FNR (Ͻ5%) resulted in 990 proteins passing threshold. We identified 21 proteins that had duplicate ID; however, they were not included in further analyses because none of them passed threshold. There were no missing data with significant proteins. A further filtering of proteins in Peptide Shaker using default confidence scores resulted in a final list of 673 proteins that were classified as "confident." Protein classification. The 673 confident proteins were classified by the Panther protein class tool into 23 different protein classes of which oxidoreductase, cytoskeletal proteins, and hydrolase had over 40 incorporated proteins each (Fig. 4A) (60). Within this classification, there were 101 mitochondrial proteins (24 were Complex I) and 15 oxidoreductase proteins, including one catalase, four PRDX (1, 2, 5, 6), two thioredoxin, four glutathione, and four superoxide dismutase proteins. GO Biological pathway analysis of the confident proteins identified two primary biological pathways that involved either metabolic or cellular processes each containing Ͼ175 proteins. There were 10 other biological pathways that contained 75 or fewer proteins in MFM and control horse muscle (Fig. 4B). GO Molecular pathway analysis identified two major molecular pathways that involved catalytic activity or binding activity that contained Ͼ100 proteins and six other pathways that contained Ͻ40 proteins each (Fig. 4C).
DE of proteins. Statistical analysis identified three significantly DE proteins between MFM and control horses after controlling for multiple comparisons. The individual spectra were manually checked for each sample/animal to ensure there were no missing data for these three proteins. The three proteins were an antioxidant peroxiredoxin 6 (PRDX6, 2log 2 4.14 MFM); an enzyme that transports fatty acids across mitochondrial membranes, carnitine palmitoyl transferase (CPT1B 1log 2 3.49 MFM); and sarcomeric protein tropomyosin beta chain (TPM2 2log 2 3.24 MFM) (Fig. 5A, Table 5). Biological processes were evaluated through the additional inclusion of proteins with Ͼ2 log 2 fold change (11 proteins, adjusted P values between P ϭ 0.113 and 0.182) (Fig.  5A, Table 5). The major processes included thiol-specific antioxidants, oxidative metabolism, supramolecular fiber organization, proteasome degradation, apoptosis, immune function, and ribosomal/purine synthesis.

Overlap between Significantly DE Genes and Identified Proteins
Out of the total 673 confident proteins, 543 had encoding gene IDs present in both transcriptomic and proteomic data sets. These confident proteins were identified but not DE in MFM and control horses. When compar-ing the proteomic data set with DE genes in resting MFM vs. control muscle, we found seven significantly DE genes that had proteins identified by proteomic analysis (Fig. 5B). These included adenylhomocysteinase (LOC100054381), two proteins involved in ketone body metabolism (OXCT1, PGFS), and extraocular muscle myosin MYH13. A B C Fig. 4. Classification of pathways for the 673 confident proteins in MFM and control horse muscle obtained at rest. Proteins used in the analysis were identified by a Search GUI protein identification protocol and had Uni-ProtKB IDs. A: the largest number of proteins in equine muscle were classified as oxidoreductase, hydrolase, or cytoskeletal proteins by Panther analysis using Protein Plus tools. B: metabolic and cellular processes were the two major biological pathways for proteins in equine muscle according to GO biological pathway classification. C: catalytic activity and binding were the two major molecular pathways for proteins in equine muscle according to GO molecular pathway classification.
When we compared DE genes in postexercise MFM vs. control muscle with all the identified proteins in the resting muscle proteomic data set, two proteins, PDLIM1 (cytoskeletal protein) and HSPB7, which act to prevent protein misfolding, overlapped with the DE genes (Fig. 5B). There were five proteins identified in the proteomic data set of MFM and control horses that had encoding gene transcripts DE in MFM vs. control both at rest and after exercise (Fig. 5B). These included a cysteine endopeptidase (CAPN1), phospholipase C (PLCD4), a vesicle transport protein (NIPSNAP3B), and four proteins involved in the extracellular matrix (COL1A1, COL1A2, COL6A1, LAMB2). None of these five proteins were DE in MFM horses, and none of the DE genes had Ͼ1 log 2 fold DE.

Availability of Data
The data sets generated and/or analyzed during the current study are available in the National Center for Biotechnology Information's Gene Expression Omnibus (GEO) (17a) and are accessible through GEO Series accession number GSE104388. The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE (66) partner repository with the data set identifier "PXD009362".

DISCUSSION
The present study represents the first study of the combined proteome and exercise transcriptome in horses and the first to utilize these analyses to probe muscle disease pathogenesis in horses. MFM Arabian horses showed DE of genes and proteins involved in pathways of sarcomere and cytoskeletal structure, oxidoreductase activity, fatty acid and amino acid metabolism, proteosomal degradation, and apoptosis. These alterations in expression do not appear to reflect a general myopathic process, because microarray studies of equine polysaccharide storage myopathy (PSSM1) and equine recurrent exertional rhabdomyolysis (RER) do not show overlap with DE genes in the present study, particularly when evaluating the top 10 DE genes in these studies (4,5). Future RNA-Seq and proteomic analyses of other equine myopathies will help to elucidate which specific proteins and genes are DE as secondary effects of muscle disease and which are specific to MFM in Arabian horses.
Increases in connective tissue and muscle fiber regeneration are common features of many myopathies. It was interesting to note that four collagen and laminin proteins involved in the extracellular matrix had expressed peptides and DE genes in resting MFM vs. control muscle consistent with increased connective tissue. There was also a Ͼ2-fold log 2 increased expression of the basement membrane-associated gene USH2A in resting MFM vs. control muscle. Histologically, however, only a slight increase in endomysial connective tissue was apparent in two of the most severely affected MFM horses (most desmin aggregates). Regeneration was not apparent in MFM muscle based on HE and MHCd stains, and notably, DLK1 was significantly downregulated in MFM vs. control horses. Ablation of DLK1 in murine models results in impaired muscle regeneration (67).
A high proportion of high oxidative type 1 and 2A fibers in gluteal muscle is a feature of successful endurance racing Arabian horses and was found in the horses in the present study (49). This fiber type composition serves to delay fatigue by enhancing oxidation of free fatty acids during aerobic exercise and sparing muscle glycogen. Compared with control horses, MFM horses had a higher capacity to transport fatty acids into mitochondria for ␤-oxidation as indicated by the significantly higher content of carnitine palmitoyl transferase (13.5 log 2 fold MFM), the enzyme responsible for transport of long chain fatty acids into mitochondria. ␤-oxidation of free fatty acids during aerobic exercise though energy efficient, is also a primary generator of reactive oxygen species (ROS). ROS are generated primarily via mitochondrial Complex I and III of the electron transport chain and must be reduced by specific enzymes to prevent cellular damage. In MFM horses, the content of two out of 24 identified Complex I subunit proteins was Ͼ2 log 2 fold lower than in control horses (36). The potential to generate more ROS in MFM vs. control horses could arise from enhanced uptake of fatty acids into mitochondria in MFM horses, increased subsequent flux of electrons  through Complex I and increased ROS produced because of an altered Complex I subunit composition (Fig. 6). Notably, oxidoreductase activity was one of the primary pathways dysregulated in MFM in the present study. To compound the effect of ROS generation further, MFM horses appeared to have a lower capacity to neutralize ROS than controls, with a significantly lower content of the antioxidant PRDX6 (24 log 2 fold MFM) (Fig. 6). PRDX6 is a widely distributed cytoplasmic cytoprotective nonseleno peroxidase that reduces hydrogen peroxides and quite uniquely, phospholipid hydroperoxides, by using a cysteine-containing active site (23). Glutathione S transferase, a cysteine selenoprotein that reduces PRDX6, was among the DE genes in MFM vs. control muscle (20.5 log 2 fold change MFM) (23). While increased PRDX6 is described in a number of neurodegenerative diseases and some cancers, only one disease, alcoholic steatitis, reports decreased PRDX6 protein content with a much lesser decrement in PRDX6 than that found for MFM horses (1.6-fold vs. 18-fold) (47). A decrease in the amount of PRDX6 could arise from decreased synthesis or from irreversible cysteine oxidation to sulfinic (SO 2 H) and sulfonic acid (SO 3 H) with subsequent ubiquitination and degradation by the proteasome (42). Enhanced proteasomal degradation was suggested in proteomic data by 2.4 log 2 fold  6. Schematic diagram of interrelationship among selected identified proteins and genes differentially expressed in MFM vs. control horse muscle at rest and after exercise. MFM horses have significantly higher activity of CPTB1 (3.5 log2 fold), which transports fat into mitochondria and lower expression of 2/24 subunits of mitochondrial Complex 1, which is a site of free radical generation. The 4.14 log2 fold lower content of cysteine-dependent peroxiredoxin 6 (PRDX6) in MFM vs. control muscle at rest and lower expression of enzymes required to synthesize cysteine following exercise indicated that the oxidoreductase pathway plays a central role in the pathogenesis of MFM. We hypothesize that protein misfolding, myofibrillar disarray, and aggregation of desmin in MFM horses arise from protein oxidation due to excessive generation of reactive oxygen species (H2O2), decreased postexercise cysteine synthesis, and depletion of cysteine-dependent antioxidants such as PRDX6.
upregulation of the 26S proteasome subunit in MFM vs. control muscle. Notably the only gene with Ͼ1.0-fold log 2 upregulation comparing rest to after exercise in MFM horses was GADD45G, a gene that regulates DNA repair in response to stimuli such as oxidative stress (70). Following exercise, a remarkably focused lower expression of genes involved in cysteine transport and biosynthesis was identified in MFM compared with control muscle (Fig. 6). Cysteine beta synthase (CBS) and a cysteine transporter had Ͼ1.8 log 2 fold lower expression in postexercise MFM compared with postexercise control muscle. In addition, adenylhomocysteinase, a key enzyme involved in synthesizing cysteine from methionine, was one of only five genes/proteins common to the transcriptomic and proteomic data sets of MFM vs. control horses. Cysteine residues are among the most sensitive to oxidation relative to other amino acids, and alterations in the redox state of cysteine residues can modify protein structure and enzyme activity (14,26). Cysteine oxidation can result in a diverse array of modifications, some of which are reversible in vivo through interactions with glutathione, PRDX, or thioredoxin (14). Reversible cysteine oxidative modifications include intramolecular disulfides, S-nitrosylation, S-glutathionylation, and sulfonic acid, whereas irreversible cysteine oxidative modifications include sulfinic and sulfonic acid (14). One of the functions of PRXD6 and glutathione is to prevent irreversible oxidation of cysteine (Fig. 6). A deficiency of intracellular cysteine and PRDX6 in skeletal muscle following exercise in MFM horses could impair synthesis of cysteinedependent antioxidants such as PRDX6 and glutathione, which could further impair PRDX6 availability and lead to irreversible oxidation of cysteine, an essential protein in MFM muscle (Fig. 6). The possibility that MFM horses have an underlying defect in cysteine biosynthesis and the redox system is intriguing. Clearly, alterations in antioxidant function, cysteine uptake, and metabolism are worthy of further investigation in MFM horses.
Our initial focus in the present study was desmin and sarcomeric proteins because desmin aggregates and disruption of myofilaments and the Z-disc typify MFM in both equine and human cases (34,53). Although desmin aggregates were evident in histologic sections of equine MFM muscle, significant differences in desmin protein or gene expression were not detected (63). One synonymous and several intronic variants were identified in the desmin gene; however, there was not a mutation that was consistently present in all MFM horses and absent in control horses. Thus, a mutation in desmin does not appear to be responsible for Arabian MFM. Myotilin, filamin C, LIM protein, and ␣␤-crystallin protein content were similar between MFM and control horses. Furthermore, only one protein, a 26S proteosomal subunit, overlapped with proteins identified in a proteomic analysis of human patients with mutations in desmin and filamin C (40). The lack of overlap with sarcomere proteins in equine vs. human MFM could be because the present study used whole muscle samples vs. laser capture isolation of myofibers with desmin aggregates in the human study (40). Alternatively, the findings in the present study could suggest that the basis for equine MFM differs from desmin or filamin C-related MFM.
Desmin is a major target for oxidation and nitrosylation, which is well described in human MFM (31). In fact, the antioxidants N-acetyl cysteine and ␣-tocopherol have been shown to prevent desmin aggregation in cultured cells that express a desmin mutation (8,54). The role of oxidation in desmin-related myopathies is also highlighted by the fact that a myopathy once classified as MFM based on desmin aggregation, selenoprotein 1-related myopathy, was found to be caused by mutations in SEPN1, which is involved in oxidative homeostasis (3,9,21). Thus, there are several precedents for disturbances in oxidoreductase activity causing the type of desmin aggregation noted in horses with MFM.
Tropomyosin was the third significantly altered protein in MFM compared with control muscle. Lower tropomyosin (23.2 log 2 fold MFM) content in MFM muscle could be the result of protein oxidation or alternatively, could reflect a primary underlying disorder (30,35,36). Decreased tropomyosin content is known to occur in limb girdle muscular dystrophy 2A and hereditary inclusion body myositis (24), (55), which have different histopathologic appearances compared with MFM in horses (63,64). Mutations in the Z-disc protein MYOZ2 can cause myofibrillar disarray and hypertrophic cardiomyopathy and gene expression of MYOZ2 was significantly decreased in resting MFM vs. control horses (45). Muscle from horses with MFM also had significantly higher DE of a desmin cross-linking gene, desmoplakin, primarily expressed in cardiac muscle. Mutations in genes encoding other members of the plakin family, plectin and dystonin, disrupt the sarcomere in a fashion similar to that seen with MFM (37). Two elements involving the cytoskeleton (PDLIM1) and protein misfolding (HSPB7) were also common to the proteomic data set and the MFM vs. control DE transcripts after exercise. Deficiency of HSBP7 can cause myofibrillar disarray (32). Thus, tropomyosin, myozenin 2, desmoplakin, and HSPB7 are also candidates for further investigation in equine MFM.
The power of combining RNA-Seq and iTRAQ proteomic analysis in the presented study was highlighted by our ability to detect distinct but overlapping pathways such as DE of cysteine transporters and cysteine synthetic enzymes in RNA-Seq and cysteine-dependent PRDX6 content in proteomics. qRT-PCR was used to validate DE of those genes with Ͼ2-fold log 2 expression and revealed similar direction of changes to RNA-Seq data. With iTRAQ, we identified Ͼ600 confidently expressed proteins in equine muscle, which is more than the 400 proteins identified in human myopathies (55) and more than the 540 proteins found in studies of porcine muscle (29). Of the 673 proteins identified, 543 had expression in both the transcriptome and proteome. Little association was found, however, between DE of the same gene ID in the transcriptome and general expression in the proteome, similar to previous studies of myopathies such as murine muscular dystrophy (50). This may reflect the fact that long-term changes in gene expression are required to alter protein expression and that altered protein expression can arise not only from increased synthesis but also from enhanced or abnormal degradation (28).
Identification of many proteins and genes in the present study was limited by the lack of annotation provided in the current EquCab2.86 reference genome. Identification of small proteins was also potentially limited in the present study by use of whole muscle homogenates containing substantial amounts of large contractile proteins that could impair detection of less abundant smaller proteins. We used homogenates so that we were able to examine sarcomere and cytoskeletal proteins, which were of interest due to their importance in MFM. Previous iTRAQ analyses have utilized t-tests without correction for multiple testing to compare protein concentrations (48,71). We utilized stringent correction for multiple testing, resulting in identification of only three highly significant proteins. We also included those proteins with Ͼ2 log 2 fold changes to fully evaluate pathways associated with MFM.

Conclusions
In conclusion, this study defines for the first time a potential pathophysiological basis for an exertional myopathy in Arabian horses. Transcriptomic and proteomic profiles at rest and with exercise indicate that equine MFM is characterized by differential regulation of pathways involved in fatty acid metabolism, oxidoreductase activity, cysteine metabolism, sarcomere and cytoskeleton structure, as well as tissue regeneration and cell differentiation. A deficiency of cysteine-containing antioxidants and enhanced generation of ROS in MFM horses could lead to oxidative stress during exercise, oxidation of key proteins such as desmin and altered cysteine synthesis following endurance exercise.