Research Article

Transcriptomic Analysis for Pork Color—The Ham Halo Effect in Biceps Femoris

Authors: , , , , , ,

Abstract

Pork color is a major indicator of product quality that guides consumer purchasing decisions. Recently, industryhas received an increase in consumer complaints about the lightness and nonuniformity of ham color, primarily lighter colorin the periphery termed “ham halo” that is not caused by manufacturing procedures. This effect is seen in fresh and processed hams and the outer lighter muscle is associated with lower myoglobin concentration, pH, and type I fibers. The objective of this study was to identify differences in gene expression profiles between light and normal-colored portions of the biceps femoris muscle from pork hams. RNA-sequencing (RNA-seq) was performed for paired light and normal-colored muscle samples from 10 animals showing the ham halo effect. Over 50 million paired-end reads (2×75 bp) per library were obtained. An average of 99.74% of trimmed high-quality reads was mapped to the Sscrofa 11.1 genome assembly. Differentially expressed genes (DEGs) were identified using both the DESeq2 and GFOLD software packages. A total of 14,049 genes were expressed in the biceps femoris; 13,907 were expressed in both light and normal muscle, while 56 and 86 genes were only expressed in light and normal muscle, respectively. Analysis with DESeq2 identified 392 DEGs with 359 genes being more highly expressed in normal-colored muscle. A total of 61 DEGs were identified in the DESeq2 analysis and identified in at least 7 of the 10 individual animal analyses. All 61 of these DEGs were up-regulated in normal-colored muscle. Gene Ontology enrichment analysis of DEGs identified the transition between fast and slow fibers and skeletal muscle adaptation and contraction as the most significant biological process terms. The evaluation of gene expression by RNA-seq identified DEGs between regions of the biceps femoris with the ham halo effect that are associated with variation in pork color.

Keywords: pork color, ham halo, transcriptome

How to Cite: Nonneman, D. J. , Keel, B. N. , Lindholm-Perry, A. K. , Rohrer, G. A. , Wheeler, T. L. , Shackelford, S. D. & King, D. (2022) “Transcriptomic Analysis for Pork Color—The Ham Halo Effect in Biceps Femoris”, Meat and Muscle Biology. 6(1). doi: https://doi.org/10.22175/mmb.13050

Introduction

Uniform pink color is a primary determinant of consumer acceptance of cured ham products. Recently, the industry has received increased consumer feedback about the nonuniformity of ham color, primarily lighter color in the periphery termed “ham halo.” This effect, not caused by manufacturing procedures, is seen in fresh and processed hams. The outer, lighter muscle tissue (lower a* and b* and higher L* values) is associated with lower myoglobin concentration, lower pH, and lower percentage type I fibers (King et al., 2018). Color variation in hams has been correlated with L*, pH, and drip loss across muscles (McKeith and Pringle, 2013). These traits were predictive of the whole ham color, but other factors appear to be associated with color variation in individual muscles (Stufft et al., 2017). Other efforts to identify biological determinants of pork color have used proteomic contrasts of individual muscles varying in color and fiber type composition (Kim et al., 2018), transcriptomic profiling of different breeds (Liu et al., 2016), transcriptomic profiling in homologous muscles from individual animals differing in muscle fiber types (Ropka-Molik et al., 2017), or transcriptomic profiling across different muscles varying in muscle fiber type and color (Li et al., 2016; Zhu et al., 2016). These approaches have identified genes and pathways related to muscle fiber type, metabolism, and mitochondrial function (Li et al., 2016; Kim et al., 2018; Ropka-Molik et al., 2017; Liu et al., 2016). However, this study examined such variation across regions of the same muscle from the same animal.

Developments in high-throughput sequencing methods have provided new strategies for identifying genes and pathways associated with complex traits. The aim of this study was to compare gene expression profiles from light and normal-colored biceps femoris in order to identify potential biomarkers associated with the ham halo defect and meat color traits.

Materials and Methods

Tissue collection and RNA isolation

Postmortem ham muscles were selected from carcasses in the chilling cooler of a USDA-inspected processing facility; thus, animal care and use approval was not obtained for this experiment. Ham muscle samples were collected and frozen in liquid nitrogen at 2 h postmortem, as described by King et al. (2020), with the exception that samples for the current experiment were collected from the opposite carcass side than those collected for King et al. (2020). At 2 h postmortem, in the 10 hams sampled for sequencing, the outer lighter portion of the muscle (n = 10) had greater L* values (53.6 ± 0.93 vs. 48.3 ± 1.18; P < 0.001), lesser a* values (4.8 ± 0.62 vs. 9.8 ± 0.97; P = 0.005) and similar b* values compared to the inner normal-colored portion. Samples were stored at −80° C. Approximately 100 mg each of light and normal-colored biceps femoris muscle was taken 2.54 cm apart from each of the 10 hams for RNA isolation. RNA was extracted using TRIzol (Thermo Fisher Scientific, Lenexa, KS) and chloroform. Genomic DNA was removed from the total RNA with the RNeasy Mini Kit and RNase-Free DNase Kit (Qiagen, Germantown, MD) according to the manufacturer’s protocol. RNA concentration was determined with a Nanodrop ND-1000 Spectrophotometer (Wilmington, DE). The average 260/280 ratio was 2.01, with a range of 1.94 to 2.07. An Agilent 2100 Bioanalyzer (Santa Clara, CA) was used to determine the RNA integrity number (RIN). RIN ranged from 8.0 to 8.9, with an average of 8.52. The libraries were prepared for RNA-sequencing (RNA-seq) with the Illumina TruSeq Stranded mRNA High Sample kit following the manufacturer’s protocol (Illumina Inc., San Diego, CA). Libraries were quantified with qRT-PCR using the NEBNext Library Quant Kit for Illumina (NEB, Ipswich, MA) kit on a Bio-Rad CFX384 Real Time System (Hercules, CA) thermal cycler, and the quality of each library was determined with an Agilent Bioanalyzer DNA kit. All 20 libraries were paired-end sequenced with 150 cycle high output sequencing kits using the Illumina NextSeq 500 platform.

Processing RNA-sequencing data

Alignment of RNA-seq reads was carried out as follows. First, the quality of the raw paired-end sequence reads in individual fastq files was assessed using FastQC (Version 0.11.5; www.bioinformatics.babraham.ac.uk/projects/fastqc), and reads were trimmed to remove adapter sequences and low-quality bases using the Trimmomatic software (Version 0.35; Bolger et al., 2014). The remaining reads were mapped to the Sscrofa 11.1 genome assembly using Hisat2 (Version 2.1.0; Kim et al., 2015) with its default parameters. The StringTie software (Pertea et al., 2015) was used to calculate raw read counts for each of the 29,651 annotated genes in the NCBI Sscrofa 11.1 reference annotation (Release 106). Genes with low read counts, <15 reads in at least 10 samples, were removed resulting in a set of 14,049 genes to be used in downstream analyses.

Identification of differentially expressed genes

Differential gene expression was tested using DESeq2 (Version 1.26; Love et al., 2014) and GFOLD (Version 1.1.4; Feng et al., 2012), which are designed for use with data that have multiple biological replicates and no biological replicates, respectively. DESeq2 was used to compare expression differences between the light-colored tissues and normal-colored tissues. A gene was considered as differentially expressed at FDR-adjusted P ≤ 0.05.

The GFOLD software package was used to identify differentially expressed genes (DEGs) between the light and normal tissue samples from each individual animal. GFOLD assigns reliable statistics for expression changes based on the posterior distribution of log fold change. When there are no biological replicates, the one assumption underlying the GFOLD algorithm is that the read count of a gene follows a Poisson distribution. DEGs were identified with GFOLD configured for a 99% confidence interval. The generalized fold change (GFOLD value), which can be considered a reliable log2 fold change, was used to assess whether a given gene was differentially expressed (Feng et al., 2012). GFOLD was run using the DESeq2 normalization scheme. A gene was considered differentially expressed when |GFOLDvalue|0.5.

Gene Ontology enrichment analysis

Enrichment analysis of gene function was performed using the Advaita Bio iPathwayGuide software (Advaita Bio, http://advaitabio.com/ipathwayguide, accessed 26 March 2019). For each Gene Ontology (GO) term, the number of DEGs annotated to the term is compared against the number of genes expected by chance. iPathwayGuide uses an overrepresentation test, based on a hypergeometric distribution, to compute the statistical significance of observing more than the expected number of DEGs. Significance of GO terms was assessed using the default Mus musculus GO annotation as background. A GO term was considered statistically significant at FDR-corrected P ≤ 0.05.

Pathway analysis

The Advaita Bio iPathwayGuide tool was used to analyze potential Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways impacted by our set of DEGs. Briefly, differential expression results were uploaded to the iPathwayGuide tool, and gene identifiers were mapped onto the Mus musculus gene set. Pathway overrepresentation analysis was performed by comparing the number of DEGs associated with a pathway between light and normal tissue groups. The P value for each pathway was considered statistically significant at FDR-corrected P ≤ 0.05.

Identification of differentially expressed genes in mitochondrial complexes and pork quality quantitative trait locus regions

Differentially expressed mitochondrial complex genes were identified from the HGNC Database (https://www.genenames.org/data/genegroup/#!/group/639, accessed 30 August 2021). DEGs residing in single nucleotide polymorphism (SNP) association quantitative trait locus (QTL) were identified as follows. Swine QTL from the Sscrofa 11.1 genome build were downloaded from the Animal QTL database (Hu et al., 2005; http://www.animalgenome.org/cgi-bin/QTLdb/SS/index, accessed 27 April 2021), and SNP associations for pork quality within 100 kb of the genes were considered potentially associated with the DEGs. Only trait associations related to color (a*, b*, L* and meat color score [mcolor]), fiber type (FIB2BP), drip loss (DRIPL), cooking loss (COOKL), water-holding capacity (WHC), pH, and conductivity (COND) were considered.

Results

Sequence read mapping and quantification of gene expression

A total of 20 RNA samples from light and normal-colored biceps femoris tissue from 10 animals were sequenced on a NextSeq 500 instrument, resulting in over one billion raw paired-end reads. The total number of sequence reads generated for each sample ranged from 46,081,506 to 90,404,370, with a median of 54,871,505 [Table S1(A)]. Adapter sequences and low-quality bases were trimmed, and the resulting high-quality reads were mapped to the Sscrofa11.1 genome assembly with an average 99.7% overall mapping rate per library. The sequence data are deposited in the NCBI SRA database under project number PRJNA753065.

Among the 29,651 genes in the NCBI reference annotation, 14,049 were expressed in the biceps femoris transcriptome [Table S1(B)]. A total of 13,907 genes were expressed in both the light and normal-colored tissue, while 56 genes were uniquely expressed in the light tissue [mean read counts >15 in light tissue only; Table S1(C)], and 86 were uniquely expressed in the normal tissue [mean read counts >15 in normal tissue only; Table S1(D)]. Because the myoglobin (MB) transcript is not annotated in the current build, we were unable to detect expression of MB.

Differentially expressed genes between light and normal-colored biceps femoris

To identify potential candidate genes related to the halo condition within the biceps femoris muscle, we identified DEGs between light (outer) and normal (inner) biceps femoris tissue from 10 animals using 2 different methods, DESeq2 and GFOLD. DESeq2 was used to analyze expression differences between light and normal tissues across all animals. Of the 14,049 expressed genes in the biceps femoris tissue, this analysis identified a total of 392 DEGs between the normal and light groups (Table S2). The sign of the log2(fold-change) was used to classify the DEGs as up- or down-regulated in the normal. A total of 359 DEGs were found to be up-regulated in normal, while the remaining 33 were down-regulated (Figure 1 and Table S2). Only 6 uniquely expressed genes (only expressed in normal or light), LMX1B, GPAT2, LOC100524035, LOC102161099, LOC110261129, and PAQR9, were significantly differentially expressed (Table S2).

Figure 1.
Figure 1.

(A) Mean normalized expression and Log2 fold change of transcripts identified in Biceps femoris. Significantly differentially expressed genes (DEGs) are shown in red. (B) Volcano plot of significantly differentially expressed genes in normal (BfN) vs. light (BfL) colored Biceps femoris by DESeq2 and GFOLD. BfN<BfL DEGs are shown in red, BfN>BfL DEGs are shown in blue and green, and GFOLD DEGs are shown in green.

To eliminate potential differences due to individual animal variation, GFOLD was used to identify DEGs between the light and normal-colored tissue within each animal. The number of DEGs identified in a single animal ranged from 302 to 1,418, with an average of 694.5 genes. There were 61 DEGs that were identified by DESeq2 in at least 7 of the 10 individual animal analyses (Table S2). In this set of genes, all were up-regulated in normal-colored biceps femoris.

Gene Ontology enrichment analysis

GO enrichment analysis of the 392 DEGs identified by DESeq2 identified 231 biological process (BP) terms, 105 molecular function (MF) terms, and 86 cellular component (CC) terms as significantly enriched (Table S3). The most abundant terms were metabolic process (BP; 248 DEGs), catalytic activity (MF; 184 DEGs), and intracellular organelles (CC; 321 DEGs). The 5 most significant BP terms were oxidation-reduction process (P = 1.00E-24), small molecule process (P = 1.00E-24), ATP metabolic process (P = 1.00E-24), cellular respiration (P = 1.00E-24), and purine nucleoside monophosphate metabolic process (P = 1.00E-24). MF terms oxidoreductase activity (P = 1.00E-24), cofactor binding (P = 5.70E-18), and hydrogen ion transmembrane transporter activity (P = 4.91E-14) were the most significant. Nineteen CC terms shared the most significant P value (P = 1.00E-24). These terms are predominantly related to mitochondrial function and are reflected by an overrepresentation of more highly expressed mitochondrial complex genes in all mitochondrial complexes I, II, III, IV, and V in normal-colored biceps femoris (Table 1, Table S2).

Table 1.

Differentially expressed mitochondrial complex genes by DESeq2

Complex Complex Name # Genes # DEGs Function
I NADH:ubiquinone oxidoreductase subunits 44 26 Electron transport
II Succinate dehydrogenase subunits 4 3 TCA cycle and electron transport
III CoQ:Ubiquinol cytochome C reductase complex subunits 10 7 Oxidative-phos and electron transport
IV Cytochrome C oxidase subunits 19 11 Electron transport
V ATP synthase subunits 20 16 ATP synthesis
  • DEG = differentially expressed gene.

GO enrichment analysis of the set of 61 DEGs that overlapped between the DESeq2 and the individual animal GFOLD analyses identified 194 BP terms, 5 MF terms, and 20 CC terms (Table S4). The most abundant terms were multicellular organismal process (BP; 28 DEGs), cytoskeletal protein binding (MF; 13 DEGs), and cytoskeleton (CC; 15 DEGs). The most significant BP terms were muscle contraction (P = 1.02E-07), striated muscle contraction (P = 1.02E-07), transition between fast and slow fiber (P = 1.02E-07), regulation of skeletal muscle adaptation (P = 1.17E-06), and muscle system process (P = 1.34E-06), while MF terms cytoskeletal protein binding (P = 0.003), actin binding (P = 0.003), and troponin T binding (P = 0.013) were most significant. The top CC terms included myofibril (P = 8.57E-07), contractile fiber (P = 8.57E-07), actin cytoskeleton (P = 8.57E-07), and myosin complex (P = 5.59E-06).

Pathway analysis

Pathway analysis was performed using the iPathwayGuide tool in order to characterize the potential pathways impacted by gene expression differences between light and normal-colored biceps femoris tissue. Fifteen pathways were identified as significantly enriched for the set of 61 DEGs (Table 2). The most significant pathways included cardiac muscle contraction (P = 8.79E-08), hypertrophic cardiomyopathy (P = 2.35E-05), adrenergic signaling in cardiomyocytes (P = 2.60E-05), and dilated cardiomyopathy (P = 4.10E-05).

Table 2.

Significant KEGG pathways identified from the set of DEGs identified by DESeq2 and GFOLD

Pathway P Valuea,b # DEGs in Pathway Genes in Pathways
Cardiac muscle contraction(KEGG: 04260) 8.79E-08 7 MYL2, TNNC1, MYH7, ATP2A2, MYL3, COX1
Hypertrophic cardiomyopathy(KEGG: 05410) 2.35E-05 6 MYL2, TNNC1, MYH7, TPM3, ATP2A2, MYL3
Adrenergic signaling in cardiomyocytes (KEGG: 04261) 2.60E-05 7 MYL2, PLN, TNNC1, MYH7, TPM3, ATP2A2, MYL3
Dilated cardiomyopathy (KEGG: 05414) 4.10E-05 7 MYL2, PLN, TNNC1, MYH7, TPM3, ATP2A2, MYL3
Parkinson’s disease (KEGG: 05012) 0.026 4 PPIF, ATP8, COX1, CYCS
Thyroid hormone signaling pathway (KEGG: 04919) 0.026 4 PLN, MED12L, ATP2A2, RCAN2
Alzheimer’s disease (KEGG: 05010) 0.029 5 ATP2A2, LPL, ATP8, COX1, CYCS
cGMP-PKG signaling pathway (KEGG: 04022) 0.029 4 PLN, MYH7, ATP2A2, PPIF
PPAR signaling pathway (KEGG: 03320) 0.029 3 PLIN5, FABP3, LPL
Huntington’s disease (KEGG: 05016) 0.030 4 PPIF, ATP8, COX1, CYCS
Tight junction(KEGG: 04530) 0.030 3 MYL6B, MYH7B, MYL2
Citrate cycle (TCA cycle) (KEGG: 00020) 0.037 2 MDH1, PDHA1
HIF-1 signaling pathway (KEGG: 04066) 0.037 3 TFRC, VEGFA, PDHA1
Pyruvate metabolism(KEGG: 00620) 0.047 2 MDH1, PDHA1
Carbon metabolism(KEGG: 01200) 0.050 3 MDH1, PDHA1
  • FDR-corrected P value.

  • In order of log(fold-change), highest to lowest.

  • DEG = differentially expressed gene; KEGG = Kyoto Encyclopedia of Genes and Genomes.

Identification of DEGs in pork quality quantitative trait locus regions

A total of 66 unique SNP associations with pork color (a*, b*, L*, and mcolor), COND, fiber type, pH, and DRIPL (also COOKL and WHC) were detected within 100 kb of 47 DEGs (Figure 2, Table S2). The most prevalent trait associations associated with DEGs were color (25), pH (18), and DRIPL (9). Most DEG alignments with SNP associations were found on chromosomes 1 (4), 3 (5), 6 (16), and 12 (5) and alignments reflect the number of QTL associations for each trait and number of QTL for each chromosome (https://www.animalgenome.org/cgi-bin/QTLdb/SS/summary).

Figure 2.
Figure 2.

Manhattan plot of pork quality SNP associations residing in differentially expressed genes.

Discussion

Consumer complaints regarding a lack of uniformity of cured ham color led to the identification of the ham halo condition. The condition was observed in raw materials from many sources and was determined to be present in live pigs (King et al., 2018). Our previous efforts to characterize the ham halo condition indicate that the light tissue on the periphery of the muscle is much lighter (lower L*) and less red (lower a*) with a concomitant reduction in myoglobin concentration relative to the unaffected deep portion of the muscle (King et al., 2018). The changes in color and pigment also were associated with dramatically increased proportions of white muscle fibers. This condition has been identified in all of the ham muscles examined regardless of genetics or stage of production. However, the severity of the color defect differs across sire groups (King et al., 2020). Thus, genetic variation exists in the halo condition and could be used to mitigate the condition.

A comparative transcriptomic analysis of light and normal-colored biceps femoris revealed genes and pathways related to muscle fiber type and possibly color that is supportive of previous reports (King et al., 2018). The normal-colored portion of the muscle had more DEGs more highly expressed than the light portion. These results are comparable with what has been previously reported for differences across different muscles of contrasting fiber type and color (Zhu et al., 2016) and in homologous muscles differing in fiber type composition (Ropka-Molik et al., 2017). Because the initial characterization of this defect found differences in myoglobin concentrations, color scores, pH, and fiber type in severely affected hams (King et al., 2018), we focused on relating these traits with DEGs.

Several DEGs in significant KEGG pathways were more highly expressed in the normal-colored portion such as TTNC1, TNNI1, TNNI3, and TNNT1, which are all associated with slow red skeletal muscle fibers (Sheng and Jin, 2016; Wei and Jin, 2016). Also, myosin heavy chains MYBPC1, MYH7, MYH7B, and MYOZ2 and light chains MYL2, MYL3, and MYL6B were more highly expressed in the normal-colored portion and are more highly expressed in cardiac and red slow twitch muscle (Stuart et al., 2016; Schiaffino, 2018). Actin 3 (ACTN3) and tropomyosins TPM1 and TPM3 were differentially expressed, with ACTN3 and TPM1 being more lowly expressed and TPM3 more highly expressed in the normal-colored portion, as would be expected with the differences in fiber type (Oe et al., 2016). The expression of these canonical genes that are related to muscle fiber type is under epigenetic regulation presumably by promoter hypermethylation (Oe et al., 2021). Also more highly expressed in the normal-colored portion, HADHB and SDHB proteins were differentially expressed in weanling pigs fed different leucine diets associated with the developmental transition of slow to fast muscle fibers (Fan et al., 2017). These differences in expression could lead to insights for the development of more accurate assessments of muscle fiber type composition and metabolic state.

Differential expression of mitochondrial complex genes is also indicative of changes in muscle fiber type (Liu et al., 2016; Kim et al., 2018; Lang et al., 2018). Differential gene expression in the 2 portions of the muscle indicate a substantial shift in oxidative capacity of the muscle based on color change and fiber type. While gene expression measurements of mitochondrial complex genes were greater in the normal-colored portion, it is not known if these changes are due to increased expression or increased mitochondrial number (Picard et al. 2012; Liu et al 2017; Lee et al., 2018).

Some of the DEGs that coaligned with SNP associations for pork color, pH, drip loss, or related traits have been proposed as candidate genes for pork quality. SNPs in DHRS4 and MDH1 (that reside within pH associations) have been associated with 45-min and ultimate pH (Choi et al., 2012; Hwang et al., 2017). Haplotypes in the COQ9 gene have been associated with color and drip loss (Brunner et al., 2012), and polymorphisms in the TNNI1 gene have been associated with pH and pork color (Yang et al., 2010). While gene expression is associated with phenotypic differences, these genetic associations suggest that variation within these genes may affect these traits as well.

The halo condition is a color defect within individual ham muscles. Numerous genes that were differentially expressed between the halo and normal portions of the biceps femoris muscle lie in close proximity to quantitative trait loci previously identified for pork quality traits. Most of the DEGs from the present experiment were more highly expressed in the normal tissue (i.e., more red muscle fibers and consequently, darker, more red color, and greater myoglobin concentration). Thus, greater expression of these genes should improve pork quality. The observation of QTL for pork quality traits corresponding to DEGs in the present experiment may uncover novel candidate genes and identification of functional variants that could potentially be used for genetic selection related to pork quality.

These results support previous work to characterize this condition (King et al., 2018, 2020) and indicate that RNA-seq provides a global, comprehensive view of the transcriptome of phenotypic differences within the same muscle, candidate gene identification for related traits, and insights into the metabolic pathways responsible for the phenotypic differences in lean color of the ham halo defect.

Acknowledgements

The authors acknowledge Brett Brockman and Linda Flathman for library preparation, the US Meat Animal Research Center DNA Core Facility for data collection, and Stephanie Schmidt for manuscript preparation.

Mention of trade names or commercial products in this publication is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the USDA.

The USDA prohibits discrimination in all its programs and activities on the basis of race, color, national origin, age, disability, and where applicable, sex, marital status, familial status, parental status, religion, sexual orientation, genetic information, political beliefs, reprisal, or because all or part of an individual’s income is derived from any public assistance program (not all prohibited bases apply to all programs.). Persons with disabilities who require alternative means for communication of program information (Braille, large print, audiotape, etc.) should contact USDA’s TARGET Center at (202) 720-2600 (voice and TDD). To file a complaint of discrimination, write to USDA, Director, Office of Civil Rights, 1400 Independence Avenue, SW, Washington, DC 20250-9410, or call (800) 795-3272 (voice) or (202) 720-6382 (TDD). USDA is an equal opportunity provider and employer.

Literature Cited

Advaita Bio. iPathwayGuide. http://advaitabio.com/ipathwayguide. (Accessed 26 March 2019.)

Bolger, A. M., M. Lohse, and B. Usadel. 2014. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 30:2114–2120. doi: https://doi.org/10.1093/bioinformatics/btu170.

Brunner, R. M., T. Srikanchai, E. Murani, K. Wimmers, and S. Ponsuksili. 2012. Genes with expression levels correlating to drip loss prove association of their polymorphism with water holding capacity of pork. Mol. Biol. Rep. 39:97–107. doi: https://doi.org/10.1007/s11033-011-0714-5.

Choi, I., R. O. Bates, N. E. Raney, J. P. Steibel, and C.W. Ernst. 2012. Evaluation of QTL for carcass merit and meat quality traits in a US commercial Duroc population. Meat Sci. 92:132–138. doi: https://doi.org/10.1016/j.meatsci.2012.04.023.

Fan, Q., B. Long, G. Yan, Z. Wang, M. Shi, X. Bao, J. Hu, X. Li, C. Chen, Z. Zheng, and X. Yan. 2017. Dietary leucine supplementation alters energy metabolism and induces slow-to-fast transitions in longissimus dorsi muscle of weanling piglets. Brit. J. Nutr. 117:1222–1234. doi: https://doi.org/10.1017/S0007114517001209.

Feng, J., C. A. Meyer, Q. Wang, J. S. Liu, X. S. Liu, and Y. Zhang. 2012. GFOLD: a generalized fold change for ranking differentially expressed genes from RNA-seq data. Bioinformatics. 28:2782–2788. doi: https://doi.org/10.1093/bioinformatics/bts515.

HGNC Database. HUGO Gene Nomenclature Committee (HGNC). www.genenames.org. (Accessed 30 August 2021.)

Hu, Z.-L., S. Dracheva, W. Jang, D. Maglott, J. Bastiaansen, M. F. Rothschild, and J. M. Reecy. 2005. A QTL resource and comparison tool for pigs: PigQTLDB. Mamm. Genome. 16:792–800. doi: https://doi.org/10.1007/s00335-005-0060-9.

Hwang, J. H., S. M. An, S. G. Kwon, D. H. Park, T. W. Kim, D. G. Kang, G. E. Yu, I.-S. Kim, H. C. Park, J. Ha, and C. W. Kim. 2017. Associations of the polymorphisms in DHRS4, SERPING1, and APOR genes with postmortem pH in Berkshire pigs. Anim. Biotechnol. 28:288–293. doi: https://doi.org/10.1080/10495398.2017.1279171.

Kim, D., B. Langmead, and S. L. Salzberg. 2015. HISAT: a fast spliced aligner with low memory requirements. Nat. Methods. 12:357–360. doi: https://doi.org/10.1038/nmeth.3317.

Kim, G.-D., H.-S. Yang, and J.-Y. Jeong. 2018. Intramuscular variations of proteome and muscle fiber type distribution in semimembranosus and semitendinosus muscles associated with pork quality. Food Chem. 244:143–152. doi: https://doi.org/10.1016/j.foodchem.2017.10.046.

King, D. A., S. D. Shackelford, D. Nonneman, G. A. Rohrer, and T. L. Wheeler. 2020. Sire variation in the severity of the ham halo condition. Meat Muscle Biol. 4:1–15. doi: https://doi.org/10.22175/mmb.9743.

King, D. A., S. D. Shackelford, T. Schnell, L. Pierce, and T. L. Wheeler. 2018. Characterizing the ham halo condition: A color defect in fresh pork biceps femoris muscle. Meat Muscle Biol. 2:205–213. doi: https://doi.org/10.22175/mmb2018.02.0001.

Lang, F., S. Khaghani, C. Türk, J. L. Wiederstein, S. Hölper, T. Piller, L. Nogara, B. Blaauw, S. Günther, S. Müller, T. Braun, and M. Krüger. 2018. Single muscle fiber proteomics reveals distinct protein changes in slow and fast fibers during muscle atrophy. J. Proteome Res. 17:3333–3347. doi: https://doi.org/10.1021/acs.jproteome.8b00093.

Lee, H., K. Kim, B. Kim, J. Shin, S. Rajan, J. Wu, X. Chen, M. D. Brown, S. Lee, and J.-Y. Park. 2018. A cellular mechanism of muscle memory facilitates mitochondrial remodelling following resistance training. Journal of Physiology. 596:4413–4426. doi: https://doi.org/10.1113/JP275308.

Li, B., C. Dong, P. Li, Z. Ren, H. Wang, F. Yu, C. Ning, K. Liu, W. Wei, R. Huang, J. Chen, W. Wu, and H. Liu. 2016. Identification of candidate genes associated with porcine meat color traits by genome-wide transcriptome analysis. Sci. Rep.-UK. 6:35224. doi: https://doi.org/10.1038/srep35224.

Liu, X., N. Trakooljul, E. Muráni, C. Krischek, K. Schellander, M. Wicke, K. Wimmers, and S. Ponsuksili. 2016. Molecular changes in mitochondrial respiratory activity and metabolic enzyme activity in muscle of four pig breeds with distinct metabolic types. J. Bioenerg. Biomembr. 48:55–65. doi: https://doi.org/10.1007/s10863-015-9639-3.

Liu, X., N. Trakooljul, F. Hadlich, E. Murani, K. Wimmers, and S. Ponsuksili. 2017. Mitochondrial-nuclear crosstalk, haplotype and copy number variation distinct in muscle fiber type, mitochondrial respiratory and metabolic enzyme activities. Sci. Rep.-UK. 7:14024. doi: https://doi.org/10.1038/s41598-017-14491-w.

Love, M. I., W. Huber, and S. Anders. 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15:550. doi: https://doi.org/10.1186/s13059-014-0550-8.

McKeith, R. O., and T. D. Pringle. 2013. Quality attributes and color characteristics in three-piece boneless hams. Meat Sci. 95:59–63. doi: https://doi.org/10.1016/j.meatsci.2013.04.020.

Oe, M., K. Ojima, and S. Muroya. 2021. Difference in potential DNA methylation impact on gene expression between fast- and slow-type myofibers. Physiol. Genomics. 53:69–83. doi: https://doi.org/10.1152/physiolgenomics.00099.2020.

Oe, M., K. Ojima, I. Nakajima, K. Chikuni, M. Shibata, and S. Muroya. 2016. Distribution of tropomyosin isoforms in different types of single fibers isolated from bovine skeletal muscles. Meat Sci. 118:129–132. doi: https://doi.org/10.1016/j.meatsci.2016.04.013.

Pertea, M., G. M. Pertea, C. M. Antonescu, T.-C. Chang, J. T. Mendell, and S. L. Salzberg. 2015. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33:290–295. doi: https://doi.org/10.1038/nbt.3122.

Picard, M., R. T. Hepple, and Y. Burelle. 2012. Mitochondrial functional specialization in glycolytic and oxidative muscle fibers: tailoring the organelle for optimal function. Am J Physiol Cell Physiol. 302:C629–C641. doi: https://doi.org/10.1152/ajpcell.00368.2011.

Ropka-Molik, K., A. Bereta, K. Żukowski, K. Piórkowska, A. Gurgul, and G. Żak. 2017. Transcriptomic gene profiling of porcine muscle tissue depending on histological properties. Anim. Sci. J. 88:1178–1188. doi: https://doi.org/10.1111/asj.12751.

Schiaffino, S. 2018. Muscle fiber type diversity revealed by anti-myosin heavy chain antibodies. FEBS J. 285:3688–3694. doi: https://doi.org/10.1111/febs.14502.

Sheng, J.-J., and J.-P. Jin. 2016. TNNI1, TNNI2 and TNNI3: Evolution, regulation, and protein structure-function relationships. Gene. 576:385–394. doi: https://doi.org/10.1016/j.gene.2015.10.052.

Stuart, C. A., W. L. Stone, M. E. Howell, M. F. Brannon, H. K. Hall, A. L. Gibson, and M.H. Stone. 2016. Myosin content of individual human muscle fibers isolated by laser capture microdissection. Am. J. Physiol.-Cell Ph. 310:C381–C389. doi: https://doi.org/10.1152/ajpcell.00317.2015.

Stufft, K., J. Elgin, B. Patterson, S. K. Matarneh, R. Preisser, H. Shi, E. M. England, T. L. Scheffler, E. W. Mills, and D. E. Gerrard. 2017. Muscle characteristics only partially explain color variations in fresh hams. Meat Sci. 128:88–96. doi: https://doi.org/10.1016/j.meatsci.2016.12.012.

Wei, B., and J.-P. Jin. 2016. TNNT1, TNNT2, and TNNT3: Isoform genes, regulation, and structure–function relationships. Gene. 582:1–13. doi: https://doi.org/10.1016/j.gene.2016.01.006.

Yang, H., Z. Y. Xu, M. G. Lei, F. E. Li, C. Y. Deng, Y. Z. Xiong, and B. Zuo. 2010. Association of 3 polymorphisms in porcine troponin I genes (TNNI1 and TNNI2) with meat quality traits. J. Appl. Genet. 51:51–57. doi: https://doi.org/10.1007/BF03195710.

Zhu, J., X. Shi, H. Lu, B. Xia, Y. Li, X. Li, Q. Zhang, and G. Yang. 2016. RNA-seq transcriptome analysis of extensor digitorum longus and soleus muscles in large white pigs. Mol. Genet. Genomics. 291:687–701. doi: https://doi.org/10.1007/s00438-015-1138-z.

Supplemental Material

The sequence data is deposited in the NCBI SRA database under accession number PRJNA753065. Supplemental material for this article can be found online at https://doi.org/10.22175/mmb.13050.

Table S1. Sequence and mapping summary.

Table S2. DESeq2 and GFOLD DEGs.

Table S3. DESeq2 DEG GO Terms.

Table S4. GFOLD DEG GO Terms.