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

,


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.

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 jGFOLD valuej ≥ 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.

Sequence read mapping and quantification of gene expression
A total of 20 RNA samples from light and normalcolored 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 log 2 (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).
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.

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).

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(King et al., , 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.

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.

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

Table 2 .
Significant KEGG pathways identified from the set of DEGs identified by DESeq2 and GFOLD a FDR-corrected P value.b In order of log(fold-change), highest to lowest.DEG = differentially expressed gene; KEGG = Kyoto Encyclopedia of Genes and Genomes.