lation (C), oxidation (M), and deamidation (NQ). Protein identification could be verified when at least two peptides derived from one protein showed statistically significant identity. The Mascot score was 33 with a False Discovery Rate of 1%, and the *p*-value was adjusted at a given threshold (0.0093).

## 2.5. *Bioinformatics Analyses*

### 2.5.1. Protein-Protein Interactions (PPI)

The protein-protein interactions between the putative protein biomarkers were analysed using the STRING web service database (https://string-db.org/, accessed on 28 November 2020). Default settings were used, i.e., medium confidence of 0.4 and 4 criteria for linkage: co-occurrence, experimental evidence, existing databases, and text mining. As the bovine Gene Ontology (GO) had limits, orthologous human Uniprot IDs, following the procedure by Gagaoua et al. [14], were used for this analysis to take advantage of the most complete annotations available.

### 2.5.2. Gene Ontology and Pathway and Process Enrichment Analyses

The pathway and Gene Ontology analyses were performed using two web-based tools, First, ProteINSIDE (http://www.proteinsinside.org/, accessed on 28 November 2020) was used to investigate GO terms for potential functions and molecular mechanisms [23]. For this analysis, the top 20 GO enrichment terms (*p*-value, Benjamini–Hochberg < 0.05) were considered and covered Biological Process (BP), Molecular Function (MF), and Cellular Component (CC) categories. The Metascape® (https://metascape.org/, accessed on 28 November 2020) web service tool was further used to investigate the pathway and process enrichment analyses using the list of 34 differential proteins. The statistically significant enriched ontology terms were displayed based on the hypergeometric test and Benjamini–Hochberg *p*-value correction algorithm [24].

## 2.6. *Statistical Analyses*

Statistical analyses of protein abundance were performed with XLSTAT 2018.2 (AddinSoft, Paris, France), as well as the online tools NormalyzerDE and MetaOmGraph, mainly for data standardisation. Raw data were scrutinised for data entry errors, any missing data, or outliers. Log2 transformation and mean normalisation were performed on protein abundance among replicate samples. For the comparison of protein abundance between the tender (low WBSF values) and tough meat samples (high WBSF values), a one-way analysis of variance was performed for each protein. Differences in protein abundance between the tender and tough groups were considered significant at *p* < 0.05, and significant proteins were considered as candidate protein biomarkers. Pearson correlations were computed between the individual WBSF values and protein abundances for those proteins significant following ANOVA. Correlations were considered significant at *p* < 0.05. To get an overview of the main proteins related to WBSF variability, Partial Least Squares (PLS) regressions on standardised data were conducted to generate explanatory models using the list of the candidate protein biomarkers and identify the most influential proteins based on the variable importance in projection (VIP) filter set at both VIP > 1.0 and >0.8, as described by Gagaoua et al. [9]. Moreover, a stepwise regression analysis was used to explain WBSF using the 34 differential proteins (as independent variables, x). The absence of collinearity was systematically tested [25], specifically, the variable was identified as collinear if it possessed a high condition index > 10. The regression model allowed the entry of no more than 3 explanatory variables based on the parsimony principle.

# 3. Results

## 3.1. *Differential Proteins between Extreme Groups of High and Low WBSF Values*

According to Huffman et al. [26], for beef cooked to 70 °C, meat with Warner–Bratzler shear force values of 4.1 kg (40.18 N) or less was correlated with high levels (98%) of consumer acceptability, while beef prepared under the same conditions with shear force
values of 5.8 kg (56.84 N) or greater remained unacceptable. Two groups of beef samples with a large difference in shear force were selected from a panel of 107 beef animals collected and profiled under similar conditions. The mean shear force value in the lower shear force group was 33.21 N, while the mean shear force value for the other group was 63.96 N. These groups were classified as tender and tough, respectively. Putative protein biomarkers of beef tenderness that significantly differed in abundance in muscle samples were identified from these divergent groups (Table 1 and details in Table S1).

<table><caption>**Table 1.** Warner–Bratzler shear force (WBSF) values of the *Longissimus thoracis* muscles used in this trial.</caption><thead><tr><th>Quality Traits</th><th>Min</th><th>Max</th><th>Mean</th><th>SD</th><th>CV (%)</th></tr></thead><tbody><tr><td>WBSF (N) (n = 9)</td><td>27.70</td><td>38.85</td><td>33.21</td><td>3.24</td><td>9.75</td></tr><tr><td>WBSF (N) (n = 9)</td><td>59.25</td><td>71.40</td><td>63.96</td><td>3.98</td><td>6.22</td></tr></tbody></table>

N, Newtons; SD, Standard Deviation; CV, Coefficient of Variation.

A total of 34 proteins were different (*p* < 0.05) in their abundance between the tender and tough groups (Table 2). These 34 proteins belonged to five major biological pathways (Table 2), these being: (i) muscle contraction, structure, and associated proteins (n = 17; 50%); (ii) energy metabolism and associated pathways (n = 5; 15%); (iii) heat shock proteins (n = 4; 12%); (iv) oxidative stress (n = 2; 6%); and (v) other pathways including regulation of cellular processes, binding, apoptotic, and transport proteins (n = 6; 17%). The 34 proteins were then compared with a database of beef tenderness biomarkers by Gagaoua et al. [13], of which 23 overlapped with the database (Table 2).

<table><caption>**Table 2.** List of the 34 differential proteins organised by biological family, identified to significantly differ among the two WBSF (tenderness) groups.</caption><thead><tr><th rowspan="2">Uniprot ID</th><th rowspan="2">Gene Name</th><th rowspan="2">Full Protein Name</th><th colspan="2">Differences</th><th>Pearson Correlations <sup>a</sup></th><th rowspan="2">Overlap with Gagaoua et al. Database [13]</th></tr><tr><th>Fold Change (Log2)</th><th>p-Value</th><th>WBSF</th></tr></thead><tbody><tr><td colspan="7" style="text-align: left; font-weight: bold;">Muscle contraction, structure and associated proteins (n = 17)</td></tr><tr><td>Q08DI7</td><td>MYOZ3 <sup>b</sup></td><td>Myozenin 3</td><td>-0.53</td><td>0.002</td><td>0.741 ***</td><td>✓</td></tr><tr><td>E1BNG8</td><td>BIN1</td><td>Bridging Integrator-1</td><td>-0.25</td><td>0.003</td><td>0.616 **</td><td></td></tr><tr><td>Q148F1</td><td>CFL2</td><td>Cofilin-2</td><td>-0.47</td><td>0.003</td><td>0.670 **</td><td></td></tr><tr><td>E1BIN0</td><td>FHOD1</td><td>Formin homology 2 domain containing 1</td><td>-0.40</td><td>0.005</td><td>0.714 ***</td><td></td></tr><tr><td>A6QLZ8</td><td>CORO6</td><td>Coronin</td><td>-0.57</td><td>0.007</td><td>0.607 **</td><td></td></tr><tr><td>Q0P571</td><td>MYLPF</td><td>Myosin regulatory light chain 2</td><td>-0.71</td><td>0.009</td><td>0.613 **</td><td>✓</td></tr><tr><td>Q3SX40</td><td>PDLIM7 <sup>b</sup></td><td>PDZ and LIM domain protein 7</td><td>-0.47</td><td>0.010</td><td>0.642 **</td><td>✓</td></tr><tr><td>Q0VC48</td><td>TMOD4</td><td>Tropomodulin-4</td><td>-0.42</td><td>0.010</td><td>0.640 **</td><td>TMOD1</td></tr><tr><td>Q2KJH4</td><td>WDR1</td><td>WD repeat-containing protein 1</td><td>-0.49</td><td>0.019</td><td>0.542 *</td><td>✓</td></tr><tr><td>P60712</td><td>ACTB</td><td>Actin, cytoplasmic 1</td><td>-0.98</td><td>0.023</td><td>0.524 *</td><td>✓</td></tr><tr><td>A0JNJ5</td><td>MYL1</td><td>Myosin light chain 1/3, skeletal muscle isoform</td><td>-0.47</td><td>0.024</td><td>0.554 *</td><td>✓</td></tr><tr><td>P02453</td><td>COL1A1</td><td>Collagen alpha-1(I) chain</td><td>0.89</td><td>0.025</td><td>-0.536 *</td><td>✓</td></tr><tr><td>Q3SYZ8</td><td>PDLIM3</td><td>PDZ and LIM domain protein 3</td><td>-0.51</td><td>0.029</td><td>0.512 *</td><td>PDLIM7/ PDLIM1</td></tr><tr><td>Q0III9</td><td>ACTN3 <sup>c</sup></td><td>Alpha-actinin-3</td><td>-0.34</td><td>0.034</td><td>0.529 *</td><td>✓</td></tr><tr><td>A4FV78</td><td>KLHL41</td><td>KBTBD10 protein</td><td>-0.38</td><td>0.037</td><td></td><td>✓</td></tr><tr><td>F1N789</td><td>VCL</td><td>Vinculin</td><td>-0.21</td><td>0.040</td><td></td><td>✓</td></tr><tr><td>Q32LP2</td><td>RDX</td><td>Radixin</td><td>-0.26</td><td>0.043</td><td>0.476 *</td><td>✓</td></tr><tr><td colspan="7" style="text-align: left; font-weight: bold;">Energy metabolism (n = 5)</td></tr><tr><td>Q5E956</td><td>TP11</td><td>Triosephosphate isomerase</td><td>-0.35</td><td>0.018</td><td>0.631 **</td><td>✓</td></tr><tr><td>Q3ZBY4</td><td>ALDOC</td><td>Fructose-bisphosphate aldolase</td><td>-0.33</td><td>0.018</td><td>0.621 **</td><td>✓</td></tr><tr><td>A5D984</td><td>PKM</td><td>Pyruvate kinase</td><td>-0.39</td><td>0.029</td><td>0.540 *</td><td>✓</td></tr><tr><td>A6QLL8</td><td>ALDOA</td><td>Fructose-bisphosphate aldolase</td><td>-0.29</td><td>0.029</td><td>0.554 *</td><td>✓</td></tr><tr><td>A3KN12</td><td>ADSL</td><td>Adenylosuccinate lyase</td><td>0.27</td><td>0.019</td><td>-0.629 **</td><td></td></tr></tbody></table>
<table><caption>Table 2. *Cont.*</caption><thead><tr><th rowspan="2">Uniprot ID</th><th rowspan="2">Gene Name</th><th rowspan="2">Full Protein Name</th><th colspan="2">Differences</th><th>Pearson Correlations <sup>a</sup></th><th rowspan="2">Overlap with<br>Gagaoua et al.<br>Database [13]</th></tr><tr><th>Fold Change<br>(Log2)</th><th><i>p</i>-Value</th><th>WBSF</th></tr></thead><tbody><tr><td colspan="7"><strong>Heat shock proteins (n = 4)</strong></td></tr><tr><td>P19120</td><td>HSPA8 <sup>b</sup></td><td>Heat shock cognate 71 kDa protein</td><td>-0.41</td><td>0.006</td><td>0.590 **</td><td>✓</td></tr><tr><td>P31081</td><td>HSPD1</td><td>60 kDa heat shock protein, mitochondrial</td><td>0.84</td><td>0.016</td><td>-0.567 *</td><td></td></tr><tr><td>Q3T149</td><td>HSPB1</td><td>Heat shock protein beta-1</td><td>-0.44</td><td>0.038</td><td>0.531 *</td><td>✓</td></tr><tr><td>Q3ZBZ8</td><td>STIP1 <sup>b,c</sup></td><td>Stress-induced-phosphoprotein 1</td><td>-0.29</td><td>0.040</td><td>0.484 *</td><td>✓</td></tr><tr><td colspan="7"><strong>Oxidative stress (n = 2)</strong></td></tr><tr><td>P35705</td><td>PRDX3</td><td>Thioredoxin-dependent peroxide reductase</td><td>0.32</td><td>0.038</td><td>-0.488 *</td><td>PRDX6/PRDX1/PRDX2</td></tr><tr><td>Q5E946</td><td>PARK7</td><td>Protein/nucleic acid deglycase DJ-1</td><td>-0.47</td><td>0.046</td><td>0.501 *</td><td>✓</td></tr><tr><td colspan="7"><strong>Other pathways (n = 6)</strong></td></tr><tr><td>P11116</td><td>LGALS1</td><td>Galectin-1</td><td>-0.56</td><td>0.006</td><td>0.639 **</td><td>✓</td></tr><tr><td>E1BE77</td><td>TRIM72</td><td>Tripartite motif containing 72</td><td>-0.45</td><td>0.008</td><td>0.620 **</td><td>✓</td></tr><tr><td>P19879</td><td>OGN <sup>b</sup></td><td>Mimecan</td><td>-0.64</td><td>0.017</td><td>0.589 *</td><td></td></tr><tr><td>Q2HJF7</td><td>CAMK2D</td><td>Calcium/calmodulin-dependent protein kinase</td><td>-0.20</td><td>0.038</td><td>0.510 *</td><td></td></tr><tr><td>Q6EWQ7</td><td>EIF5A</td><td>Eukaryotic translation initiation factor 5A-1</td><td>-0.45</td><td>0.046</td><td>0.574 *</td><td></td></tr><tr><td>Q3SYR3</td><td>APOBEC2</td><td>Probable C->U-editing enzyme APOBEC-2</td><td>-0.66</td><td>0.047</td><td>0.530 *</td><td></td></tr></tbody></table>

The regression model built for WBSF is presented in Table 3. The model explained 79% of the variability in WBSF (*p* < 0.01), including the abundance of three proteins: MYOZ3 (Myozenin 3), BIN1 (Bridging Integrator-1), and OGN (Mimecan), which were all positively correlated with WBSF (negatively with tenderness). It should be highlighted that MYOZ3 alone explained 52% of the variability. In this model, the correlation of MYOZ3 with WBSF values is depicted in Figure 1.

<table><caption>Table 3. Best regression equation of WBSF based on the list of the significant differential proteins from Table 2.</caption><thead><tr><th>R-Squared <sup>a</sup></th><th>S.E</th><th>Entered Independent Variable <sup>b</sup></th><th>Partial R-Squared</th><th>Regression Coefficient</th><th>t-Value</th><th><i>p</i>-Value</th></tr></thead><tbody><tr><td rowspan="3">0.79 **</td><td>0.125</td><td>MYOZ3</td><td>0.52</td><td>0.486</td><td>3.875</td><td>0.002</td></tr><tr><td>0.116</td><td>BIN1</td><td>0.17</td><td>0.454</td><td>3.907</td><td>0.002</td></tr><tr><td>0.121</td><td>OGN</td><td>0.1</td><td>0.347</td><td>2.868</td><td>0.012</td></tr></tbody></table>

From the list of the putative protein biomarkers, 30 were negatively correlated with tenderness (positively with WBSF), from which MYOZ3, CFL2, and BIN1 were the most highly significantly correlated proteins. In addition, 4 proteins (COL1A1, ADSL, HSPD1, and PRDX3) were positively correlated with tenderness (negatively with WBSF; Figure 2a and Table 2). From the correlation analyses, Myozenin (MYOZ3) was strongly and significantly correlated with WBSF (Figure 1). No significant correlation was found between WBSF with KLHL41 and VCL (Table 2).

<footnote><sup>a</sup> Significance of the correlations: * *p* < 0.05; ** *p* < 0.01; *** *p* < 0.001. <sup>b</sup> Proteins identified as Quantitative Trait Loci (QTL) of shear force using ProteQTL tool included in ProteINSIDE (http://www.proteinside.org/, accessed on 28 November 2020) from the Animal QTL Database (https://www.animalgenome.org/QLLdb/, accessed on 28 November 2020). <sup>c</sup> Proteins identified as QTL of sensory tenderness.</footnote>

<footnote><sup>a</sup> Significance of the models: ** *p* < 0.01. <sup>b</sup> Variables are shown in order of their entrance, in a stepwise manner, in the regression model.</footnote>
<figure><img src="image_24.png" /><figcaption>**Figure 1.** Example of significant correlations between the abundance of Myozenin 3 (MYOZ3) and WBSF values. The tender samples are shown by triangles (Δ) and the tough samples by circles (○). The R-square of the correlation is given.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 2.** Bioinformatics and statistical analyses of the proteins identified to be differential between the tough and tender *Longissimus thoracis* muscle steaks. (a) Volcano plot of the differential proteins in terms of their abundance, with a total of 34 proteins that were significantly different between the two tenderness groups shown in red (negative direction with tenderness) and green colour (positive direction with tenderness). The other proteins that had a tendency or were not significant were in grey and black colour, respectively. (b) Networks of pathways and process enrichment cluster analysis based on the 34 differential proteins using Metascape®(https://metascape.org/, accessed on 28 November 2020). (c) Functional enrichment analysis based on the list of significant 17 Gene Ontology (GO) terms ranked by their *p*-value.</figcaption></figure>
## 3.2. Partial Least Squares to Explain the Variability of WBSF Values

Based on the VIP filter (Table 4), the WBSF PLS regression model retained 32 proteins, of which 16 proteins (MYOZ3, FHOD1, CFL2, PDLIM7, TMOD4, LGALS1, TPI1, ADSL, ALDOC, TRIM72, BIN1, MYLPF, CORO6, HSPA8, OGN, EIF5A) had a VIP > 1.0. The other 16 proteins (HSPD1, MYL1, ALDOA, WDR1, PKM, COL1A1, HSPB1, APOBEC2, ACTN3, ACTB, PDLIM3, CAMK2D, PARK7, PRDX3, STIP1, RDX) had a VIP between 0.8 and 1.0; KLHL41 and VCL were the only two proteins whose VIP values were under 0.8. Combined with the results of the correlation analyses, the 32 proteins were identified as related to WBSF regardless of the statistical method. In addition, MYOZ3 was the first ranked protein, with the highest VIP.

<table><caption>**Table 4.** Partial Least Squares (PLS) prediction of beef tenderness (WBSF) using the list of the 34 putative protein biomarkers based on their variable importance in the projection (VIP).</caption><thead><tr><th>Proteins</th><th>VIP</th><th>Direction (+ or -)</th></tr></thead><tbody><tr><td>MYOZ3: Myozenin 3</td><td>1.291</td><td>-</td></tr><tr><td>FHOD1: Formin homology 2 domain containing 1</td><td>1.243</td><td>-</td></tr><tr><td>CFL2: Cofilin-2</td><td>1.168</td><td>-</td></tr><tr><td>PDLIM7: PDZ and LIM domain protein 7</td><td>1.119</td><td>-</td></tr><tr><td>TMOD4: Tropomodulin-4</td><td>1.115</td><td>-</td></tr><tr><td>LGALS1: Galectin-1</td><td>1.113</td><td>-</td></tr><tr><td>TPI1: Triosephosphate isomerase</td><td>1.099</td><td>-</td></tr><tr><td>ADSL: Adenylosuccinate lyase</td><td>1.095</td><td>+</td></tr><tr><td>ALDOC: Fructose-bisphosphate aldolase</td><td>1.082</td><td>-</td></tr><tr><td>TRIM72: Tripartite motif containing 72</td><td>1.080</td><td>-</td></tr><tr><td>BIN1: Bridging Integrator-1</td><td>1.073</td><td>-</td></tr><tr><td>MYLPF: Myosin regulatory light chain 2, skeletal muscle isoform</td><td>1.067</td><td>-</td></tr><tr><td>CORO6: Coronin</td><td>1.056</td><td>-</td></tr><tr><td>HSPA8: Heat shock cognate 71 kDa protein</td><td>1.027</td><td>-</td></tr><tr><td>OGN: Mimercan</td><td>1.025</td><td>-</td></tr><tr><td>EIF5A: Eukaryotic translation initiation factor 5A-1</td><td>1.000</td><td>-</td></tr><tr><td>HSPD1: 60 kDa heat shock protein, mitochondrial</td><td>0.987</td><td>+</td></tr><tr><td>MYL1: Myosin light chain 1/3, skeletal muscle isoform</td><td>0.965</td><td>-</td></tr><tr><td>ALDOA: Fructose-bisphosphate aldolase</td><td>0.965</td><td>-</td></tr><tr><td>WDR1: WD repeat-containing protein 1</td><td>0.943</td><td>-</td></tr><tr><td>PKM: Pyruvate kinase</td><td>0.940</td><td>-</td></tr><tr><td>COL1A1: Collagen alpha-1(I) chain</td><td>0.934</td><td>+</td></tr><tr><td>HSPB1: Heat shock protein beta-1</td><td>0.924</td><td>-</td></tr><tr><td>APOBEC2: Probable C-&gt;U-editing enzyme APOBEC-2</td><td>0.923</td><td>-</td></tr><tr><td>ACTN3: Alpha-actinin-3</td><td>0.922</td><td>-</td></tr><tr><td>ACTB: Actin, cytoplasmic 1</td><td>0.913</td><td>-</td></tr><tr><td>PDLIM3: PDZ and LIM domain protein 3</td><td>0.892</td><td>-</td></tr><tr><td>CAMK2D: Calcium/calmodulin-dependent protein kinase type II subunit delta</td><td>0.889</td><td>-</td></tr><tr><td>PARK7: Protein/nucleic acid deglycase DJ-1</td><td>0.873</td><td>-</td></tr><tr><td>PRDX3: Thioredoxin-dependent peroxide reductase, mitochondrial</td><td>0.850</td><td>+</td></tr><tr><td>STIP1: Stress-induced-phosphoprotein 1</td><td>0.843</td><td>-</td></tr><tr><td>RDX: Radixin</td><td>0.830</td><td>-</td></tr><tr><td>KLHL41: KBTBD10 protein</td><td>0.771</td><td>-</td></tr><tr><td>VCL: Vinculin</td><td>0.708</td><td>-</td></tr></tbody></table>

## 3.3. Protein-Protein Interactions (PPI)

The protein-protein interaction network highlighted the importance of the structural and contractile pathways in beef tenderisation (Figure 3). In the network, ACTB (Actin) had the most interactions with other pathways, including energy metabolism and heat shock proteins, while ACTN3 (Alpha-actinin-3) and MYLPF (Myosin regulatory light chain 2) had
more involvement within the muscle structure pathway. The next most dominant pathways were cellular processes, binding, apoptosis, and transport proteins, which showed multiple interactions with the energy metabolism (TPI1, ALDOA, and PKM), heat shock proteins (HSPD1), and muscle contraction (PDLIM3). It should be noted that the proteins in the heat shock pathway had a close interaction with each other.

<figure><img src="image_26.png" /><figcaption>**Figure 3.** Protein-Protein interaction network built using the 34 differentially abundant proteins. The interaction map was generated from the web-based search STRING database (https://string-db.org/, accessed on 28 November 2020). Default settings of confidence of 0.6 and 4 criteria for linkage: co-occurrence, experimental evidence, existing databases, and text mining were used.</figcaption></figure>

### 3.4. Pathway and Process Enrichment Analysis

The Gene Ontology (GO) results are given in Table 5. Canonical glycolysis (GO:0061621), glycolytic process (GO:0006096), and muscle contraction (GO:0006936) were the top three Gene-Ontology (GO)-enriched terms identified from the list of the 34 differential proteins (Table 5), while Cellular Component (CC), cytosol (GO:0005829), extracellular exosome (GO:0070062), and cytoplasm (GO:0005737) were the most important three CC terms. It should be noted that a considerable number of proteins were classified as proteins binding (GO:0005515) in molecular function (Table 5). From the Metascape analysis, 17 top and significantly enriched terms were validated and allowed to construct process enrichment networks of the pathways (Figure 2b,c). The top six enriched term clusters were highlighted, including supramolecular fibre organisation (GO:0097435), muscle contraction (GO:0006936), muscle structure development (GO:0061061), glucose catabolic process (GO:0006007), striated muscle contraction (GO:0006941), and homotypic cell-cell adhesion (GO:0034109). The most dominant pathway was supramolecular fibre organisation and muscle contraction, which was consistent with the PPI data confirming their pivotal role in beef tenderisation of young Limousin bull beef (Figure 3).
<table><caption>**Table 5.** Top20 Gene Ontology (GO) terms computed using the list of the 34 putative protein biomarkers.</caption><thead><tr><th>GO</th><th>Function</th><th>Gene Name</th><th>GO Frequency within the Dataset (%)</th><th>GO Frequency within the Genome (%)</th><th>p-Values</th></tr></thead><tbody><tr><td colspan="6" style="text-align: center; font-weight: bold;">Biological Process (BP)</td></tr><tr><td>GO:0061621</td><td>canonical glycolysis</td><td>ALDOC TPI1 PKM<br>ALDOA</td><td>11.76</td><td>14.81</td><td>2.13 × 10<sup>-9</sup></td></tr><tr><td>GO:0006096</td><td>glycolytic process</td><td>ALDOC ALDOA PKM<br>TPI1</td><td>11.76</td><td>10.26</td><td>5.4 × 10<sup>-9</sup></td></tr><tr><td>GO:0006936</td><td>muscle contraction</td><td>TRIM72 MYL1 VCL<br>MYLPF TMOD4</td><td>14.71</td><td>2.35</td><td>2.86 × 10<sup>-8</sup></td></tr><tr><td>GO:0043312</td><td>neutrophil degranulation</td><td>VCL HSPA8 ALDOA<br>ALDOC PKM</td><td>14.71</td><td>1.03</td><td>1.28 × 10<sup>-6</sup></td></tr><tr><td>GO:0070527</td><td>platelet aggregation</td><td>ACTB VCL HSPB1</td><td>8.82</td><td>7.14</td><td>2.06 × 10<sup>-6</sup></td></tr><tr><td>GO:0006094</td><td>gluconeogenesis<br>response to unfolded<br>protein<br>maintenance of<br>blood-brain barrier</td><td>ALDOC TPI1 ALDOA<br>HSPB1 HSPA8 HSPD1</td><td>8.82</td><td>6.82</td><td>2.27 × 10<sup>-6</sup></td></tr><tr><td>GO:0035633</td><td>fructose 1,6-bisphosphate<br>metabolic process<br>actin filament<br>depolymerisation</td><td>VCL ACTB</td><td>5.88</td><td>66.67</td><td>5.99 × 10<sup>-6</sup></td></tr><tr><td>GO:0030388</td><td>apical junction assembly</td><td>ALDOA ALDOC</td><td>5.88</td><td>28.57</td><td>1.98 × 10<sup>-5</sup></td></tr><tr><td>GO:0030042</td><td>platelet degranulation</td><td>WDR1 CFL2</td><td>5.88</td><td>25</td><td>2.34 × 10<sup>-5</sup></td></tr><tr><td>GO:0002576</td><td>positive regulation of actin<br>filament depolymerisation</td><td>VCL WDR1<br>ALDOA WDR1 VCL</td><td>5.88<br>8.82</td><td>25<br>2.44</td><td>2.34 × 10<sup>-5</sup><br>3.11 × 10<sup>-5</sup></td></tr><tr><td>GO:0030836</td><td>fructose metabolic process</td><td>CFL2 WDR1</td><td>5.88</td><td>15.38</td><td>4.82 × 10<sup>-5</sup></td></tr><tr><td>GO:0006000</td><td>actin filament organisation</td><td>ALDOA ALDOC</td><td>5.88</td><td>13.33</td><td>5.97 × 10<sup>-5</sup></td></tr><tr><td>GO:0007015</td><td>protein refolding</td><td>ALDOA CORO6 TMOD4<br>HSPD1 HSPA8</td><td>8.82<br>5.88</td><td>1.54<br>9.52</td><td>9.6 × 10<sup>-5</sup><br>9.98 × 10<sup>-5</sup></td></tr><tr><td>GO:0042026</td><td>myofibril assembly</td><td>KLHL41 TMOD4</td><td>5.88</td><td>7.14</td><td>0.000162</td></tr><tr><td>GO:0030239</td><td>adherens junction<br>assembly</td><td>ACTB VCL</td><td>5.88</td><td>5.88</td><td>0.000221</td></tr><tr><td>GO:0043066</td><td>negative regulation of<br>apoptotic process</td><td>HSPD1 HSPB1 PRDX3<br>PARK7</td><td>11.76</td><td>0.49</td><td>0.000221</td></tr><tr><td>GO:0086091</td><td>regulation of heart rate by<br>cardiac conduction</td><td>BIN1 CAMK2D</td><td>5.88</td><td>5.56</td><td>0.000239</td></tr><tr><td colspan="6" style="text-align: center; font-weight: bold;">Cellular Component (CC)</td></tr><tr><td>GO:0005829</td><td>cytosol</td><td>VCL TPI1 MYLPF HSPA8<br>STIP1 EIF5A CAMK2D<br>ADSL MYL1 ACTN3<br>ALDOC HSPD1 BIN1<br>PRDX3 PARK7 PDLIM3<br>ALDOA HSPB1 WDR1<br>PKM FHOD1 PDLIM7<br>ACTB KLHL41<br>ALDOC HSPD1 ACTN3<br>ACTB PARK7 TPI1 PKM<br>VCL RDX CFL2 LGALS1<br>OGN ALDOA HSPA8<br>HSPB1 WDR1</td><td>70.59</td><td>0.5</td><td>3.03 × 10<sup>-21</sup></td></tr><tr><td>GO:0070062</td><td>extracellular exosome</td><td>PARK7 FHOD1 COL1A1<br>EIF5A PKM BIN1 HSPD1<br>HSPA8 KLHL41 PRDX3<br>TRIM72 APOBEC2<br>CAMK2D ACTB HSPB1<br>CFL2 LGALS1</td><td>47.06</td><td>0.58</td><td>1.67 × 10<sup>-15</sup></td></tr><tr><td>GO:0005737</td><td>cytoplasm</td><td></td><td>50</td><td>0.41</td><td>3 × 10<sup>-14</sup></td></tr></tbody></table>
<table><caption>Table 5. *Cont.*</caption><thead><tr><th>GO</th><th>Function</th><th>Gene Name</th><th>GO Frequency within the Dataset (%)</th><th>GO Frequency within the Genome (%)</th><th>p-Values</th></tr></thead><tbody><tr><td>GO:0005615</td><td>extracellular space</td><td>HSPD1 RDX ALDOA CFL2 HSPA8 TPI1 OGN COL1A1 HSPB1 LGALS1 ACTB VCL HSPB1 RDX HSPA8 ACTN3 ACTB PDLIM7 CFL2 MYOZ3 BIN1</td><td>32.35</td><td>0.77</td><td>3 × 10<sup>-12</sup></td></tr><tr><td>GO:0005925</td><td>focal adhesion</td><td>PDLIM7 ACTB ALDOA VCL KLHL41 HSPB1 BIN1</td><td>20.59</td><td>1.82</td><td>2.33 × 10<sup>-10</sup></td></tr><tr><td>GO:0015629</td><td>actin cytoskeleton</td><td>TMOD4 ACTB FHOD1 ALDOC PDLIM7 MYOZ3 BIN1 CFL2 PDLIM3 ALDOC HSPA8 PKM ALDOA VCL STIP1 CAMK2D EIF5A</td><td>17.65</td><td>3.14</td><td>2.54 × 10<sup>-10</sup></td></tr><tr><td>GO:0005856</td><td>cytoskeleton</td><td>BIN1 FHOD1 ALDOA PKM HSPB1 APOBEC2 ACTB HSPA8 PARK7 TPI1 PKM ALDOC WDR1</td><td>23.53</td><td>1.08</td><td>3.05 × 10<sup>-10</sup></td></tr><tr><td>GO:0030018</td><td>Z disc</td><td>LGALS1 ALDOA HSPA8 VCL COL1A1 OGN ALDOC PKM VCL ALDOA HSPA8</td><td>14.71</td><td>4.2</td><td>2.57 × 10<sup>-9</sup></td></tr><tr><td>GO:1904813</td><td>ficolin-1-rich granule lumen</td><td>ALDOC HSPA8 PKM ALDOA VCL STIP1 CAMK2D EIF5A BIN1 FHOD1 ALDOA PKM HSPB1 APOBEC2 ACTB HSPA8 PARK7 TPI1 PKM ALDOC WDR1</td><td>14.71</td><td>4.03</td><td>2.92 × 10<sup>-9</sup></td></tr><tr><td>GO:0005634</td><td>nucleus</td><td>LGALS1 ALDOA HSPA8 VCL COL1A1 OGN ALDOC PKM VCL ALDOA HSPA8</td><td>38.24</td><td>0.26</td><td>4.61 × 10<sup>-9</sup></td></tr><tr><td>GO:0005576</td><td>extracellular region</td><td>ALDOC PKM VCL ALDOA HSPA8</td><td>26.47</td><td>0.49</td><td>7.91 × 10<sup>-9</sup></td></tr><tr><td>GO:0034774</td><td>secretory granule lumen</td><td>ALDOC PKM VCL ALDOA HSPA8</td><td>14.71</td><td>1.56</td><td>2 × 10<sup>-7</sup></td></tr><tr><td>GO:0031674</td><td>I band</td><td>ALDOC CFL2 BIN1 PARK7 ACTB PDLIM3 PDLIM7 VCL</td><td>8.82</td><td>13.64</td><td>3.96 × 10<sup>-7</sup></td></tr><tr><td>GO:0005912</td><td>adherens junction</td><td>PARK7 ACTB PDLIM3 PDLIM7 VCL</td><td>14.71</td><td>1.04</td><td>1.28 × 10<sup>-6</sup></td></tr><tr><td>GO:0030864</td><td>cortical actin cytoskeleton</td><td>RDX CFL2 WDR1</td><td>8.82</td><td>8.57</td><td>1.28 × 10<sup>-6</sup></td></tr><tr><td>GO:0001725</td><td>stress fibre</td><td>PDLIM3 PDLIM7 FHOD1</td><td>8.82</td><td>6.25</td><td>2.62 × 10<sup>-6</sup></td></tr><tr><td>GO:0030424</td><td>axon</td><td>PARK7 HSPA8 BIN1 ACTB PARK7 CAMK2D ACTB</td><td>11.76</td><td>1.65</td><td>3.52 × 10<sup>-6</sup></td></tr><tr><td>GO:0005654</td><td>nucleoplasm</td><td>KLHL41 HSPA8 FHOD1 PDLIM7 STIP1 HSPA8</td><td>20.59</td><td>0.24</td><td>3.81 × 10<sup>-5</sup></td></tr><tr><td>GO:0101031</td><td>chaperone complex</td><td>VCL HSPA8 BIN1 HSPD1 ACTB RDX WDR1 KLHL41</td><td>5.88</td><td>15.38</td><td>4.82 × 10<sup>-5</sup></td></tr><tr><td>GO:0005886</td><td>plasma membrane</td><td>VCL HSPA8 BIN1 HSPD1 ACTB RDX WDR1 KLHL41</td><td>23.53</td><td>0.18</td><td>5.47 × 10<sup>-5</sup></td></tr><tr><td colspan="6" style="text-align: center; font-weight: bold;">Molecular Function (MF)</td></tr><tr><td>GO:0005515</td><td>protein binding</td><td>PRDX3 PARK7 ALDOA OGN MYOZ3 ACTN3 FHOD1 CAMK2D PKM HSPD1 TRIM72 COL1A1 CFL2 TMOD4 STIP1 ALDOC CORO6 BIN1 LGALS1 VCL EIF5A HSPA8 HSPB1 PDLIM3 TPI1 KLHL41 RDX ACTB PDLIM7 ALDOC PRDX3 TRIM72 PARK7 ACTB CAMK2D</td><td>85.29</td><td>0.45</td><td>9.5 × 10<sup>-24</sup></td></tr><tr><td>GO:0042802</td><td>identical protein binding</td><td>ADSL BIN1 HSPB1 APOBEC2 ACTN3 FHOD1 COL1A1 HSPB1 HSPD1 LGALS1 EIF5A APOBEC2 HSPA8 PKM ALDOA STIP1 RDX</td><td>38.24</td><td>0.92</td><td>2.98 × 10<sup>-15</sup></td></tr><tr><td>GO:0003723</td><td>RNA binding</td><td>ADSL BIN1 HSPB1 APOBEC2 ACTN3 FHOD1 COL1A1 HSPB1 HSPD1 LGALS1 EIF5A APOBEC2 HSPA8 PKM ALDOA STIP1 RDX</td><td>29.41</td><td>0.62</td><td>2.33 × 10<sup>-10</sup></td></tr></tbody></table>
<table><caption>Table 5. *Cont.*</caption><thead><tr><th>GO</th><th>Function</th><th>Gene Name</th><th>GO Frequency within the Dataset (%)</th><th>GO Frequency within the Genome (%)</th><th>p-Values</th></tr></thead><tbody><tr><td>GO:0045296</td><td>cadherin binding</td><td>PARK7 HSPA8 VCL RDX ALDOA PKM</td><td>17.65</td><td>2.03</td><td>2.44 × 10<sup>-9</sup></td></tr><tr><td>GO:0051015</td><td>actin filament binding</td><td>BIN1 FHOD1 WDR1 CFL2 CORO6</td><td>14.71</td><td>3.57</td><td>4.61 × 10<sup>-9</sup></td></tr><tr><td>GO:0003779</td><td>actin binding</td><td>PDLIM3 VCL ALDOA PDLIM7 MYOZ3 RDX</td><td>17.65</td><td>1.48</td><td>1.01 × 10<sup>-8</sup></td></tr><tr><td>GO:0008307</td><td>structural constituent of muscle</td><td>ACTN3 MYL1 MYLPF</td><td>8.82</td><td>6.52</td><td>2.48 × 10<sup>-6</sup></td></tr><tr><td>GO:0004332</td><td>fructose-bisphosphate aldolase activity</td><td>ALDOA ALDOC</td><td>5.88</td><td>66.67</td><td>5.99 × 10<sup>-6</sup></td></tr><tr><td>GO:0031625</td><td>ubiquitin protein ligase binding</td><td>HSPA8 HSPD1 VCL TPI1</td><td>11.76</td><td>1.36</td><td>6.77 × 10<sup>-6</sup></td></tr><tr><td>GO:0051087</td><td>chaperone binding</td><td>BIN1 HSPD1 HSPA8</td><td>8.82</td><td>3.3</td><td>1.44 × 10<sup>-5</sup></td></tr><tr><td>GO:0048156</td><td>tau protein binding</td><td>BIN1 ACTB</td><td>5.88</td><td>20</td><td>3.26 × 10<sup>-5</sup></td></tr><tr><td>GO:0023026</td><td>MHC class II protein complex binding</td><td>PKM HSPA8</td><td>5.88</td><td>11.76</td><td>7.34 × 10<sup>-5</sup></td></tr><tr><td>GO:0051371</td><td>muscle alpha-actinin binding</td><td>PDLIM3 PDLIM7</td><td>5.88</td><td>11.11</td><td>7.98 × 10<sup>-5</sup></td></tr><tr><td>GO:0044183</td><td>protein folding chaperone</td><td>HSPA8 HSPB1</td><td>5.88</td><td>7.14</td><td>0.000162</td></tr><tr><td>GO:0042803</td><td>protein homodimerisation activity</td><td>PARK7 CAMK2D TPI1 HSPB1</td><td>11.76</td><td>0.53</td><td>0.000173</td></tr><tr><td>GO:0003697</td><td>single-stranded DNA binding</td><td>HSPD1 PARK7</td><td>5.88</td><td>2.06</td><td>0.001335</td></tr><tr><td>GO:0019901</td><td>protein kinase binding</td><td>HSPB1 ACTB PRDX3</td><td>8.82</td><td>0.54</td><td>0.001393</td></tr><tr><td>GO:0044325</td><td>ion channel binding</td><td>ACTN3 CAMK2D</td><td>5.88</td><td>1.82</td><td>0.001617</td></tr><tr><td>GO:0002020</td><td>protease binding (S)-2-(5-amino-1-(5-phospho-D-ribosyl)imidazole-4-carboxamido)succinate</td><td>COL1A1 BIN1</td><td>5.88</td><td>1.61</td><td>0.001708</td></tr><tr><td>GO:0070626</td><td>AMP-lyase (fumarate-forming) activity</td><td>ADSL</td><td>2.94</td><td>100</td><td>0.001708</td></tr></tbody></table>

## 4. Discussion

The beef industry is consistently confronted with challenges in supplying beef with consistent eating qualities. Tenderness is one of the most important palatability traits of beef that affects the repurchase decisions of consumers. The pathways underpinning beef tenderness determination are complex and not fully elucidated, although a recent integromics meta-analysis by Gagaoua et al. [13] on the molecular signatures shed light on some of them. Thus, it was valuable to identify putative protein biomarkers of beef tenderness from two tenderness groups with a strong difference in shear force: tender (33.21 N) vs. tough (63.96 N; Table 1).

This study on Irish Limousin-cross cattle allowed us to get more insights and validate the association of certain proteins with tenderness and propose new ones that will further increase our knowledge and progress in the pipeline of beef tenderness discovery of biomarkers. This study allowed us also to (i) propose preliminary explanatory models of tenderness using multiple regression and partial least squares; (ii) compare the list of putative protein biomarkers identified in this trial with previous studies to verify the robustness of the discovered proteins; and finally, (iii) increase our knowledge on the biological pathways involved in the variation of beef tenderness evaluated in this study using WBSF at an end-point cooking temperature of 71 °C. The relationship between tenderness and the list of candidate proteins was discussed in the following sections.
## 4.1. *The Best Explanatory Proteins in the Regression Model of WBSF*

The best regression model built with MYOZ3, BIN1, and OGN proteins explained 79% of the observed variability in WBSF (p < 0.01). MYOZ3 is mainly expressed in skeletal muscle and enriched in fast-twitch muscle fibres. MYOZ3 belongs to the myzenen family, of which three other members were previously proposed as tenderness biomarkers [13]. MYOZ3 acts as an intracellular binding protein to link with Z-disc proteins such as alpha-actinin and gamma-filamin and transmit calcineurin signalling to the sarcomere [27]. Due to the capacity to bind multiple proteins, the relationship between MYOZ3 and meat quality, specifically tenderness, could be through regulating the Z-disc structure and signal transduction, influencing muscle fibre differentiation [28]. Consistent with our findings, a previous study reported a negative association between MYOZ3 and the shear force of *M. longissimus thoracis* in heifers [29].

BIN1, also known as Bridging Integrator-1, was identified in the present study for the first time to have a potential association with beef tenderness. BIN1 plays an important role in the regulation of endocytosis and has other roles as a central regulator of cell proliferation and apoptosis [30]. While no evidence was present in the literature on a specific relationship with tenderness, BIN1 was associated with another important beef production attribute, residual feed intake [31], which was previously associated with meat quality. OGN, which is also called mimecan, belongs to a secreted protein family of small leucine-rich proteoglycans located in the extracellular matrix [32]. OGN was negatively correlated with beef tenderness, and both MYOZ3 and OGN genes were located within a Quantitative Trait Loci (QTL) for shear force on chromosome 8 (Table 2). Interestingly, when protein profiles were compared between Japanese Black cattle and Holstein cattle, a higher abundance of OGN protein (mimecan) was found in the Holstein breed known to have lower fat content [33]. An important function of OGN is in collagen fibrillogenesis [32]. For this reason, it could be hypothesised that the greater abundance of OGN protein observed for tougher beef animals may be related to a higher abundance of connective tissue content in the muscle of tough beef [34], although we did not measure the connective tissue content in the present study.

## 4.2. *Dominant Pathway Related to WBSF of Young Limousin-Sired Bulls*

Muscle contraction and structure were identified as the most important pathway associated with WBSF in this study. Most of the proteins from this pathway were localised in the sarcomere. Of these, compared with the database of beef tenderness biomarkers of Gagaoua [13], 13 were already identified, and 4 proteins (BIN1, FHOD1, CORO6, and RDX) were reported for the first time in this study.

Myosin and actin were critically important to textural changes in muscle that occurred post-mortem during meat ageing through the weakening of the actin/myosin complex in the myofibril [22]. As the major component of the thick filaments of the myofibril, molecular myosin consisted of two heavy and four light chains. This study revealed, for example, MYLPF (Myosin regulatory light chain 2) and MYL1 (Myosin light chain 1/3) to be negative biomarkers of beef tenderness. Myosin light chains were wrapped around the head/rod junction of the myosin heavy chain in skeletal muscle myosin [35]. The MYLPF and MYL1 proteins were demonstrated [13] to correlate with beef tenderness; however, the direction of their relationships with this trait lacks consistency across studies; this phenomenon was well known to vary depending on the breed and muscle type [5,25,36], and was suggested to be related to post-translational modifications of the proteins [37]. Myosin light chain proteins were highly expressed in fast-twitch fibres. It was noteworthy that phosphorylation of MYL might play an essential role in proteolysis and onset of apoptosis in post-mortem muscle, which was favourable to the degradation of large molecules and final tenderisation of aged meat [38]. As for the second, most abundant myofibrillar protein in muscle, actin was interrelated with the apoptosis of the cytoskeleton during meat tenderisation [39]. In our study, ACTB, ACTN3, and TMOD4 were negatively correlated with beef tenderness, which is consistent with previous studies [9,29].
The collagen alpha-1(I) chain is an abundant connective tissue protein with an important function of support in the muscle tissue and bone in the body, and is encoded by *COL1A1* gene. Several studies showed a close relationship between collagen content and variation in meat tenderness [40,41]. Interestingly, in a previous study by Bjarnadóttir et al. [42], COL1A1 and COL1A2 were found to have lower abundance in tender beef muscle, which was opposite to our findings. However, there was also evidence of a positive relationship between COL1A1 abundance and intramuscular fat content, which could have an effect in promoting beef tenderness [43]. Thus, there was no consistent conclusion regarding the direct influence of COL1A1 on meat quality.

PDLIM3 and PDLIM7 are two members of the PDZ and LIM domain (PDLIM) family, participating in multifunctional protein-protein interaction, cytoskeleton, and signal transduction pathways [44]. PDLIM family proteins contain a PDZ domain in the N-terminal portion and the LIM domain in the C-terminal portion [45]. PDLIM1 and PDLIM7 were previously identified as negative biomarkers of beef tenderness [13], which was consistent with our results stating that PDLIM3 and PDLIM7 were positively correlated with WBSF and the negative relationship of this protein family with beef quality.

Of the putative biomarkers identified for the first time in this study, FHOD1 is an actin regulator which played an important role in the stabilisation of filamentous (F)-actin bundles by selectively covering and binding their barbed ends to actin filaments, thus protecting actin filaments in cytoskeletal structures [46]. Likewise, CORO6 is also an actin-binding protein that is mainly expressed in the heart and skeletal muscle [47]. RDX (Radixin) is referred to as a member of ERM (Ezrin/Radixin/Moesin) proteins that help maintain cytoskeletal organisation by binding specific membrane proteins to the actin cytoskeleton [48]. FHOD1, CORO6, and RDX were all positively correlated with WBSF (and negatively with tenderness), which could be related to their protective effect on the integrity of the cytoskeleton.

### 4.3. Candidate Protein Biomarkers of WBSF from the Energy Metabolism Pathway

Energy metabolism comprises a series of interconnected pathways that can function in the presence or absence of oxygen to generate adenosine triphosphate (ATP), which is an end-product of the processes of oxidative phosphorylation [49]. Of the five proteins identified from this pathway (TPI1, ALDOA, PKM, ADSL, and ALDOC), the first four proteins have been previously identified as putative biomarkers of beef tenderness. In this study, these proteins were all negatively correlated with beef tenderness except ADSL. TPI1 can catalyse the conversion of dihydroxyacetone phosphate to D-glyceraldehyde 3-phosphate, meanwhile, maintaining the equilibration of the triosephosphates produced by aldolase (ALDOA) [50]. Aldolase is an enzyme that catalyses the reversible conversion of fructose-1,6-bisphosphate to glyceraldehyde 3-phosphate and dihydroxyacetone phosphate [51]. ALDOA and ALDOC are two different isoformes of aldolase. In the literature, ALDOA was both positively and negatively correlated with beef tenderness depending on the gender and muscle fibre type [13], while ALDOC was first identified as a putative negative biomarker of tenderness in the present study. The relationship between aldolase and tenderness could be explained by its participation in muscle glycolysis, which, if variable, could alter the profile and extent of pH decline, thereby further influencing the integrity of the Z-line with consequences for beef tenderness [52].

Pyruvate kinase, also known as PKM, is an enzyme that catalyses the dephosphorylation of phosphoenolpyruvate to pyruvate, generating ATP and regulating cell metabolism during glycolysis [53]. PKM1 and PKM2 are the two predominant isoforms of PKM in skeletal muscles. PKM was listed in the repertoire of beef tenderness biomarkers in *longissimus* muscle [13], and our findings provided further corroboration for its role in beef tenderness determination.
### 4.4. Heat Shock Proteins (HSPs) as Important Indicators of WBSF

Heat shock proteins (HSPs) are a family of proteins that have as their main function the protection of the organism itself and its cellular structures in response to exposure to stressful conditions [54]. The current study showed a differential abundance of HSPA8, HSPD1, HSPB1, and STIP1, three of which were already discovered to play a role in the variability of tenderness. Interestingly, three of the HSPs identified here were from three different subfamilies of HSPs, i.e., the small HSPs (HSPB1), HSP70s (HSPA8), and HSP60s (HSPD1). Among those four proteins, HSPA8, HSPB1, and STIP1 showed higher abundance in tough meat while HSPD1 was in the opposite direction.

As one important member of the large HSP70 family, HSPA8 was identified in six previous studies to be related to beef tenderness, but the mechanistic connection with tenderness was not clear because the protein was sometimes positively correlated and sometimes negatively correlated with beef tenderness [13]. The impact of HSP70 proteins on meat tenderness was thought to be mainly because they obstruct pro-apoptotic factors such as Bcl-2 in apoptotic pathways [13,55]. HSPA8 is grouped in the response to unfolded protein (GO:0006986) and protein refolding (GO:0042026) in Table 5. In this sense, HSPA8 played an important role in response to cellular stress [56]. Moreover, this protective role might be based on its interaction with structural proteins or by regulating cell signalling pathways (Figure 3). STIP1, known as a stress-induced-phosphoprotein, is a co-chaperone whose negative relationship with tenderness, found here, was consistent with the findings of Picard and Gagaoua [8].

As for the small HSPs, HSPB1 was identified as a robust biomarker of beef tenderness (referring to the database by Gagaoua et al. [13]), and it was, from that integromics study, in the top five biomarkers of beef tenderness from a list of 124 proteins. Extrinsic stressors, such as pre-slaughter or post-mortem management conditions, were sources of the intensive production of sHSPs in the muscle, which, like the larger HSPs, also play a regulatory role in delaying the apoptosis onset, the protection of myofibrillar proteins from proteolysis, and other cellular homeostasis roles [15,39]. The positive and negative relationships identified between sHSPs and tenderness might be due to interactions of factors such as animal type/breed, gender, muscle type, and pre-slaughter conditions [12,39,57]. In this study, HSPB1 was negatively correlated with beef tenderness, which would be consistent with its protective function against proteolysis in skeletal muscle.

It was notable that HSPD1, which is a member of the HSP60 family, was identified to be correlated with beef tenderness for the first time in the present study. Under stress conditions, the HSP60 family of proteins inside the mitochondrial matrix usually acts as molecular chaperones, collaborating with the co-chaperone Hsp10 to promote the correct folding of imported proteins and proper assembly of unfolded polypeptides [58]. A positive relationship between HSPD1 and tenderness might be hypothesised by its function in the energy metabolism pathway to maintain energy supply during proteolysis of myofibril proteins (Figure 3). As with HSPB1, HSPD1 was also associated with beef colour, which deepened our knowledge of the influential role of HSP proteins in post-mortem muscle events and consequences on meat quality [14].

### 4.5. Putative Biomarkers of Tenderness Related to Oxidative Stress

After slaughter, the lipid and protein fractions of muscle are targeted by various reactive oxygen species (ROS), causing structural alteration or denaturation of proteins [59]. In the context of oxidative stress, meat tenderness can be affected by cellular antioxidants, which include both enzymatic and non-enzymatic scavenger agents engaged in protecting the muscle proteins from damage by ROS, thereby further maintaining cell homeostasis. Meanwhile, it is noteworthy that meat tenderness could also be influenced by ROS damage produced by mitochondria, which play an important role in supplying energy during the conversion of muscle into meat [60,61].

In this study, two important proteins from the oxidative pathway, i.e., PARK7 and PRDX3 were identified and validated in comparison to the previous list of robust beef
tenderness biomarkers [13]. PARK7, also named DJ-1, was secreted from the cytosol to mitochondria to remove the mitochondrial H2O2 and maintain the integrity of the organelle in response to oxidative stress [62]. Consistent with the results of most previous studies, PARK7 level was negatively correlated with beef tenderness [8,63,64]. A mechanism could be deduced where PARK7 had an inhibitory effect on the pro-apoptotic factors and caspases (proteolytic enzymes) by interacting with other proteins from energy metabolism and HSPs pathways, as depicted by the PPI network (Figure 3), thus contributing to beef toughness by slowing or limiting post-mortem apoptosis of muscle cells [13,65].

PRDX3 is a member of the peroxiredoxins (Prxs), a ubiquitous family with six subgroups, and which, in bovine, contains six members [66]. PRDX3 is exclusively located in mitochondria with an oligomeric ring structure [67]. It should be highlighted that PRDX3 was previously identified to be positively related to beef tenderness in *Semimembranosus* muscle [68], which was consistent with our findings for *Longissimus thoracis* muscle. Regarding the possible mechanism, it could be assumed that the antioxidant enzyme PRDX3 could prevent the accumulation of ROS, protecting the function of proteases and the operation of the electron transport system, and thus, leading to apoptosis promotion and meat tenderisation. Other members of peroxiredoxins were also found to be associated with several beef quality traits, including PRDX1 [69] and PRDX2 [70] with tenderness and PRDX6 with tenderness [9], pH decline [71], and beef colour [25].

### 4.6. *Proteins from Other Pathways*

LGALS1 and TRIM72 were negative biomarkers of beef tenderness in this study, which was consistent with previous reports [8,9,42]. LGALS1 (Galectin-1) belongs to a family of β-galactoside-binding proteins, which may act as promoters of apoptosis and have an impact on cell proliferation and skeletal muscle differentiation [42]. However, the mechanism behind the association between Galectin-1 and meat tenderness was still obscure due to its complex functions under different conditions. As a signalling protein expressed in skeletal muscle, Tripartite motif-containing 72 (TRIM72) was considered as a sensor of oxidation on membrane damage [72]. TRIM72 may act as a scavenger of the harmful agents accumulated under the apoptotic process, leading to a limitation of apoptosis and tough meat [9]. In line with our findings, there was also a higher abundance of TRIM72 reported in tough beef, hence showing its negative role in the apoptotic pathway [68]. In addition, TRIM72 was first identified to correlate with beef colour, which confirmed the anti-oxidative properties of this protein, allowing it to be suggested as a relevant marker for multiple beef quality traits [63].

## 5. Conclusions

The results of this study allowed us to validate 23 putative biomarkers on Irish cattle (Limousin-sired bulls) and to propose 11 new proteins that increased our knowledge on the main biological pathways underpinning beef tenderness variation in the *Longissimus thoracis* muscle of young bulls. The network and gene ontology analyses allowed us to better characterise the enriched molecular pathways. This study also suggested a regression model with an R-squared of 79% using three proteins-MYOZ3, BIN1, and OGN-to explain the relationship between the abundance of these protein biomarkers and WBSF values. Further analyses would assess the robustness of the list of putative biomarkers identified in this study using accurate methods and new populations.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10.3390/foods10050952/s1, Table S1. Individual Warner-Bratzler shear force (WBSF) values of the *Longissimus thoracis* muscles from the 18 animals used in this trial.

**Author Contributions:** Conceptualisation, R.M.H., M.G., and A.M.M.; methodology, M.G., R.M.H., and Y.Z.; sampling, Y.Z., R.M.H., and J.C.; lab work, Y.Z. and M.G.; data curation, M.G., D.V., J.C. and Y.Z.; writing-original draft preparation, Y.Z. and M.G.; writing-review and editing, M.G., R.M.H., A.L.K., A.M.M., and T.S.; supervision, R.M.H., A.M.M., A.L.K., and M.G.; project administration,
R.M.H., A.M.M., and A.L.K.; funding acquisition, R.M.H., A.M.M. and T.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Teagasc and the Walsh Scholarship programme (Walsh Scholarship 2016022), the BreedQuality project (11/SF/311), which is supported by The Irish Department of Agriculture, Food, and the Marine (DAFM) under the National Development Plan 2007–2013. Mohammed Gagaoua is a Marie Skłodowska-Curie Career-FIT Fellow under project number MF20180029, grant agreement No. 713654.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data sharing is not applicable to this article.

**Acknowledgments:** The authors convey special thanks to David Sheehan and Dilip K. Rai for their help to Yao Zhu and for support in funding acquisition. We acknowledge the Irish Cattle Breeders Federation for access to samples.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Kerth, C.R.; Miller, R.K. Beef Flavor: A Review from Chemistry to Consumer. *J. Sci. Food Agric.* **2015**, *95*, 2783–2798. [CrossRef]

2. Purslow, P.P.; Gagaoua, M.; Warner, R.D. Insights on Meat Quality from Combining Traditional Studies and Proteomics. *Meat Sci.* **2021**, *174*, 108423. [CrossRef] [PubMed]

3. Henchion, M.; McCarthy, M.; Resconi, V.C.; Troy, D. Meat Consumption: Trends and Quality Matters. *Meat Sci.* **2014**, *98*, 561–568. [CrossRef] [PubMed]

4. Cuvelier, C.; Clinquart, A.; Hocquette, J.F.; Cabaraux, J.F.; Dufrasne, I.; Istasse, L.; Hornick, J.L. Comparison of Composition and Quality Traits of Meat from Young Finishing Bulls from Belgian Blue, Limousin and Aberdeen Angus Breeds. *Meat Sci.* **2006**, *74*, 522–531. [CrossRef] [PubMed]

5. Gagaoua, M.; Terlouw, C.; Richardson, I.; Hocquette, J.-F.; Picard, B. The Associations between Proteomic Biomarkers and Beef Tenderness Depend on the End-Point Cooking Temperature, the Country Origin of the Panelists and Breed. *Meat Sci.* **2019**, *157*, 107871. [CrossRef] [PubMed]

6. Cafferky, J.; Hamill, R.M.; Allen, P.; O'Doherty, J.V.; Cromie, A.; Sweeney, T. Effect of Breed and Gender on Meat Quality of M. Longissimus Thoracis et Lumborum Muscle from Crossbred Beef Bulls and Steers. *Foods* **2019**, *8*, 173. [CrossRef]

7. Nian, Y.; Kerry, J.P.; Prendiville, R.; Allen, P. The Eating Quality of Beef from Young Dairy Bulls Derived from Two Breed Types at Three Ages from Two Different Production Systems. *Ir. J. Agric. Food Res.* **2017**, *56*, 31–44. [CrossRef]

8. Picard, B.; Gagaoua, M. Muscle Fiber Properties in Cattle and Their Relationships with Meat Qualities: An Overview. *J. Agric. Food Chem.* **2020**, *68*, 6021–6039. [CrossRef] [PubMed]

9. Gagaoua, M.; Monteils, V.; Picard, B. Data from the Farmgate-to-Meat Continuum Including Omics-Based Biomarkers to Better Understand the Variability of Beef Tenderness: An Integromics Approach. *J. Agric. Food Chem.* **2018**, *66*, 13552–13563. [CrossRef]

10. Terlouw, E.M.C.; Picard, B.; Deiss, V.; Berri, C.; Hocquette, J.-F.; Lebret, B.; Lefèvre, F.; Hamill, R.; Gagaoua, M. Understanding the Determination of Meat Quality Using Biochemical Characteristics of the Muscle: Stress at Slaughter and Other Missing Keys. *Foods* **2021**, *10*, 84. [CrossRef]

11. Gagaoua, M.; Monteils, V.; Couvreur, S.; Picard, B. Identification of Biomarkers Associated with the Rearing Practices, Carcass Characteristics, and Beef Quality: An Integrative Approach. *J. Agric. Food Chem.* **2017**, *65*, 8264–8278. [CrossRef] [PubMed]

12. Picard, B.; Gagaoua, M. Meta-Proteomics for the Discovery of Protein Biomarkers of Beef Tenderness: An Overview of Integrated Studies. *Food Res. Int.* **2020**, *127*, 108739. [CrossRef] [PubMed]

13. Gagaoua, M.; Terlouw, E.M.C.; Mullen, A.M.; Franco, D.; Warner, R.D.; Lorenzo, J.M.; Purslow, P.P.; Gerrard, D.; Hopkins, D.L.; Troy, D.; et al. Molecular Signatures of Beef Tenderness: Underlying Mechanisms Based on Integromics of Protein Biomarkers from Multi-Platform Proteomics Studies. *Meat Sci.* **2021**, *172*, 108311. [CrossRef]

14. Gagaoua, M.; Bonnet, M.; Picard, B. Protein Array-Based Approach to Evaluate Biomarkers of Beef Tenderness and Marbling in Cows: Understanding of the Underlying Mechanisms and Prediction. *Foods* **2020**, *9*, 1180. [CrossRef]

15. Picard, B.; Gagaoua, M. Chapter 11—Proteomic Investigations of Beef Tenderness. In *Proteomics in Food Science*; Colgrave, M.L., Ed.; Academic Press: Cambridge, MA, USA, 2017; pp. 177–197. ISBN 978-0-12-804007-2.

16. Munekata, P.E.; Pateiro, M.; López-Pedrouso, M.; Gagaoua, M.; Lorenzo, J.M. Foodomics in Meat Quality. *Curr. Opin. Food Sci.* **2021**, *38*, 79–85. [CrossRef]

17. Gagaoua, M.; Hughes, J.; Terlouw, E.M.C.; Warner, R.D.; Purslow, P.P.; Lorenzo, J.M.; Picard, B. Proteomic Biomarkers of Beef Colour. *Trends Food Sci. Technol.* **2020**, *101*, 234–252. [CrossRef]

18. Zhu, Y.; Mullen, A.M.; Rai, D.K.; Kelly, A.L.; Sheehan, D.; Cafferky, J.; Hamill, R.M. Assessment of RNAlater® as a Potential Method to Preserve Bovine Muscle Proteins Compared with Dry Ice in a Proteomic Study. *Foods* **2019**, *8*, 60. [CrossRef]
19. Bouley, J.; Chambon, C.; Picard, B. Mapping of Bovine Skeletal Muscle Proteins Using Two-Dimensional Gel Electrophoresis and Mass Spectrometry. *Proteomics* **2004**, *4*, 1811–1824. [CrossRef]

20. Bradford, M.M. A Rapid and Sensitive Method for the Quantitation of Microgram Quantities of Protein Utilizing the Principle of Protein-Dye Binding. *Anal. Biochem.* **1976**, *72*, 248–254. [CrossRef]

21. Zhu, Y.; Gagaoua, M.; Mullen, A.M.; Viala, D.; Rai, D.K.; Kelly, A.L.; Sheehan, D.; Hamill, R.M. Shotgun Proteomics for the Preliminary Identification of Biomarkers of Beef Sensory Tenderness, Juiciness and Chewiness from Plasma and Muscle of Young Limousin-Sired Bulls. *Meat Sci.* **2021**, *176*, 108488. [CrossRef] [PubMed]

22. Gagaoua, M.; Troy, D.; Mullen, A.M. The Extent and Rate of the Appearance of the Major 110 and 30 KDa Proteolytic Fragments during Post-Mortem Aging of Beef Depend on the Glycolysing Rate of the Muscle and Aging Time: An LC-MS/MS Approach to Decipher Their Proteome and Associated Pathways. *J. Agric. Food Chem.* **2021**, *69*, 602–614. [CrossRef] [PubMed]

23. Kaspric, N.; Picard, B.; Reichstadt, M.; Tournayre, J.; Bonnet, M. ProteINSIDE to Easily Investigate Proteomics Data from Ruminants: Application to Mine Proteome of Adipose and Muscle Tissues in Bovine Foetuses. *PLoS ONE* **2015**, *10*, e0128086. [CrossRef]

24. Zhou, Y.; Zhou, B.; Pache, L.; Chang, M.; Khodabakhshi, A.H.; Tanaseichuk, O.; Benner, C.; Chanda, S.K. Metascape Provides a Biologist-Oriented Resource for the Analysis of Systems-Level Datasets. *Nat. Commun.* **2019**, *10*, 1523. [CrossRef]

25. Gagaoua, M.; Terlouw, E.M.C.; Picard, B. The Study of Protein Biomarkers to Understand the Biochemical Processes Underlying Beef Color Development in Young Bulls. *Meat Sci.* **2017**, *134*, 18–27. [CrossRef] [PubMed]

26. Huffman, K.L.; Miller, M.F.; Hoover, L.C.; Wu, C.K.; Brittin, H.C.; Ramsey, C.B. Effect of Beef Tenderness on Consumer Satisfaction with Steaks Consumed in the Home and Restaurant. *J. Anim Sci.* **1996**, *74*, 91–97. [CrossRef]

27. Ye, M.; Ye, F.; He, L.; Luo, B.; Yang, F.; Cui, C.; Zhao, X.; Yin, H.; Li, D.; Xu, H.; et al. Transcriptomic Analysis of Chicken Myozenin 3 Regulation Reveals Its Potential Role in Cell Proliferation. *PLoS ONE* **2017**, *12*, e0189476. [CrossRef]

28. Ye, M.; Ye, F.; He, L.; Liu, Y.; Zhao, X.; Yin, H.; Li, D.; Xu, H.; Zhu, Q.; Wang, Y. Molecular Cloning, Expression Profiling, and Marker Validation of the Chicken Myoz3 Gene. *Biomed. Res. Int.* **2017**, *2017*, 5930918. [CrossRef] [PubMed]

29. Boudon, S.; Ounaissi, D.; Viala, D.; Monteils, V.; Picard, B.; Cassar-Malek, I. Label Free Shotgun Proteomics for the Identification of Protein Biomarkers for Beef Tenderness in Muscle and Plasma of Heifers. *J. Proteom.* **2020**, *217*, 103685. [CrossRef] [PubMed]

30. Elliott, K.; Sakamuro, D.; Basu, A.; Du, W.; Wunner, W.; Staller, P.; Gaubatz, S.; Zhang, H.; Prochownik, E.; Eilers, M.; et al. Bin1 Functionally Interacts with Myc and Inhibits Cell Proliferation via Multiple Mechanisms. *Oncogene* **1999**, *18*, 3564–3573. [CrossRef]

31. Karisa, B.K.; Thomson, J.; Wang, Z.; Stothard, P.; Moore, S.S.; Plastow, G.S. Candidate Genes and Single Nucleotide Polymorphisms Associated with Variation in Residual Feed Intake in Beef Cattle1. *J. Anim. Sci.* **2013**, *91*, 3502–3513. [CrossRef] [PubMed]

32. Hu, S.-M.; Li, F.; Yu, H.-M.; Li, R.-Y.; Ma, Q.-Y.; Ye, T.-J.; Lu, Z.-Y.; Chen, J.-L.; Song, H.-D. The Mimecan Gene Expressed in Human Pituitary and Regulated by Pituitary Transcription Factor-1 as a Marker for Diagnosing Pituitary Tumors. *J. Clin. Endocrinol. Metab.* **2005**, *90*, 6657–6664. [CrossRef] [PubMed]

33. Ohsaki, H.; Okada, M.; Sasazaki, S.; Hinenoya, T.; Sawa, T.; Iwanaga, S.; Tsuruta, H.; Mukai, F.; Mannen, H. Proteomic Comparison between Japanese Black and Holstein Cattle by Two-Dimensional Gel Electrophoresis and Identification of Proteins. *Asian-Australas. J. Anim. Sci.* **2007**, *20*, 638–644. [CrossRef]

34. Tasheva, E.S.; Koester, A.; Paulsen, A.Q.; Garrett, A.S.; Boyle, D.L.; Davidson, H.J.; Song, M.; Fox, N.; Conrad, G.W. Mimecan/Osteoglycin-Deficient Mice Have Collagen Fibril Abnormalities. *Mol. Vis.* **2002**, *8*, 407–415. [PubMed]

35. Lowey, S.; Waller, G.S.; Trybus, K.M. Skeletal Muscle Myosin Light Chains Are Essential for Physiological Speeds of Shortening. *Nature* **1993**, *365*, 454–456. [CrossRef] [PubMed]

36. Picard, B.; Gagaoua, M.; Micol, D.; Cassar-Malek, I.; Hocquette, J.-F.; Terlouw, C.E.M. Inverse Relationships between Biomarkers and Beef Tenderness According to Contractile and Metabolic Properties of the Muscle. *J. Agric. Food Chem.* **2014**, *62*, 9808–9818. [CrossRef]

37. Mato, A.; Rodriguez-Vázquez, R.; López-Pedrouso, M.; Bravo, S.; Franco, D.; Zapata, C. The First Evidence of Global Meat Phosphoproteome Changes in Response to Pre-Slaughter Stress. *Bmc Genom.* **2019**, *20*, 590. [CrossRef] [PubMed]

38. De Rodrigues, R.T.S.; Chizzotti, M.L.; Vital, C.E.; Baracat-Pereira, M.C.; Barros, E.; Busato, K.C.; Gomes, R.A.; Ladeira, M.M.; Martins, T.D.S. Differences in Beef Quality between Angus (*Bos Taurus Taurus*) and Nellore (*Bos Taurus Indicus*) Cattle through a Proteomic and Phosphoproteomic Approach. *PLoS ONE* **2017**, *12*, e0170294. [CrossRef]

39. Ouali, A.; Gagaoua, M.; Boudida, Y.; Becila, S.; Boudjellal, A.; Herrera-Mendez, C.H.; Sentandreu, M.A. Biomarkers of Meat Tenderness: Present Knowledge and Perspectives in Regards to Our Current Understanding of the Mechanisms Involved. *Meat Sci.* **2013**, *95*, 854–870. [CrossRef] [PubMed]

40. Weston, A.R.; Rogers, R.W.; Althen, T.G. Review: The Role of Collagen in Meat Tenderness. *Prof. Anim. Sci.* **2002**, *18*, 107–111. [CrossRef]

41. Lepetit, J. Collagen Contribution to Meat Toughness: Theoretical Aspects. *Meat Sci.* **2008**, *80*, 960–967. [CrossRef]

42. Bjarnadóttir, S.G.; Hollung, K.; Høy, M.; Bendixen, E.; Codrea, M.C.; Veiseth-Kent, E. Changes in Protein Abundance between Tender and Tough Meat from Bovine Longissimus Thoracis Muscle Assessed by Isobaric Tag for Relative and Absolute Quantitation (ITRAQ) and 2-Dimensional Gel Electrophoresis Analysis1. *J. Anim. Sci.* **2012**, *90*, 2035–2043. [CrossRef]

43. Liao, H.; Zhang, X.H.; Qi, Y.X.; Wang, Y.Q.; Pang, Y.Z.; Liu, Z.B.Z.P. The Relationships of Collagen and ADAMTS2 Expressionlevels with Meat Quality Traits in Cattle. *Indian J. Anim. Res.* **2016**, *52*, 167–172. [CrossRef]
44. Cui, L.; Cheng, Z.; Hu, K.; Pang, Y.; Liu, Y.; Qian, T.; Quan, L.; Dai, Y.; Pang, Y.; Ye, X.; et al. Prognostic Value of the PDLIM Family in Acute Myeloid Leukemia. *Am. J. Transl. Res.* **2019**, *11*, 6124–6131.

45. Rios, H.; Paganelli, A.R.; Fosser, N.S. The Role of PDLIM1, a PDZ-LIM Domain Protein, at the Ribbon Synapses in the Chicken Retina. *J. Comp. Neurol.* **2020**, *528*, 1820–1832. [CrossRef]

46. Schönichen, A.; Mannherz, H.G.; Behrmann, E.; Mazur, A.J.; Kühn, S.; Silván, U.; Schoenenberger, C.-A.; Fackler, O.T.; Raunser, S.; Dehmelt, L.; et al. FHOD1 Is a Combined Actin Filament Capping and Bundling Factor That Selectively Associates with Actin Arcs and Stress Fibers. *J. Cell Sci.* **2013**, *126*, 1891–1901. [CrossRef]

47. Hemerich, D.; Pei, J.; Harakalova, M.; van Setten, J.; Boymans, S.; Boukens, B.J.; Efimov, I.R.; Michels, M.; van der Velden, J.; Vink, A.; et al. Integrative Functional Annotation of 52 Genetic Loci Influencing Myocardial Mass Identifies Candidate Regulatory Variants and Target Genes. *Circ. Genom Precis Med.* **2019**, *12*, e002328. [CrossRef]

48. Hansen, M.D.H.; Kwiatkowski, A.V. Chapter One—Control of Actin Dynamics by Allosteric Regulation of Actin Binding Proteins. In *International Review of Cell and Molecular Biology*; Jeon, K.W., Ed.; Academic Press: Cambridge, MA, USA, 2013; Volume 303, pp. 1–25.

49. Ferguson, D.M.; Gerrard, D.E. Regulation of Post-Mortem Glycolysis in Ruminant Muscle. *Anim. Prod. Sci.* **2014**, *54*, 464–481. [CrossRef]

50. Orosz, F.; Oláh, J.; Ovádi, J. Triosephosphate Isomerase Deficiency: New Insights into an Enigmatic Disease. *Biochim. Biophys. Acta (BBA) Mol. Basis Dis.* **2009**, *1792*, 1168–1174. [CrossRef]

51. Esposito, G.; Vitagliano, L.; Costanzo, P.; Borrelli, L.; Barone, R.; Pavone, L.; Izzo, P.; Zagari, A.; Salvatore, F. Human Aldolase A Natural Mutants: Relationship between Flexibility of the C-Terminal Region and Enzyme Function. *Biochem. J.* **2004**, *380*, 51–56. [CrossRef]

52. Hughes, J.; Clarke, F.; Li, Y.; Purslow, P.; Warner, R. Differences in Light Scattering between Pale and Dark Beef Longissimus Thoracis Muscles Are Primarily Caused by Differences in the Myofilament Lattice, Myofibril and Muscle Fibre Transverse Spacings. *Meat Sci.* **2019**, *149*, 96–106. [CrossRef]

53. Israelsen, W.J.; Vander Heiden, M.G. Pyruvate Kinase: Function, Regulation and Role in Cancer. *Semin. Cell Dev. Biol.* **2015**, *43*, 43–51. [CrossRef]

54. Daugaard, M.; Rohde, M.; Jäättelä, M. The Heat Shock Protein 70 Family: Highly Homologous Proteins with Overlapping and Distinct Functions. *Febes Lett.* **2007**, *581*, 3702–3710. [CrossRef]

55. Jiang, B.; Liang, P.; Deng, G.; Tu, Z.; Liu, M.; Xiao, X. Increased Stability of Bcl-2 in HSP70-Mediated Protection against Apoptosis Induced by Oxidative Stress. *Cell Stress Chaperones* **2011**, *16*, 143–152. [CrossRef]

56. Mayer, M.P. Hsp70 Chaperone Dynamics and Molecular Mechanism. *Trends Biochem. Sci.* **2013**, *38*, 507–514. [CrossRef]

57. Morzel, M.; Terlouw, C.; Chambon, C.; Micol, D.; Picard, B. Muscle Proteome and Meat Eating Qualities of Longissimus Thoracis of “Blonde d’Aquitaine” Young Bulls: A Central Role of HSP27 Isoforms. *Meat Sci.* **2008**, *78*, 297–304. [CrossRef]

58. Caruso Bavisotto, C.; Alberti, G.; Vitale, A.M.; Paladino, L.; Campanella, C.; Rappa, F.; Gorska, M.; Conway de Macario, E.; Cappello, F.; Macario, A.J.L.; et al. Hsp60 Post-Translational Modifications: Functional and Pathological Consequences. *Front. Mol. Biosci.* **2020**, 7. [CrossRef]

59. McDonagh, B.; Sheehan, D. Redox Proteomics in the Blue Mussel Mytilus Edulis: Carbonylation Is Not a Pre-Requisite for Ubiquitination in Acute Free Radical-Mediated Oxidative Stress. *Aquat. Toxicol.* **2006**, *79*, 325–333. [CrossRef]

60. Sierra, V.; Oliván, M. Role of Mitochondria on Muscle Cell Death and Meat Tenderization. *Recent Pat. Endocr. Metab. Immune Drug Discov.* **2013**. [CrossRef]

61. Lana, A.; Zolla, L. Apoptosis or Autophagy, That Is the Question: Two Ways for Muscle Sacrifice towards Meat. *Trends Food Sci. Technol.* **2015**, *46*, 231–241. [CrossRef]

62. Thomas, K.J.; McCoy, M.K.; Blackinton, J.; Beilina, A.; van der Brug, M.; Sandebring, A.; Miller, D.; Maric, D.; Cedazo-Minguez, A.; Cookson, M.R. DJ-1 Acts in Parallel to the PINK1/Parkin Pathway to Control Mitochondrial Function and Autophagy. *Hum. Mol. Genet.* **2011**, *20*, 40–50. [CrossRef]

63. Gagaoua, M.; Bonnet, M.; Ellies-Oury, M.-P.; De Koning, L.; Picard, B. Reverse Phase Protein Arrays for the Identification/Validation of Biomarkers of Beef Texture and Their Use for Early Classification of Carcasses. *Food Chem.* **2018**, *250*, 245–252. [CrossRef]

64. Jia, X.; Veiseth-Kent, E.; Grove, H.; Kuziora, P.; Aass, L.; Hildrum, K.I.; Hollung, K. Peroxiredoxin-6-a Potential Protein Marker for Meat Tenderness in Bovine Longissimus Thoracis Muscle. *J. Anim. Sci.* **2009**, *87*, 2391–2399. [CrossRef]

65. Fan, J.; Ren, H.; Jia, N.; Fei, E.; Zhou, T.; Jiang, P.; Wu, M.; Wang, G. DJ-1 Decreases Bax Expression through Repressing P53 Transcriptional Activity. *J. Biol. Chem.* **2008**, *283*, 4022–4030. [CrossRef] [PubMed]

66. Fisher, A.B. Peroxiredoxin 6: A Bifunctional Enzyme with Glutathione Peroxidase and Phospholipase A2 Activities. *Antioxid. Redox Signal.* **2010**, *15*, 831–844. [CrossRef]

67. Cao, Z.; Bhella, D.; Lindsay, J.G. Reconstitution of the Mitochondrial PrxIII Antioxidant Defence Pathway: General Properties and Factors Affecting PrxIII Activity and Oligomeric State. *J. Mol. Biol.* **2007**, *372*, 1022–1033. [CrossRef] [PubMed]

68. Grabež, V.; Kathri, M.; Phung, V.; Moe, K.M.; Slinde, E.; Skaugen, M.; Saarem, K.; Egelandsdal, B. Protein Expression and Oxygen Consumption Rate of Early Postmortem Mitochondria Relate to Meat Tenderness. *J. Anim. Sci.* **2015**, *93*, 1967–1979. [CrossRef]

69. Polati, R.; Menini, M.; Robotti, E.; Millioni, R.; Marengo, E.; Novelli, E.; Balzan, S.; Cecconi, D. Proteomic Changes Involved in Tenderization of Bovine Longissimus Dorsi Muscle during Prolonged Ageing. *Food Chem.* **2012**, *135*, 2052–2069. [CrossRef]
70. Malheiros, J.M.; Braga, C.P.; Grove, R.A.; Ribeiro, F.A.; Calkins, C.R.; Adamec, J.; Chardulo, L.A.L. Influence of Oxidative Damage to Proteins on Meat Tenderness Using a Proteomics Approach. *Meat Sci.* **2019**, *148*, 64–71. [CrossRef]

71. Gagaoua, M.; Claudia Terlouw, E.M.; Boudjellal, A.; Picard, B. Coherent Correlation Networks among Protein Biomarkers of Beef Tenderness: What They Reveal. *J. Proteom.* **2015**, *128*, 365–374. [CrossRef] [PubMed]

72. Cai, C.; Masumiya, H.; Weisleder, N.; Matsuda, N.; Nishi, M.; Hwang, M.; Ko, J.-K.; Lin, P.; Thornton, A.; Zhao, X.; et al. MG53 Nucleates Assembly of Cell Membrane Repair Machinery. *Nat. Cell Biol.* **2009**, *11*, 56–64. [CrossRef] [PubMed]
<figure><img src="image_27.png" /></figure>

Article

# Protein Array-Based Approach to Evaluate Biomarkers of Beef Tenderness and Marbling in Cows: Understanding of the Underlying Mechanisms and Prediction

Mohammed Gagaoua *,†, Muriel Bonnet and Brigitte Picard *

National Research Institute for Agriculture, Food and Environment, Université Clermont Auvergne, VetAgro Sup, UMR Herbivores, F-63122 Saint Genès Champanelle, France; muriel.bonnet@inrae.fr

* Correspondence: gmberv2001@yahoo.fr (M.G.); brigitte.picard@inrae.fr (B.P.)

† Present Address: Food Quality and Sensory Science Department, Teagasc Ashtown Food Research Centre, Ashtown, Dublin 15, Ireland.

Received: 1 August 2020; Accepted: 21 August 2020; Published: 26 August 2020

**Abstract:** This study evaluated the potential of a panel of 20 protein biomarkers, quantified by Reverse Phase Protein Array (RPPA), to explain and predict two important meat quality traits, these being beef tenderness assessed by Warner–Bratzler shear force (WBSF) and the intramuscular fat (IMF) content (also termed marbling), in a large database of 188 Protected Designation of Origin (PDO) Maine-Anjou cows. Thus, the main objective was to move forward in the progression of biomarker-discovery for beef qualities by evaluating, at the same time for the two quality traits, a list of candidate proteins so far identified by proteomics and belonging to five interconnected biological pathways: (i) energy metabolic enzymes, (ii) heat shock proteins (HSPs), (iii) oxidative stress, (iv) structural proteins and (v) cell death and protein binding. Therefore, three statistical approaches were applied, these being Pearson correlations, unsupervised learning for the clustering of WBSF and IMF into quality classes, and Partial Least Squares regressions (PLS-R) to relate the phenotypes with the 20 biomarkers. Irrespective of the statistical method and quality trait, seven biomarkers were related with both WBSF and IMF, including three small HSPs (CRYAB, HSP20 and HSP27), two metabolic enzymes from the oxidative pathway (MDH1: Malate dehydrogenase and ALDH1A1: Retinal dehydrogenase 1), the structural protein MYH1 (Myosin heavy chain-IIx) and the multifunctional protein FHL1 (four and a half LIM domains 1). Further, three more proteins were retained for tenderness whatever the statistical method, among which two were structural proteins (MYL1: Myosin light chain 1/3 and TNNT1: Troponin T, slow skeletal muscle) and one was glycolytic enzyme (ENO3: β-enolase 3). For IMF, two proteins were, in this trial, specific for marbling whatever the statistical method: TRIM72 (Tripartite motif protein 72, negative) and PRDX6 (Peroxiredoxin 6, positive). From the 20 proteins, this trial allowed us to qualify 10 and 9 proteins respectively as strongly related with beef tenderness and marbling in PDO Maine-Anjou cows.

**Keywords:** meat tenderness; fat; proteomics; proteins; enzymes; quantification; biological pathways; chemometrics

## 1. Introduction

During the last two decades, OMICs techniques, especially proteomics, have been applied by meat scientists to understand the modifications occurring in post-mortem muscle in an attempt to explain the variation in several meat quality traits [1–3]. Further, proteomics allowed the identification of putative biomarkers (for review: [4–6]), with the objective of predicting the potential quality and also of proposing a molecular test for the beef industry [7]. Overall, two-dimensional gel electrophoresis
combined with mass spectrometry was efficiently used to map and characterize bovine muscle proteins [8]. The main protein biomarkers so far identified belong to myriad interconnected pathways, such as structure and contraction, heat stock proteins, energy metabolism including the glycolytic and oxidative pathways, oxidative stress, transport, binding and signaling, apoptosis, and proteolysis, including endogenous muscle inhibitors such as serpins, cystatins and calpastatins [4,5,9,10].

Among the most investigated beef qualities using proteomics, tenderness has gained the most attention [11,12]. Beef tenderness is considered worldwide to be one of the most critical quality attributes for consumers and for re-purchase decisions. However, with beef tenderness being a multifactorial trait, it is highly variable, and is impacted by both intrinsic and extrinsic factors measurable along the continuum from farm-to-fork [13,14]. At the carcass and muscle levels, the intramuscular fat (IMF) content, also termed marbling, impacts tenderness variation and is considered as an important driver of beef palatability [15,16]. Proteomics was also applied to identify potential biomarkers of IMF and to address differences in adiposity [17] under several factors, such as breed [18], rearing practices [19], individual variability [20–22] and distinct stages of IMF development [23,24].

Based on the work by Rifai et al. [25], we recently detailed the process of meat quality biomarker discovery that should be followed so as to identify, evaluate and validate protein biomarkers at the research/industry scale [4,26]. This process is composed of six main steps that are discovery/identification, qualification, verification, research assay optimization, industrial validation and commercialization. In fact, the products of the discovery phase are lists of 10 to 100 proteins with different abundances between two compared situations or conditions (i.e., tender versus tough, or lean versus fat, etc.). Thanks to the current developments in the high-throughput techniques, several research groups have been able to move forward in the process of beef tenderness and IMF (marbling) biomarkers discovery by assaying fast techniques for the qualification of biomarkers in a few animals, which entails the confirmation of the differential abundances of the proteins using a method different from the one used for the discovery step. Thus, some shortlisted biomarkers were tested using immune-based techniques such as western blotting [27–29], Dot-Blot [30–34] and Reverse Phase Protein Array (RPPA) [35,36], or using label-free gel mass spectrometry tools, namely Selected Reaction Monitoring (SRM) and the sequential window acquisition of all theoretical spectra (SWATH) [37]. More recently, a combination of RPPA and Parallel Reaction Monitoring (PRM) [26] was used to both qualify and verify the reliability of a set of 10 proteins for predicting tenderness and marbling. Therefore, this study was designed to check the ability of a previous list of 20 protein candidate markers quantified by the RPPA technique to discriminate both beef tenderness and IMF in a large database of PDO Maine-Anjou cows, comprising 188 animals. It further aimed to qualify the most robust candidates that are predictors of both tenderness and marbling whatever the statistical method. The data of this trial will further increase our understanding of the role played by the qualified proteins in these two important beef qualities, to propose in the future generic biomarker-based tools for the early sorting of carcasses to meet consumer and industrial expectations.

## 2. Materials and Methods

### 2.1. *Experimental Designs, Cows Handling and Slaughtering*

This trial was conducted using two replicate groups of 110 and 78 animals for a total of 188 Protected Designation of Origin (PDO) Maine-Anjou cows representative of the "Rouge des Prés" breed (for details on the experimental designs refer to Gagaoua, et al. [38] and Picard et al. [19]). The Rouge des Prés breed has, since 2004, been approved to be used in France for PDO meat production, hence it has a special economic importance for the valorization of local breeds as it is a dual-purpose cattle used for both beef and milk. PDO is the name of a geographical region or specific area that is recognized by official rules to produce certain foods with special characteristics related to location. The PDO regulation covers agricultural products and foodstuffs that are produced, processed and prepared in a given geographical area using recognized know-how in this specific zone. The certification label
PDO and other labels are required by the European Commission for assuring the authenticity of food products.

The main characteristics of the PDO Maine-Anjou breed are an age at slaughter lower than 10 years, having calved at least once and a minimal carcass weight of 380 kg. In this trial, all the cows had these qualities, with average age at slaughter and carcass weight of 67.4 ± 13.9 months and 445 ± 676 kg, respectively. The cows originated from the north-western part of France and were collected from the same cooperative of livestock farmers located in the department of Maine-et-Loire. All the animals were slaughtered in industrial slaughterhouses (n= 110 at Elivia, Lion d'Angers, France and n = 78 at Charal, Sablé sur Sarthes, France) following the same protocol. The cows had free access to water before their slaughter but food was deprived for 24 h. The exsanguination from the jugular vein was performed after electrical stunning using a captive-bolt pistol. Slaughtering was performed in compliance with French welfare and by respecting EU regulations (Council Regulation (EC) No. 1099/2009). The carcasses were dressed according to standard commercial practices and between 30-50 min post exsanguination the carcasses were split in half then chilled for 24 h at 2-4 °C. None of the carcasses were electrically stimulated. Ultimate pH was recorded at 24 h post-mortem for all the carcasses using a pH meter equipped with a glass electrode, and none of them had pHu > 6.0 (the benchmark used to sort DFD carcasses).

## 2.2. Muscle and Meat Steaks Sampling

*Longissimus thoracis* (LT) muscle samples, known as a mixed fast oxido-glycolytic muscle, were taken from the 5th rib of the left-hand side of each carcass 24 h post-mortem. This sampling time is in line with the outcome of the discovery and validation of protein biomarkers of bee quality traits aiming at a management/prediction of the potential quality of the carcasses early post-mortem. Therefore, in this study and for three reasons we used muscle samples taken at 24 h post-mortem. First, this is in line with previous proteomics investigations from our team for the PDO Maine-Anjou breed and other laboratories for other animal types. Second, the sampling at industrial level was only possible at this time due to the company facilities and technical considerations. Third, at this time, we expect that the dynamic properties of the muscle expressed by changes in few but important proteins belonging to the energy metabolic or heat shock proteins pathways are less variable. Thus, the first part, free of connective tissue, was subsequently frozen in liquid nitrogen and stored at -80 °C until protein extractions for the quantification of the protein biomarkers using the RPPA technique. The second part of the sample was cut into 1-2 cm cross-section pieces, vacuum packed, and stored at -20 °C for intramuscular fat (IMF) content determination. The third part was cut into 20 mm thick steaks and placed in sealed plastic bags under vacuum and aged for 14 days at 4 °C (the usual and standard condition of ageing). Each loin sample was then frozen and stored at -20 °C awaiting Warner-Bratzler shear force (WBSF) measurements.

## 2.3. Intramuscular Fat Content Determination

The amount of IMF on each muscle sample was determined using a Dionex ASE 200 Accelerated Solvent Extractor (Dionex Corporation, Sunnyvale, CA, USA) as previously described [38]. Briefly, 1 ± 0.001 g of meat powder was placed in a 22 mL extraction cell initially prepared with a cellulose filter and silicon balls. Then, petroleum ether at a temperature and pressure of 125 °C and 103 bars respectively was used for the extraction. The slurry containing both fat and petroleum ether was collected and transferred into a previously weighed evaporation vial (±0.001 g). After 15 min of evaporation, the vial was placed in a drying oven at 105 °C for 17 h and then weighed (±0.001 g) to determine the amount of IMF in the sample. The results were expressed as the % of IMF in fresh meat.

## 2.4. Meat Tenderness Measurement by Warner-Bratzler Shear Force (WBSF)

For objective beef tenderness determination, Warner-Bratzler shear force (WBSF), known as a negative proxy of sensory tenderness, was measured according to Lepetit and Culoli [39] using an
INSTRON 5944 testing machine. Briefly, the frozen steaks firstly thawed for 48 h at 4 °C, then were placed for 4 h in a thermostated bath at 18 °C before cooking on a double grooved plate griddle (SOFRACA, Morangis, France) set at 300 °C until the end-point temperature of 55 °C, which is the usual cooking temperature in France [40]. A temperature probe (Type K-Thermocouple, HI 98704, HANNA Instruments, Newark, NJ, USA) at the geometric center of the steak was used to control the end-point cooking temperature. After cooking, each steak sample was used to prepare five cores (1 cm × 1 cm × 4 cm) parallel to the longitudinal orientation of the muscle fiber. WBSF was assessed 2 or 3 times per core in order to obtain around 10 repetitions per sample. A 1 kN load cell and a 60 mm/min crosshead speed were used (universal testing machine, MTS, Synergie 200H). The force at the rupture during shear compression testing was expressed in N/cm².

## 2.5. Protein Extraction and Quantification

The muscle proteins were extracted by homogenization of the samples in "Precellys 24" tissue homogenizer (Bertin technologies, Saint Quentin-en-Yvelines, France) following the previously described protocol [41]. Briefly, 80 mg of frozen muscle were ground using 1.4 mm ceramic beads in an extraction buffer containing 50 mM Tris (pH 6.8), 5% glycerol, 2% SDS, 2 mM DTT, 2.5 mM EDTA, 2.5 mM EGTA, 1 × HALT Phosphatase inhibitor, Protease inhibitor cocktail complete MINI EDTA-free, 2 mM Na₃VO₄ and 10 mM NaF. The extracts were then boiled for 10 min at 100 °C, sonicated to reduce viscosity and centrifuged for 10 min at 25,000× g. The supernatants were collected and stored at -80 °C until use for protein assay and biomarkers quantification by RPPA.

Protein concentrations were determined with a commercial protein assay (Pierce BCA reducing agent compatible kit, ref. 23252, Thermo Scientific, Waltham, MA, USA) with bovine serum albumin (BSA) as standard.

## 2.6. Reverse Phase Protein Array (RPPA) for Protein Biomarkers Quantification

The relative abundances of 20 protein biomarkers of tenderness and/or IMF (Table 1) were quantified by the Reverse Phase Protein Array [42] following exactly the same protocol recently detailed by our group on bovine muscle [26,35,36,41]. The 20 proteins belong to five biological pathways that are:

(i) *Energy metabolic enzymes (n = 7)*: Malate dehydrogenase (MDH1), β-enolase 3 (ENO3), Retinal dehydrogenase 1 (ALDH1A1), Triosephosphate isomerase (TPI1), Phosphoglycerate kinase 1 (PGK1), Fructose-bisphosphate aldolase (ALDOA) and Glycogen phosphorylase (PYGB);

(ii) *Heat shock proteins (n = 5)*: αB-crystallin (CRYAB), Hsp20 (HSPB6), Hsp27 (HSPB1), Hsp40 (DNAJA1) and Hsp70-1A (HSPA1A);

(iii) *Oxidative stress proteins (n = 1)*: Peroxiredoxin6 (PRDX6);

(iv) *Structural proteins (n = 5)*: MLC-1F (MYL1), Myosin heavy chain-IIx (MYH1), Troponin T, slow skeletal muscle (TNNT1), Titin (TTN) and Tubulin alpha-4A chain (TUBA4A);

(v) *Cell death and protein binding (n = 2)*: Tripartite motif protein 72 (TRIM72) and Four and a half LIM domains 1 (FHL1).

After protein quantification by RPPA and for the determination of the relative abundance of each protein, the raw data were normalized using NormaCurve following the method described by Troncale et al. [43]. This is a SuperCurve-based method that simultaneously quantifies and normalizes RPPA data for fluorescent background per spot, the total protein stain and the potential spatial bias on the slide. Then, each RPPA slide was median-centered and scaled (divided by median absolute deviation). Further corrections to the sample loadings effects were performed individually for each array by correcting the dependency of the data for individual arrays with the median value of each sample over all 20 arrays using a linear regression.
<table><caption>**Table 1.** List of the 20 protein biomarkers quantified using the Reverse Phase Protein Array (RPPA) technique <fnref n="1" />.</caption><thead><tr><th>Protein Biomarkers Name (Gene)</th><th>Uniprot ID</th><th>Monoclonal (Mo) or Polyclonal (Po) Antibodies References</th><th>Antibody Dilutions</th></tr></thead><tbody><tr><td colspan="4"><strong>Energy metabolic enzymes</strong></td></tr><tr><td>Malate dehydrogenase (MDH1)</td><td>Q3T145</td><td>Mo. anti-pig Rockland 100-601-145</td><td>1/1000</td></tr><tr><td>β-enolase 3 (ENO3)</td><td>Q3ZC09</td><td>Mo. anti-human Abnova Eno3 (M01), clone 5D1</td><td>1/30,000</td></tr><tr><td>Retinal dehydrogenase 1 (ALDH1A1)</td><td>P48644</td><td>Po. anti-bovine Abcam ab23375</td><td>1/500</td></tr><tr><td>Triosephosphate isomerase (TPH1)</td><td>Q5E956</td><td>Po. anti-human Novus NBP1-31470</td><td>1/50,000</td></tr><tr><td>Phosphoglycerate kinase 1 (PGK1)</td><td>Q3T0P6</td><td>Po. anti-human Abcam ab90787</td><td>1/5000</td></tr><tr><td>Fructose-bisphosphate aldolase (ALDOA)</td><td>A6QLL8</td><td>Po. anti-human Sigma AV48130</td><td>1/4000</td></tr><tr><td>Glycogen phosphorylase (PYGB)</td><td>Q3B7M9</td><td>Po. anti-human Santa Cruz SC-46347</td><td>1/250</td></tr><tr><td colspan="4"><strong>Heat shock proteins</strong></td></tr><tr><td>Alpha-crystallin B chain (CRYAB)</td><td>P02510</td><td>Mo. anti-bovine Assay Designs SPA-222</td><td>1/1000</td></tr><tr><td>Heat shock protein beta-6, Hsp20 (HSPB6)</td><td>Q148F8</td><td>Mo. anti-human Santa Cruz HSP20-11:SC51955</td><td>1/500</td></tr><tr><td>Heat shock protein beta-1, Hsp27 (HSPB1)</td><td>Q3T149</td><td>Mo. anti-human Santa Cruz HSP27 (F-4):SC13132</td><td>1/3000</td></tr><tr><td>DnaJ homolog subfamily A member 1, Hsp40 (DNAJA1)</td><td>Q5E954</td><td>Mo. anti-human Santa Cruz HSP40-4 (SPM251):SC-56400</td><td>1/250</td></tr><tr><td>Heat shock 70 kDa protein 1A, Hsp70-1A (HSPA1A)</td><td>Q27975</td><td>Mo. anti-human RD Systems MAB1663</td><td>1/1000</td></tr><tr><td colspan="4"><strong>Oxidative stress proteins</strong></td></tr><tr><td>Peroxiredoxin-6 (PRDX6)</td><td>O77834</td><td>Mo. anti-human Abnova PRDX6 (M01), clone 3A10-2A11</td><td>1/500</td></tr><tr><td colspan="4"><strong>Structural proteins</strong></td></tr><tr><td>Myosin light chain 1/3 (MYL1)</td><td>A0JNJ5</td><td>Po. anti-human Abnova MYL1 (A01)</td><td>1/1000</td></tr><tr><td>Myosin heavy chain-IIx (MYH1)</td><td>Q9BE40</td><td>Mo. anti-bovine Biocytex 8F4</td><td>1/500</td></tr><tr><td>Troponin T, slow skeletal muscle (TNNT1)</td><td>Q8MKH6</td><td>Po. anti-human Sigma SAB2102501</td><td>1/4000</td></tr><tr><td>Titin (TTN)</td><td>Q8WZ42</td><td>Mo. anti-human Novocastra NCL-TITIN</td><td>1/100</td></tr><tr><td>Tubulin alpha-4A chain (TUBA4A)</td><td>P81948</td><td>Mo. anti-human Sigma T6074</td><td>1/1000</td></tr><tr><td colspan="4"><strong>Cell death and protein binding</strong></td></tr><tr><td>Tripartite motif protein 72 (TRIM72)</td><td>E1BE77</td><td>Po. anti-human Sigma SAB2102571</td><td>1/2000</td></tr><tr><td>Four and a half LIM domains 1 (FHL1)</td><td>F1MR86</td><td>Po. anti-human Sigma AV34378</td><td>1/5000</td></tr></tbody></table>

## 2.7. Statistical Analyses

The statistical analyses were carried out using XLSTAT 2018.3 (AddinSoft, Paris, France). Raw data were scrutinized for data entry errors and outliers using Smirnov–Grubb's outlier test at a significance level of 5%. Then, all the data were normalized for replicates (experiment) and the factor related to the rearing practices of the animals [19,44]. This step was based on Z-scores, which represent the number of standard deviations for each observation relative to the mean of the corresponding replicate/condition. Therefore, after this transformation, the data had a mean of 0 and standard deviation of 1. Following this first step, three main statistical approaches were applied to the whole database to predict/explain each quality trait and evaluate the potential of each biomarker to contribute to its associated variation.

**Correlations:** Pearson correlation coefficients, based on the Z-scores, at the level of 5% were computed between WBSF values and IMF content with the 20 protein biomarkers.

**Clustering into WBSF and IMF classes:** three unsupervised learning methods, which were (i) hierarchical cluster analysis (HCA), (ii) *k*-means and (iii) partitioning around medoids (PAM), were tested as previously described [45] to create meat quality classes of WBSF and IMF. For both quality traits, *k*-means gave the best results based on the average silhouette width (Si) criterion (Euclidean distance), allowing in the two cases three clusters or classes (*k* = 3) that we named tender, medium and tough for WBSF, and for IMF (marbling) fat, medium and lean.

<footnote n="1">¹ The suppliers and conditions for each primary antibody after western blotting validation are given as in Picard et al. [36].</footnote>
The value of *Si* ranges from -1 to +1, with observations that have a positive large *Si* being very well clustered. Those close to 0 (low *Si*) means that the observation lies between two clusters and those with negative *Si* are partitioned in the wrong cluster [46]. Afterward, variance analyses (ANOVA) were used to compare the protein abundances among the classes for each beef quality trait. Significant differences were performed using Tukey's test at a significance level of 5% and were presented using heatmaps. Subsequently, principal component analyses (PCA) for each beef quality trait were carried out using the significant differential proteins (*p* < 0.05) to illustrate in a more complete picture the separation of the classes, and thus of the individuals and the distribution of the variables.

*Partial Least Squares regressions (PLS-R)*: to deepen our understanding of the mechanisms and identify the most robust explanatory protein biomarkers, Partial Least Squares regressions (PLS-R) were performed per beef quality trait to generate explanatory models using the optimal number of components in each case [47]. This is an appropriate tool to include all the 20 biomarkers in one model and identify those that had a biological and relevant significance using the criterion of the variable's importance in the projection (VIP). This filter method based on VIP scores estimates the importance of each protein in the projection used in a PLS model. A protein with a VIP > 1.0 is considered important, thus highly influential in a given model, 0.8 < VIP < 1.0 is considered moderately influential, and any protein with VIP < 0.8 is less influential thus considered weak and rejected. For the selection of the variables, the jack-knife method was included in the PLS-R as a selective parameter. In this study, all proteins for which the VIP scores were above a threshold of 1.0 and 0.8 (highly and moderately influential proteins) were considered and then compared with those selected from correlation and variance analyses to be used in the future for validation on meat from PDO Maine-Anjou cows.

**3. Results**

3.1. *Pearson Correlation Analyses between the Biomarkers and Meat Quality Traits*

The correlation coefficients computed between WBSF values and IMF content with the 20 protein biomarkers and for all the 188 cows are given in Table 2.

For WBSF, 11 proteins were significantly (*p* < 0.05) correlated with this texture trait. Four proteins were from the energy metabolism pathway: MDH1, ENO3, PGK1 (all negative) and ALDH1A1 (positive). These were followed by three small heat shock proteins (sHSP): CRYAB, HSP20 and HSP27 (all positive). Three were structural proteins: MYL1 and MYH1 (both negative) and TNNT1 (positive). One protein from the last pathway, which was FHL1, was positively correlated with WBSF.

For IMF content, 10 proteins were significantly (*p* < 0.05) correlated, of which 2 were energy metabolic enzymes that are MDH1 (negative) and ALDH1A1 (positive), 3 were sHSP (CRYAB, HSP20 and HSP27) that were all positive, and 2 were structural proteins (MYH1 (negative) and TNNT1 (positive)). The PRDX6 from the oxidative pathway was positive, and two proteins from the last pathway, TRIM72 and FHL1, were respectively negatively and positively correlated with IMF content (Table 2).

From the above, eight proteins, MDH1, ALDH1A1, CRYAB, HSP20, HSP27, MYH1, TNNT1 and FHL1, were common to the two traits. Interestingly, sHSPs seemed from this analysis to be important biomarkers for both WBSF and IMF.
<table><caption>**Table 2.** Pearson correlation analyses between the 20 protein biomarkers with WBSF values and Intramuscular fat (IMF).</caption><thead><tr><th>Protein Biomarkers <sup>1</sup></th><th>WBSF</th><th>IMF</th></tr></thead><tbody><tr><td colspan="3"><strong>Energy metabolic enzymes</strong></td></tr><tr><td>MDH1</td><td>-0.29 ***</td><td>-0.18 *</td></tr><tr><td>ENO3</td><td>-0.21 **</td><td>-</td></tr><tr><td>ALDH1A1</td><td>+0.34 ***</td><td>+0.38 ***</td></tr><tr><td>TP11</td><td>-</td><td>-</td></tr><tr><td>PGK1</td><td>-0.15 *</td><td>-</td></tr><tr><td>ALDOA</td><td>-</td><td>-</td></tr><tr><td>PYGB</td><td>-</td><td>-</td></tr><tr><td colspan="3"><strong>Heat shock proteins</strong></td></tr><tr><td>CRYAB</td><td>+0.32 ***</td><td>+0.37 ***</td></tr><tr><td>HSP20 (HSPB6)</td><td>+0.21 **</td><td>+0.31 ***</td></tr><tr><td>HSP27 (HSPB1)</td><td>+0.28 ***</td><td>+0.22 **</td></tr><tr><td>HSP40 (DNAJA1)</td><td>-</td><td>-</td></tr><tr><td>HSP70-1A</td><td>-</td><td>-</td></tr><tr><td colspan="3"><strong>Oxidative stress proteins</strong></td></tr><tr><td>PRDX6</td><td>-</td><td>+0.20 *</td></tr><tr><td colspan="3"><strong>Structural proteins</strong></td></tr><tr><td>MYL1</td><td>-0.18 *</td><td>-</td></tr><tr><td>MYH1</td><td>-0.26 ***</td><td>-0.26 **</td></tr><tr><td>TNNT1</td><td>+0.22 **</td><td>+0.15 *</td></tr><tr><td>TTN</td><td>-</td><td>-</td></tr><tr><td>TUBA4A</td><td>-</td><td>-</td></tr><tr><td colspan="3"><strong>Cell death and protein binding</strong></td></tr><tr><td>TRIM72</td><td>-</td><td>-0.33 ***</td></tr><tr><td>FHL1</td><td>+0.18 *</td><td>+0.22 **</td></tr></tbody></table>

## 3.2. Discriminant Biomarkers of WBSF and Marbling

In the attempt to evaluate discriminant biomarkers for both WBSF and IMF classes, the *k*-means algorithm, as the best clustering method, allowed the identification of three classes for each beef quality trait.

For WBSF, the *k*-means clustering of the 188 steaks of the PDO Maine-Anjou cows categorized them into tender (*n* = 93), medium (*n* = 71) and tough (*n* = 24) samples (Figure 1). The tender class has a mean value of 32.96 ± 3.99 N/cm², a coefficient of variation (CV) of 12% and WBSF values ranging between 23.05 and 38.76 N/cm². The medium tenderness class has a mean value of 44.74 ± 3.69 N/cm², a CV of 8% and WBSF values ranging between 39.00 and 52.22 N/cm². The tough class has a mean value of 61.18 ± 7.87 N/cm², a CV of 13% and WBSF values ranging between 53.03 and 81.49 N/cm². Comparison of the protein abundances of the 20 biomarkers, based on variance analysis among the three clusters, highlighted that 11 proteins were significantly different at the level of 5% (Figure 1a). Among them, three proteins were from the energy metabolism pathway, MDH1 and ENO3 being highly abundant in the tender class compared to the tough samples, and the inverse was found for ALDH1A1 (*p* < 0.001). The three sHSPs (CRYAB, HSP20 and HSP27) were all highly abundant in the tough compared to the tender class. Among the structural proteins, three proteins were discriminant, which were MYL1 and MYH1 (highly abundant in tender meat), and TNNT1 was found to be more abundant in the tough class. Finally, TRIM72 and FHL1 from the last pathway were respectively highly and less abundant in the tender class, compared to the tough class (Figure 1a). The main significant discriminant proteins were MDH1, ALDH1A1, CRYAB and HSP27 (*p* < 0.001), followed by MYH1, TNNT1 and FHL1 (*p* < 0.01), and then ENO3, HSP20, MYL1 and TRIM72 (*p* < 0.05). The projection of these WBSF discriminant proteins on a PCA allowed acceptable separation of the three WBSF classes (Figure 1b). The first two principal components (PC) explained around 50% of the WBSF variability, with most variation being explained by the first PC (34.1%). The cows characterized by tender meat,

<footnote n="1">¹ The correlation coefficients (q = 8) in bold font highlight the protein biomarkers common to WBSF and IMF. WBSF: Warner–Bratzler shear force; IMF: Intramuscular fat. Significance: * p < 0.05; ** p < 0.01; *** p < 0.001.</footnote>
thus being in the tender class, were all loaded on the left, the medium were in the center, and the tough were on the right. A total of six proteins characterize the tough class and the remaining five proteins were higher in the tender class (Figure 1b).

<figure><img src="image_28.png" /><figcaption>**Figure 1.** Protein biomarkers differing among the tenderness classes (Tender (n = 93), Medium (n = 71) and Tough (n = 24)). (a) Heatmap comparing the protein abundances among the three WBSF (tenderness) classes. Significance—ns: not significant; \*: p < 0.05; \*: p < 0.01; \*:*: p < 0.001. The proteins are given by their biological family following the legend. (b) Principal component analysis highlighting the distribution of the individuals of each tenderness class based on the 11 discriminant protein biomarkers. Individuals belonging to the same class are encircled in clusters using the corresponding schematic colors. The descriptive statistics of the three tenderness classes are as follows—**Tender class**: mean value of 32.96 ± 3.99 N/cm² (CV, 12%), Min = 23.05 and Max = 38.76 N/cm². **Medium class**: mean value of 44.74 ± 3.69 N/cm² (CV, 8%), Min = 39.00 and Max = 52.22 N/cm². **Tough class**: mean value of 61.18 ± 7.87 N/cm² (CV, 13%), Min = 53.03 and Max = 81.49 N/cm².</figcaption></figure>

For IMF, the 188 steaks of the PDO Maine-Anjou cows were categorized into three marbling classes, as they were for WBSF, (Figure 2) namely fat (n = 28), medium (n = 69) and lean (n = 87). The fat class has a mean value of 7.72 ± 1.58%, a CV of 20% and values ranging between 6.34 and
13.82%. The medium fat class has a mean value of 4.72 ± 0.63%, a CV of 13% and IMF values ranging between 3.76 and 6.11%. The lean class has a mean value of 2.72 ± 0.62%, a CV of 23% and IMF values ranging between 0.45 and 3.69%. Eleven proteins were significantly different (*p* < 0.05) among the three marbling classes (Figure 2a). The main discriminant biomarkers were from the HSP superfamily, with a total of four proteins (CRYAB, HSP20, HSP27 and HSP40) that were all highly abundant in the fat class compared to the others (Figure 2a). This is followed by the energy metabolism pathway, with three proteins: MDH1 (high in lean class), ALDH1A1 (high in fat class) and PYGB (high in medium class). PRDX6 from the oxidative stress pathway and MYH1 from the structural pathway were respectively high and low in the fat class (Figure 2a,b). From the last family group, the two proteins TRIM72 and FHL1 were both different in their abundance among the marbling classes, being respectively low and high in the fat class. FHL1 was in this trial more abundant in the intermediate fat group. The main most significant discriminant proteins were ALDH1A1, CRYAB and TRIM72 (*p* < 0.001), followed by HSP20, MYH1 and FHL1 (*p* < 0.01), and then MDH1, PYGB, HSP27, HSP40 and PRDX6 (*p* < 0.05). The projection of these 11 discriminant proteins on a PCA allowed for separating efficiently the three marbling classes, especially the lean and fat samples (Figure 2b). The first two PC explained around 44% of the variation, with 30.5% in the first PC.

From the above, the clustering analysis allowed us to observe that eight proteins (MDH1, ALDH1A1, CRYAB, HSP20, HSP27, MYH1, TRIM72 and FHL1) are common for the two traits, and clearly delineate the beef quality classes. As for the correlation analyses, sHSPs seemed also to be important biomarkers for both WBSF and IMF.

### 3.3. Partial Least Squares for the Prediction of WBSF and Marbling Using the Panel of 20 Protein Biomarkers

The investigation using PLS-R of the relationships between the 20 protein biomarkers and the two beef quality traits evaluated on the 188 PDO Maine-Anjou cows generated explanatory models with the main drivers of their variation (Figure 3a,b). The WBSF model retained 10 proteins (Figure 3a). Among them, eight proteins had VIP values > 1.0 (ALDH1A1, CRYAB, MDH1, HSP27, MYH1, TNNT1, ENO3 and HSP20) and 2 had 0.8 < VIP < 1.0 (FHL1 and MYL1). Three proteins for the energy metabolism, small HSP and structural proteins pathways, respectively, were retained. From the cell death and binding protein family, one protein, FHL1, was retained for WBSF.

For IMF, the PLS-R generated an explanatory model with nine proteins (Figure 3a). Among the retained proteins, seven had a VIP value > 1.0 (ALDH1A1, CRYAB, TRIM72, HSP20, MYH1, HSP27 and FHL1) and two had 0.8 < VIP < 1.0 (PRDX6 and MDH1).

For both PLS-R explanatory models, seven proteins were common (ALDH1A1, CRYAB, MDH1, HSP20, HSP27, MYH1 and FHL1) in explaining WBSF and IMF variations. Interestingly, for both quality traits models, the first ranked proteins with VIP values > 1.5 were ALDH1A1 and CRYAB (Figure 3a,b).
<figure><img src="image_29.png" /><figcaption>**Figure 2.** Protein biomarkers differing among the marbling classes (Fat (n = 28), Medium (n = 69) and Lean (n = 87)). (a) Heatmap comparing the protein abundances among the three IMF (marbling) classes. Significance—ns: not significant; \*: *p* < 0.05; \*:*: *p* < 0.01; \*:\*: *p* < 0.001. The proteins are given by their biological family following the legend. (b) Principal component analysis highlighting the distribution of the individuals of each marbling class based on the 11 discriminant protein biomarkers. Individuals belonging to the same class are encircled in clusters using the corresponding schematic colors. The descriptive statistics of the three marbling classes are as follows—**Fat class:** mean value of 7.72 ± 1.58% (CV, 20%), Min = 6.34 and Max = 13.82%. **Medium class:** 4.72 ± 0.63% (CV, 13%), Min = 3.76 and Max = 6.11%. **Lean class:** 2.72 ± 0.62% (CV, 23%), Min = 0.45 and Max = 3.69%.</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>**Figure 3.** Partial least squares highlighting the protein biomarkers retained to explain (a) tenderness evaluated by WBSF and (b) IMF content (marbling) based on their variable importance in the projection (VIP). The proteins retained in positive and negative directions are shown in blue and red colors, respectively. For (a) WBSF, a total of 10 proteins were retained from which 8 had VIP > 1.0, and for (b) IMF, 9 proteins were retained, from which 7 had VIP > 1.0. A total of 7 proteins were common (MDH1, ALDH1A1, CRYAB, HSP20, HSP27, MYH1 and FHL1) in the two models to explain both WBSF and IMF variation.</figcaption></figure>

### 3.4. Summary of the Putative Common Protein Biomarkers from the Three Statistical Methods

The summary of the proteins retained from the 20 panel biomarkers, based on the results of the three statistical methods presented above to explain/predict WBSF and IMF content, is given in Figure 4. Overall, irrespective of the statistical method and quality trait and based on the variable importance of the proteins in the models, their accuracy in discriminating the quality groups and their significant associations with WBSF and IMF, seven biomarkers were common, including three small HSPs (CRYAB, HSP20 and HSP27), two energy metabolic enzymes from the oxidative pathway (MDH1 and ALDH1A1), the structural protein MYH1 and the multifunctional protein FHL1. Further, for WBSF three more proteins were retained whatever the statistical method, among which two were structural proteins (MYL1 and TNNT1) and one was a metabolic enzyme (ENO3). For IMF, two proteins were retained whatever the statistical method, these being TRIM72 and PRDX6.
<figure><img src="image_31.png" /><figcaption>**Figure 4.** Summary of the evaluation of the 20 protein biomarkers quantified by RPPA using the three statistical methods (Pearson correlations, K-means clustering and Partial Least Squares regressions (PLS-R)) to explain/predict WBSF and IMF content on the *Longissimus thoracis* muscle of the 188 PDO Maine-Anjou cows.</figcaption></figure>

## 4. Discussion

This trial aimed to evaluate the potential of 20 protein biomarkers, previously identified by proteomics to be potential markers of beef tenderness and marbling [4,11,12,17,19,22,26,34,35,48–50] and belonging to five interconnected biological pathways ((i) energy metabolic enzymes, (ii) heat shock proteins (HSPs), (iii) oxidative stress, (iv) structural proteins and (v) cell death and protein binding) to explain/predict two important beef quality traits for both the consumers and the meat sector (meat tenderness measured by WBSF and the marbling evaluated by the percentage of IMF content). The list of proteins was selected for qualification in this trial based on two main criteria: (i) association of the protein with the quality traits from previous studies, and (ii) validation of an antibody against the protein for its quantification by RPPA. Therefore, we used the immunological
RPPA technique for their semi-quantification, and applied three statistical methods to explain the variability of each quality trait, being Pearson correlations to assess the type of associations with the quality traits, unsupervised learning to perform a clustering of the quality traits and determine the main protein splitters, and Partial Least Squares regressions (PLS-R) to propose the first overall regressions models and identify the main predictor proteins based on their importance in the model. This firstly correlates the relative abundances of the candidate protein markers with WBSF and IMF, and secondly describes the consistencies and differences for the two traits.

This study is the first to use such a large database of 188 cows to perform this qualification/evaluation step on the selected 20 protein biomarkers, hence allowing us to move forward in the progression of biomarkers discovery for beef qualities [4] by evaluating at the same time two beef quality traits. From the evaluated list of proteins, we revealed (regardless of statistical method) that certain biomarkers are robust for both beef qualities of PDO Maine-Anjou cows. A robust biomarker is the protein that is identified in this study to be related with both beef quality traits whatever the statistical method. In the following sections, the proteins common to the two traits are discussed together, and those that were trait-dependent are presented separately. The biological pathways behind the associations identified between the two phenotypes and the proteins were further briefly presented. The relationships between the robust evaluated protein biomarkers and the two beef quality traits showed that the most tender meat of PDO Maine-Anjou cows had higher abundances of glycolytic enzymes, such as ENO3, and of fast contractile proteins such as MYL1 and MYH1, while they had lower abundances of slow contractile proteins such as TNNT1, and lower abundances of small HSPs with higher abundances of FHL1. These relationships are consistent with each other and are in accordance with the results of Couvreur et al. [51] based on the contractile and metabolic properties of the *Longissimus thoracis* muscle. Furthermore, the contractile properties of the *Longissimus thoracis* muscle from PDO Maine-Anjou cows have been associated with a specific fiber type composition, as the muscle contains a very low proportion of IIX fibers and a higher proportion of IIA fibers compared to the French beef breeds [38,51–53]. Consequently, its glycolytic metabolism is very low [51,52]. This demonstrates that for a slow oxidative type of muscle, the most tender are the less slow oxidative and the most fast glycolytic, as observed for Aberdeen Angus [54] or Chianina breeds [49].

## 4.1. Common Biomarkers Explaining the Variation in WBSF (Tenderness) and IMF (Marbling)

The data of the present trial showed that the common proteins between tenderness and IMF are usually inversely associated, except for FHL1 which was positive. WBSF was weakly and negatively correlated with IMF ($r = +0.16; p = 0.032$). The fatter meat samples are those containing higher proportions of proteins of the slow oxidative type, such as TNNT1 and small HSPs, and low proportions of the fast glycolytic type, such as MYH1, which is consistent with previous results about biomarkers of adiposity from the same breed [22,26]. From these data, it seems that the fatter meat samples do not lead to the tenderest beef. Thus, irrespective of the quality trait and the statistical method, seven proteins were robustly related with the WBSF and IMF content of PDO Maine-Anjou cows (Figure 4).

Without any surprise and in line with the results presented above, the superfamily of heat shock proteins, specifically the small HSP members, seemed to be just as good and important biomarkers for both beef qualities. Indeed, CRYAB, HSP20 and HSP27 were positively related in this trial with WBSF (negatively with tenderness), and so with toughness and IMF content. An inverse relationship between the abundance of small HSPs and tenderness was reported in the *Longissimus thoracis* muscle of the Angus breed [54]. It is worthwhile to note that the *Longissimus thoracis* muscles of the PDO Maine-Anjou breed investigated in this study, as well as those of Aberdeen Angus, were described as having more oxidative metabolisms with high amounts of fat content [18,38,52,54]. The slow oxidative type fibers are known to contain high levels of sHSP proteins [52], as also evidenced by Golenhofen et al. [55] in rats, showing more HSPs in slow oxidative muscles.

The involvement of these chaperones in post-mortem muscle, and thus in its conversion into meat, was reported by many earlier studies (for review: [5,11,12,56]). This is in agreement with the onset
of apoptosis, the first phase in the conversion of muscle into meat, involving major biochemical and
structural changes that influence not only the meat tenderization process but also the homeostasis of
the post-mortem muscle [9,57]. Small HSP proteins were thought to stabilize and ensure the correct
folding of newly synthesized proteins, or help refold proteins altered by cell stress to protect them
against metabolic disorders and ischemia [58]. For example, small HSPs can bind to myofibrils [59–61],
thereby protecting skeletal muscle through structural protein complexes, which partly explain the
high abundance observed in the tough meat class. Accordingly and during the process of apoptosis,
CRYAB is able to negatively regulate Cytochrome c and Caspase-8, hence inhibiting the executor
caspase 3 [62,63]. For instance, any increase in the level of CRYAB (maybe also the other small HSPs)
leads to less muscle structure degradation, which would produce tough meat [48]. In support of their
important roles in meat tenderization, small HSPs were identified by several previous proteomics and
under different conditions as good biomarkers for the two qualities investigated in this trial [6,17],
or other sensory qualities of beef such as color [5] and pH decline [52,64,65]. CRYAB, which was ranked
second in the PLS-R models, is a well-known beef tenderness biomarker [4,29,50,66–69], also recently
identified [22] using gel-based and gel-free proteomics methods for IMF of the same breed using
PRM and antibody-based proteomics [26]. HSP27 was reported by proteomics in several studies for
beef tenderness [4,27,29,35,66,67,70], and as a consistent biomarker of adiposity (marbling) in the few
proteomics that investigated this trait [18,20,22,24,71]. An absolute quantification of HSP27 in PDO
Maine-Anjou cows confirmed a higher abundance in the less tender *Longissimus thoracis* muscle [26],
as observed in this study. For HSP20, this study is the first to identify it as a biomarker for marbling,
but it is already connected by several authors to beef tenderness [4,34,35,49,66–68]. Worthy of note is
that in humans and among the HSP protein members, HSP20 was found by DeLany et al. [72] to be
the most up-regulated chaperone during the differentiation of human adipose-derived stem cells into
mature adipocytes.

In line with the above, ALDH1A1, a metabolic enzyme that is ranked first in the PLS-R models regardless of the quality trait, and is also significantly correlated with both traits and discriminates groups of tenderness and IMF, appeared as the most robust biomarker in this study with CRYAB. ALDH1A1 has been already identified as a biomarker of adiposity of PDO Maine-Anjou [22,26]. The authors proposed that any increase in this protein would mediate an increase in CRYAB and maybe in other small HSPs thanks to the retinoic acid content [73]. Further, the same authors postulated a mechanism in PDO Maine-Anjou [22,26], suggesting that a higher abundance of CRYAB resulting from a high proportion of slow fibers may sustain the elevated oxidative metabolism of marbled muscle, which may partly explain the important role that ALDH1A1 and CRYAB play in the determination of both WBSF and IMF content. ALDH1A1 was identified as a good biomarker of beef tenderness in different breeds and muscles [35,68,74,75], and was also connected to several color parameters [5], especially redness (a*) [41] and metmyoglobin-reducing activity [76]. These last parameters are associated with the oxidative metabolism resulting from the fiber composition of the *Longissimus thoracis* muscle from Rouge des Près cows (with a very low proportion of fast glycolytic fibers [52]). It is worthy of mention that ALDH1A1 is able to protect cells against the cytotoxic effects of various aldehydes [77], which are generated in the cytosol by lipid peroxidation [78] and not cleared due to the cessation of blood flow.

The second metabolic enzyme, the cytosolic MDH1, is a member of the malate dehydrogenase
enzymes very important for both gluconeogenesis and the Krebs cycle, and therefore plays a crucial
role in energy and cellular metabolism [79]. In adipocytes, including muscular adipocytes, MDH1 is
involved in the reduced nicotinamide adenine dinucleotide phosphate (NADPH) supply for de novo
fatty acid synthesis, and is considered as a lipogenic enzyme. It was retained in this study in agreement
with earlier studies. For instance, it was reported (as a biomarker of beef tenderness in bulls) to be
positive in two studies [28,54] and negative in another [4]. MDH1 was also shown as a good predictor
of *Longissimus thoracis* tenderness variability (in a few animals) when quantified by PRM [26]. Moreover,
malic enzyme activity (MDH1) and IMF quantitative values in the *Rectus abdominis* and *Semitendinosus*
muscles from Limousin, Angus and Japanese Black cross Angus steers were shown to be positively correlated [80]. Similar positive correlations between MDH1 and IMF were also observed when the abundance of MDH1 was assayed by absolute PRM quantification in both *Longissimus thoracis* and *Semitendinosus* muscles, but not by RPPA [26]. In the present study, using a larger number of animals and using the RPPA technique, an expected link between MDH1 and IMF was revealed by three mathematical methods. However, an unexpected negative correlation between MDH1 abundance and IMF probably arose from the relative and normalized abundance of MDH1 in the present study, making difficult the comparison between the current and previous results. This warrants further investigation in order to understand and describe the link between the MDH1 protein's abundance and IMF depending on the methodological specificities of both the protein quantification method and data processing.

MYH1 is the only structural protein that was related to both WBSF and IMF content whatever the statistical method (Figure 4). Changes in the cytoskeletal proteins have been shown and investigated for decades to play a role in meat tenderization [81], and are also evidenced by proteomics [4,6,11,82]. Further, myosin fibers play pivotal roles from fetal life to the slaughter of cattle (for review: [83]). MYH1 is the gene that encodes the fast glycolytic IIX fibers. In this study, the low abundance of MYH1 in the lean marbling class is in agreement with the earlier data reported for breeds of French origin, characterized by high proportions of fast glycolytic IIX fibers [84]. Cytoskeletal proteins were proposed to participate in IMF deposition [22,85,86], and MYH1 has already been identified by proteomics as a negative biomarker of adiposity [71]. MYH1 was revealed in several studies and under different factors as a good biomarker of beef tenderness [4,34,35,54,67,87,88]. As reported by Picard et al. [54] and Gagaoua et al. [34], the direction of its relationships with tenderness depends on muscle type, breed, end-point cooking temperature and origin of the panelists. In fact, in muscles with a low proportion of fast glycolytic fibers, such as the *Longissimus thoracis* muscle in breeds like Rouge des Prés or Aberdeen Angus, MYH1 was already reported to be positively related with tenderness. On the contrary, in French beef breeds in which the *Longissimus thoracis* muscle contains a high proportion of IIX fibers, the relationship with tenderness is negative. The identification of MYH1 as an important biomarker agrees with the theory stating that muscles with high proportions of fast fiber types are more susceptible to early post-mortem proteolytic degradation [89]. This can be further explained by their susceptibility to post-mortem glycolysis, hence leading to more tender than tough meat [83].

The seventh and last protein retained in this trial for the two beef quality traits was FHL1, known as a multifunctional protein regulating metabolism, cell proliferation, gene transcription and apoptosis [90]. This protein can further interact with metabolic enzymes as a response to the oxidative stress in muscle as well as hypoxia [90], thereby explaining its tendency to be projected with proteins characterizing slow oxidative properties, such as ALDH1A1 and TNNT1. To our knowledge, FHL1 was reported as a biomarker of beef tenderness in the *Longissimus thoracis* muscle in four previous studies [35,91–93], and of marbling in three studies [22,26,94]. We suggest that it is mainly via the regulation of calcium homeostasis [95] that FHL1 plays a role in beef tenderness determination. Indeed, Ca²⁺ ions contribute to the regulation of the energy metabolism pathways [96], as they affect the enzymatic speed of several crucial metabolic enzymes [11,12]. This is consistent with our results showing modifications of the abundances of the proteins ALDH1A1, MDH1 and ENO3. Calcium concentrations have also been implicated as initiators of apoptosis via some signaling pathways in skeletal muscle [12,97]. For example, apoptosis was documented to affect the integrity of the skeletal muscle, through the modification of Ca²⁺ flux during ageing and its consequences on protein proteolysis involving ultra-structural modifications. This can be supported in this trial by the robust association of TNNT1, MYL1 and MYH1, which are all proteins of structure. Further, FHL1 belongs to this last group of proteins and was described as an activator of the myostatin signaling pathway in skeletal muscle by promoting muscle atrophy [98]. This would explain the variation in the different muscles fibers and the involvement of MYH1 (further details, see [83]). For marbling, the mechanism would also be partly through the activation of myostatin signaling [98]. It is worthy of
note that in beef, FHL1 was further related negatively with the lightness (L\*) and positively with the redness (a\*) color parameters [41].

## 4.2. Biomarkers Specific to WBSF

Irrespective of the statistical method, three proteins were found to be robust and specific to WBSF (Figure 4), these being two structural proteins (MYL1 and TNNT1) and one glycolytic enzyme (ENO3). The positive association of MYL1 with tenderness (negative with WBSF) is in agreement with several earlier studies on cows [4,38], steers [68,99] and young bulls [54,66,70]. MYL1 is considered as an indicator of proteolysis [81] and is prone to phosphorylation, a reaction playing a pivotal role in muscle to meat conversion [82,100]. It is worthy to consider that this phosphorylation is induced through sarcoplasmic reticulum Ca²⁺ release in a concentration-dependent manner. TNNT1 is a slow isoform of the troponins complex that is involved in the regulation of muscle contraction [81]. The release of TNNT members, including TNNT1, has been extensively studied and they have been considered as important substrates of the endogenous muscle proteolytic systems [11,12,82]. TNNT members were thought to be easily degraded by calpains during the aging period of muscle. In fact, several proteomics trials identified TNNT1 as a good biomarker of beef tenderness [4,66,99,101]. Its negative association with tenderness in this trial agrees with what we know from the literature [102], as recently evidenced by a meta-proteomics on different muscles and genders [4], highlighting inverse relationships. In the present study, the identification of TNNT1 is consistent with the positive relation observed between the fast glycolytic isoform of MYH1 and the glycolytic enzyme ENO3, and tenderness.

ENO3 is the last protein we identified as robustly related to tenderness. Its identification as a positive biomarker of beef tenderness in this study on PDO Maine-Anjou cows and several previous proteomics on beef [4,35,67,87,92,103] points to its importance, as it is especially described as being a key moonlighting enzyme associated with hypoxic conditions and stress [104]. Enolase is a cytosolic enzyme responsible for the conversion of 2-phosphoglycerate into phosphoenolpyruvate, thereby playing an important role in pH decline and post-mortem metabolism [41]. Indeed, ENO3 induces glucose metabolism under hypoxic conditions [105], hence enacting a cellular stress response to the deprivation of oxygen supply and glucose levels.

## 4.3. Biomarkers Specific to IMF

Two proteins, Peroxiredoxin 6 (PRDX6) and Tripartite motif protein 72 (TRIM72), appeared to be robustly related as biomarkers of marbling (Figure 4). PRDX6 was already identified as a good biomarker of beef tenderness [35,49,50,68,92], color [5] and pH decline [33]. PRDX6 is a bi-functional protein with both phospholipase A2 (PLA2) and glutathione peroxidase activities, which is expressed in nearly all tissues and protects cells against oxidative stress [106]. Earlier studies reported the high abundance of PRDX6 in Aberdeen Angus (young bulls and cross sired steers), known as a marbled breed compared to young Limousin bulls [32] and Belgian Blue sired steer [18]. PRDX6 might further play a role through its PLA₂ activity via the catalysis of the hydrolysis of the acyl group at the sn-2 position of glycerophospholipids, with a specific link to phosphatidylcholine to release free fatty acids and a lysophospholipid [106], hence partly explaining its relation with IMF. Further, using PRDX6 knockout mice, Arriga et al. [107] demonstrated that PRDX6 modulates the link between glycemic and lipogenic metabolisms, thereby playing a pivotal role in fat deposition.

TRIM72, also termed MG53 (Mitsugumin 53), is a signaling protein that acts as a sensor of oxidation on membrane damage [108], including playing a crucial role in the muscle membrane repair process. Its negative relationship with IMF would be partly explained by its implication in the clearance of harmful agents collected under the cell death process and lipid oxidation. TRIM72 was reported as a negative biomarker of beef tenderness [4,35,74], and thus one can suggest that a reduced cell death phase in tough meat occurred. TRIM72 was further identified as a biomarker of beef color Lightness (L\*) [41]. The identification of TRIM72 as a biomarker of marbling in PDO Maine-Anjou would be supported by its involvement in the negative feedback regulation of myogenesis, by targeting the
insulin receptor substrate-1 [109]. Further investigations are warranted in order to clarify the exact role
of TRIM72 in the muscle to meat conversion, including the role it plays in marbling.

**5. Conclusions**

This study allowed us to qualify, on *Longissimus thoracis* muscle and using the immune-based RPPA
technique and different statistical methods, the potential of a list of 20 protein biomarkers to explain
the variation of two important beef quality traits of PDO Maine-Anjou cows. These were tenderness,
measured by the instrumental method WBSF, and the marbling of the carcasses as evaluated by IMF
content. This study is the first to propose 10 and 9 proteins as robust candidate biomarkers of WBSF
and IMF, respectively, in Rouge des Prés cows regardless the statistical method. Seven proteins are of
specific interest, as they are related to both traits. They are, in the order of importance, ALDH1A1,
CRYAB, HSP27, HSP20, MYH1, FHL1 and MDH1. These proteins belong to the superfamily of heat
shock proteins, energy metabolism (especially the oxidative pathway), and the structural proteins
(including FHL1) playing roles in cell death, metabolism and the regulation of calcium homeostasis.
Our findings further highlight that similarities exist in the biological pathways underpinning tenderness
and marbling determination.

**Author Contributions:** Conceptualization, M.G., M.B. and B.P.; software, M.G.; validation, M.G., M.B., and B.P.; formal analysis, M.G.; investigation, M.G., M.B. and B.P.; resources, M.B. and B.P.; data curation, M.G., M.B. and B.P.; writing—original draft preparation, M.G.; writing—review and editing, M.G., M.B. and B.P.; visualization, M.G.; funding acquisition, M.B. and B.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the regional council of Pays de Loire (Quafima 1 & 2).

**Acknowledgments:** The authors thank the SICA Rouge des Prés, especially A. Valais and G. Aminot for animal handling and sampling. We further thank S. Couvreur and G. Le Bec (Ecole Supérieure d'Agriculture (ESA) Angers) for muscle sampling, as well as IMF content and shear force assays. The authors thank Nicole Dunoyer for extraction of proteins and convey special thanks to the Institut Curie RPPA platform for acquisition of RPPA data under the supervision of Leanne De Koning and the support of Aurélie Cartier and Bérengère Ouine.

**Conflicts of Interest:** The authors declared no conflict of interest.

**References**

1. Jia, X.; Hollung, K.; Therkildsen, M.; Hildrum, K.I.; Bendixen, E. Proteome analysis of early post-mortem changes in two bovine muscle types: M. longissimus dorsi and M. semitendinosus. *Proteomics* **2006**, *6*, 936–944. [CrossRef]

2. Lametsch, R.; Karlsson, A.; Rosenvold, K.; Andersen, H.J.; Roepstorff, P.; Bendixen, E. Postmortem proteome changes of porcine muscle related to tenderness. *J. Agric. Food Chem.* **2003**, *51*, 6992–6997. [CrossRef] [PubMed]

3. Huang, H.; Lametsch, R. Application of Proteomics for Analysis of Protein Modifications in Postmortem Meat. In *Proteomics in Foods: Principles and Applications*; Toldrá, F., Nollet, L.M.L., Eds.; Springer: Boston, MA, USA, 2013; pp. 111–125. [CrossRef]

4. Picard, B.; Gagaoua, M. Meta-proteomics for the discovery of protein biomarkers of beef tenderness: An overview of integrated studies. *Food Res. Int.* **2020**, *127*, 108739. [CrossRef] [PubMed]

5. Gagaoua, M.; Hughes, J.; Terlouw, E.M.C.; Warner, R.D.; Purslow, P.P.; Lorenzo, J.M.; Picard, B. Proteomic biomarkers of beef colour. *Trends Food Sci. Technol.* **2020**, *101*, 234–252. [CrossRef]

6. Picard, B.; Gagaoua, M.; Hollung, K. Chapter 12—Gene and Protein Expression as a Tool to Explain/Predict Meat (and Fish) Quality. In *New Aspects of Meat Quality: From Genes to Ethics*; Purslow, P., Ed.; Woodhead Publishing: Sawston, UK; Cambridge, UK, 2017; pp. 321–354. [CrossRef]

7. Te Pas, M.F.; Hoekman, A.J.; Smits, M.A. Biomarkers as management tools for industries in the pork production chain. *J. Chain Netw. Sci.* **2011**, *11*, 155–166. [CrossRef]

8. Bouley, J.; Chambon, C.; Picard, B. Mapping of bovine skeletal muscle proteins using two-dimensional gel electrophoresis and mass spectrometry. *Proteomics* **2004**, *4*, 1811–1824. [CrossRef] [PubMed]
9. Gagaoua, M.; Hafid, K.; Boudida, Y.; Becila, S.; Ouali, A.; Picard, B.; Boudjellal, A.; Sentandreu, M.A. Caspases and Thrombin Activity Regulation by Specific Serpin Inhibitors in Bovine Skeletal Muscle. *Appl. Biochem. Biotechnol.* **2015**, *177*, 279–303. [CrossRef] [PubMed]

10. Gagaoua, M.; Dunoyer, N.; Picard, B. Differences in the levels of bov-SERPINA3 in *Rectus Abdominis* Muscle Explain the Variability of Beef Tenderness. In Proceedings of the 63rd International Congress of Meat Science and Technology, Cork, Ireland, 13–18 August 2017; pp. 832–833.

11. Picard, B.; Gagaoua, M. Chapter 11—Proteomic Investigations of Beef Tenderness. In *Proteomics in Food Science: From Farm to Fork*; Colgrave, M.L., Ed.; Academic Press: London, UK, 2017; pp. 177–197. [CrossRef]

12. Ouali, A.; Gagaoua, M.; Boudida, Y.; Becila, S.; Boudjellal, A.; Herrera-Mendez, C.H.; Sentandreu, M.A. Biomarkers of meat tenderness: Present knowledge and perspectives in regards to our current understanding of the mechanisms involved. *Meat Sci.* **2013**, *95*, 854–870. [CrossRef]

13. Gagaoua, M.; Picard, B.; Monteils, V. Assessment of cattle inter-individual cluster variability: The potential of continuum data from the farm-to-fork for ultimate beef tenderness management. *J. Sci. Food Agric.* **2019**, *99*, 4129–4141. [CrossRef]

14. Gagaoua, M.; Monteils, V.; Picard, B. Data from the farmgate-to-meat continuum including omics-based biomarkers to better understand the variability of beef tenderness: An integromics approach. *J. Agric. Food Chem.* **2018**, *66*, 13552–13563. [CrossRef]

15. Hocquette, J.F.; Gondret, F.; Baeza, E.; Medale, F.; Jurie, C.; Pethick, D.W. Intramuscular fat content in meat-producing animals: Development, genetic and nutritional control, and identification of putative markers. *Animal* **2010**, *4*, 303–319. [CrossRef] [PubMed]

16. Wood, J.D.; Enser, M.; Fisher, A.V.; Nute, G.R.; Sheard, P.R.; Richardson, R.I.; Hughes, S.I.; Whittington, F.M. Fat deposition, fatty acid composition and meat quality: A review. *Meat Sci.* **2008**, *78*, 343–358. [CrossRef] [PubMed]

17. Ceciliani, F.; Lecchi, C.; Bazile, J.; Bonnet, M. Proteomics Research in the Adipose Tissue. In *Proteomics in Domestic Animals: From Farm to Systems Biology*; de Almeida, A.M., Eckersall, D., Miller, I., Eds.; Springer: Cham, Switzerland, 2018; pp. 233–254. [CrossRef]

18. Keady, S.M.; Kenny, D.A.; Ohlendieck, K.; Doyle, S.; Keane, M.G.; Waters, S.M. Proteomic profiling of bovine M. longissimus lumborum from Crossbred Aberdeen Angus and Belgian Blue sired steers varying in genetic merit for carcass weight. *J. Anim. Sci.* **2013**, *91*, 654–665. [CrossRef]

19. Picard, B.; Gagaoua, M.; Al Jammas, M.; Bonnet, M. Beef tenderness and intramuscular fat proteomic biomarkers: Effect of gender and rearing practices. *J. Proteom.* **2019**, *200*, 1–10. [CrossRef]

20. Mao, Y.; Hopkins, D.L.; Zhang, Y.; Li, P.; Zhu, L.; Dong, P.; Liang, R.; Dai, J.; Wang, X.; Luo, X. Beef quality with different intramuscular fat content and proteomic analysis using isobaric tag for relative and absolute quantitation of differentially expressed proteins. *Meat Sci.* **2016**, *118*, 96–102. [CrossRef]

21. Shen, Y.N.; Kim, S.H.; Yoon, D.H.; Lee, H.G.; Kang, H.S.; Seo, K.S. Proteome Analysis of Bovine Longissimus dorsi Muscle Associated with the Marbling Score. *Asian Australas. J. Anim. Sci.* **2012**, *25*, 1083–1088. [CrossRef] [PubMed]

22. Bazile, J.; Picard, B.; Chambon, C.; Valais, A.; Bonnet, M. Pathways and biomarkers of marbling and carcass fat deposition in bovine revealed by a combination of gel-based and gel-free proteomic analyses. *Meat Sci.* **2019**, *156*, 146–155. [CrossRef]

23. Kim, N.K.; Lee, S.H.; Cho, Y.M.; Son, E.S.; Kim, K.Y.; Lee, C.S.; Yoon, D.; Im, S.K.; Oh, S.J.; Park, E.W. Proteome analysis of the m. longissimus dorsi between fattening stages in Hanwoo steer. *BMB Rep.* **2009**, *42*, 433–438. [CrossRef]

24. Zhang, Q.; Lee, H.G.; Han, J.A.; Kim, E.B.; Kang, S.K.; Yin, J.; Baik, M.; Shen, Y.; Kim, S.H.; Seo, K.S.; et al. Differentially expressed proteins during fat accumulation in bovine skeletal muscle. *Meat Sci.* **2010**, *86*, 814–820. [CrossRef]

25. Rifai, N.; Gillette, M.A.; Carr, S.A. Protein biomarker discovery and validation: The long and uncertain path to clinical utility. *Nat. Biotechnol.* **2006**, *24*, 971–983. [CrossRef]

26. Bonnet, M.; Soulat, J.; Bons, J.; Léger, S.; De Koning, L.; Carapito, C.; Picard, B. Quantification of biomarkers for beef meat qualities using a combination of Parallel Reaction Monitoring- and antibody-based proteomics. *Food Chem.* **2020**, *317*, 126376. [CrossRef] [PubMed]
27. Kim, N.K.; Cho, S.; Lee, S.H.; Park, H.R.; Lee, C.S.; Cho, Y.M.; Choy, Y.H.; Yoon, D.; Im, S.K.; Park, E.W. Proteins in longissimus muscle of Korean native cattle and their relationship to meat quality. *Meat Sci.* **2008**, *80*, 1068–1073. [CrossRef] [PubMed]

28. Jia, X.; Veiseth-Kent, E.; Grove, H.; Kuziora, P.; Aass, L.; Hildrum, K.I.; Hollung, K. Peroxiredoxin-6—A potential protein marker for meat tenderness in bovine longissimus thoracis muscle. *J. Anim. Sci.* **2009**, *87*, 2391–2399. [CrossRef] [PubMed]

29. Oh, E.; Lee, B.; Choi, Y.M. Associations of Heat-Shock Protein Expression with Meat Quality and Sensory Quality Characteristics in Highly Marbled Longissimus Thoracis Muscle from Hanwoo Steers Categorized by Warner–Bratzler Shear Force Value. *Foods* **2019**, *8*, 638. [CrossRef] [PubMed]

30. Guillemin, N.; Meunier, B.; Jurie, C.; Cassar-Malek, I.; Hocquette, J.F.; Leveziel, H.; Picard, B. Validation of a Dot-Blot quantitative technique for large scale analysis of beef tenderness biomarkers. *J. Physiol. Pharmacol.* **2009**, *60* (Suppl. 3), 91–97. [PubMed]

31. Guillemin, N.; Jurie, C.; Cassar-Malek, I.; Hocquette, J.F.; Renand, G.; Picard, B. Variations in the abundance of 24 protein biomarkers of beef tenderness according to muscle and animal type. *Animal* **2011**, *5*, 885–894. [CrossRef]

32. Gagaoua, M.; Terlouw, E.M.C.; Picard, B. The study of protein biomarkers to understand the biochemical processes underlying beef color development in young bulls. *Meat Sci.* **2017**, *134*, 18–27. [CrossRef]

33. Gagaoua, M.; Terlouw, E.M.; Micol, D.; Boudjellal, A.; Hocquette, J.F.; Picard, B. Understanding Early Post-Mortem Biochemical Processes Underlying Meat Color and pH Decline in the Longissimus thoracis Muscle of Young Blond d'Aquitaine Bulls Using Protein Biomarkers. *J. Agric. Food Chem.* **2015**, *63*, 6799–6809. [CrossRef]

34. Gagaoua, M.; Terlouw, C.; Richardson, I.; Hocquette, J.F.; Picard, B. The associations between proteomic biomarkers and beef tenderness depend on the end-point cooking temperature, the country origin of the panelists and breed. *Meat Sci.* **2019**, *157*, 107871. [CrossRef]

35. Gagaoua, M.; Bonnet, M.; Ellies-Oury, M.P.; De Koning, L.; Picard, B. Reverse phase protein arrays for the identification/validation of biomarkers of beef texture and their use for early classification of carcasses. *Food Chem.* **2018**, *250*, 245–252. [CrossRef]

36. Picard, B.; Gagaoua, M.; Al-Jammas, M.; De Koning, L.; Valais, A.; Bonnet, M. Beef tenderness and intramuscular fat proteomic biomarkers: Muscle type effect. *Peer J.* **2018**, *6*, e4891. [CrossRef] [PubMed]

37. Wu, W.; Dai, R.T.; Bendixen, E. Comparing SRM and SWATH Methods for Quantitation of Bovine Muscle Proteomes. *J. Agric. Food Chem.* **2019**, *67*, 1608–1618. [CrossRef] [PubMed]

38. Gagaoua, M.; Monteils, V.; Couvreur, S.; Picard, B. Identification of Biomarkers Associated with the Rearing Practices, Carcass Characteristics, and Beef Quality: An Integrative Approach. *J. Agric. Food Chem.* **2017**, *65*, 8264–8278. [CrossRef] [PubMed]

39. Lepetit, J.; Culioli, J. Mechanical properties of meat. *Meat Sci.* **1994**, *36*, 203–237. [CrossRef]

40. Gagaoua, M.; Micol, D.; Picard, B.; Terlouw, C.E.; Moloney, A.P.; Juin, H.; Meteau, K.; Scollan, N.; Richardson, I.; Hocquette, J.F. Inter-laboratory assessment by trained panelists from France and the United Kingdom of beef cooked at two different end-point temperatures. *Meat Sci.* **2016**, *122*, 90–96. [CrossRef] [PubMed]

41. Gagaoua, M.; Bonnet, M.; De Koning, L.; Picard, B. Reverse Phase Protein array for the quantification and validation of protein biomarkers of beef qualities: The case of meat color from Charolais breed. *Meat Sci.* **2018**, *145*, 308–319. [CrossRef] [PubMed]

42. Akbani, R.; Becker, K.F.; Carragher, N.; Goldstein, T.; de Koning, L.; Korf, U.; Liotta, L.; Mills, G.B.; Nishizuka, S.S.; Pawlak, M.; et al. Realizing the promise of reverse phase protein arrays for clinical, translational, and basic research: A workshop report: The RPPA (Reverse Phase Protein Array) society. *Mol. Cell. Proteom.* **2014**, *13*, 1625–1643. [CrossRef] [PubMed]

43. Troncale, S.; Barbet, A.; Coulibaly, L.; Henry, E.; He, B.; Barillot, E.; Dubois, T.; Hupe, P.; de Koning, L. NormaCurve: A SuperCurve-based method that simultaneously quantifies and normalizes reverse phase protein array data. *PLoS ONE* **2012**, *7*, e38686. [CrossRef]

44. Gagaoua, M.; Monteils, V.; Couvreur, S.; Picard, B. Beef Tenderness Prediction by a Combination of Statistical Methods: Chemometrics and Supervised Learning to Manage Integrative Farm-To-Meat Continuum Data. *Foods* **2019**, *8*, 274. [CrossRef]
45. Gagaoua, M.; Picard, B.; Soulat, J.; Monteils, V. Clustering of sensory eating qualities of beef: Consistencies and differences within carcass, muscle, animal characteristics and rearing factors. *Livest. Sci.* **2018**, *214*, 245–258. [CrossRef]

46. Rousseeuw, P.J. Silhouettes: A graphical aid to the interpretation and validation of cluster analysis. *J. Comput. Appl. Math.* **1987**, *20*, 53–65. [CrossRef]

47. Gagaoua, M.; Listrat, A.; Andueza, D.; Gruffat, D.; Normand, J.; Mairesse, G.; Picard, B.; Hocquette, J.-F. What are the drivers of beef sensory quality using metadata of intramuscular connective tissue, fatty acids and muscle fiber characteristics? *Livest. Sci.* **2020**. [CrossRef]

48. Bernard, C.; Cassar-Malek, I.; Le Cunff, M.; Dubroeucq, H.; Renand, G.; Hocquette, J.F. New indicators of beef sensory quality revealed by expression of specific genes. *J. Agric. Food Chem.* **2007**, *55*, 5229–5237. [CrossRef] [PubMed]

49. D'Alessandro, A.; Marrocco, C.; Rinalducci, S.; Mirasole, C.; Failla, S.; Zolla, L. Chianina beef tenderness investigated through integrated Omics. *J. Proteom.* **2012**, *75*, 4381–4398. [CrossRef] [PubMed]

50. Guillemin, N.; Bonnet, M.; Jurie, C.; Picard, B. Functional analysis of beef tenderness. *J. Proteom.* **2011**, *75*, 352–365. [CrossRef] [PubMed]

51. Couvreur, S.; Le Bec, G.; Micol, D.; Picard, B. Relationships Between Cull Beef Cow Characteristics, Finishing Practices and Meat Quality Traits of Longissimus thoracis and Rectus abdominis. *Foods* **2019**, *8*, 141. [CrossRef]

52. Gagaoua, M.; Couvreur, S.; Le Bec, G.; Aminot, G.; Picard, B. Associations among Protein Biomarkers and pH and Color Traits in Longissimus thoracis and Rectus abdominis Muscles in Protected Designation of Origin Maine-Anjou Cull Cows. *J. Agric. Food Chem.* **2017**, *65*, 3569–3580. [CrossRef]

53. Gagaoua, M.; Terlouw, E.M.C.; Micol, D.; Hocquette, J.F.; Moloney, A.P.; Nuernberg, K.; Bauchart, D.; Boudjellal, A.; Scollan, N.D.; Richardson, R.I.; et al. Sensory quality of meat from eight different types of cattle in relation with their biochemical characteristics. *J. Integr. Agric.* **2016**, *15*, 1550–1563. [CrossRef]

54. Picard, B.; Gagaoua, M.; Micol, D.; Cassar-Malek, I.; Hocquette, J.F.; Terlouw, C.E. Inverse relationships between biomarkers and beef tenderness according to contractile and metabolic properties of the muscle. *J. Agric. Food Chem.* **2014**, *62*, 9808–9818. [CrossRef]

55. Golenhofen, N.; Perng, M.D.; Quinlan, R.A.; Drenckhahn, D. Comparison of the small heat shock proteins alphaB-crystallin, MKBP, HSP25, HSP20, and cvHSP in heart and skeletal muscle. *Histochem. Cell Biol.* **2004**, *122*, 415–425. [CrossRef]

56. Lomiwes, D.; Farouk, M.M.; Wiklund, E.; Young, O.A. Small heat shock proteins and their role in meat tenderness: A review. *Meat Sci.* **2014**, *96*, 26–40. [CrossRef] [PubMed]

57. Picard, B.; Kammoun, M.; Gagaoua, M.; Barboiron, C.; Meunier, B.; Chambon, C.; Cassar-Malek, I. Calcium Homeostasis and Muscle Energy Metabolism Are Modified in HspB1-Null Mice. *Proteomes* **2016**, *4*, 17. [CrossRef] [PubMed]

58. Dreiza, C.M.; Komalavilas, P.; Furnish, E.J.; Flynn, C.R.; Sheller, M.R.; Smoke, C.C.; Lopes, L.B.; Brophy, C.M. The small heat shock protein, HSPB6, in muscle function and disease. *Cell Stress Chaperones* **2010**, *15*, 1–11. [CrossRef] [PubMed]

59. Rembold, C.M.; Foster, D.B.; Strauss, J.D.; Wingard, C.J.; Eyk, J.E. cGMP-mediated phosphorylation of heat shock protein 20 may cause smooth muscle relaxation without myosin light chain dephosphorylation in swine carotid artery. *J. Physiol. 2000*, *524 Pt 3*, 865–878. [CrossRef]

60. Lomiwes, D.; Hurst, S.M.; Dobbie, P.; Frost, D.A.; Hurst, R.D.; Young, O.A.; Farouk, M.M. The protection of bovine skeletal myofibrils from proteolytic damage post mortem by small heat shock proteins. *Meat Sci.* **2014**, *97*, 548–557. [CrossRef]

61. Ma, D.; Kim, Y.H.B. Proteolytic changes of myofibrillar and small heat shock proteins in different bovine muscles during aging: Their relevance to tenderness and water-holding capacity. *Meat Sci.* **2020**, *163*, 108090. [CrossRef]

62. Kamradt, M.C.; Chen, F.; Cryns, V.L. The Small Heat Shock Protein αB-Crystallin Negatively Regulates Cytochrome c- and Caspase-8-dependent Activation of Caspase-3 by Inhibiting Its Autoproteolytic Maturation. *J. Biol. Chem.* **2001**, *276*, 16059–16063. [CrossRef]

63. Kamradt, M.C.; Chen, F.; Sam, S.; Cryns, V.L. The Small Heat Shock Protein αB-crystallin Negatively Regulates Apoptosis during Myogenic Differentiation by Inhibiting Caspase-3 Activation. *J. Biol. Chem.* **2002**, *277*, 38731–38736. [CrossRef]
64. Pulford, D.J.; Fraga Vazquez, S.; Frost, D.F.; Fraser-Smith, E.; Dobbie, P.; Rosenvold, K. The intracellular distribution of small heat shock proteins in post-mortem beef is determined by ultimate pH. *Meat Sci.* **2008**, *79*, 623–630. [CrossRef]

65. Pulford, D.J.; Dobbie, P.; Fraga Vazquez, S.; Fraser-Smith, E.; Frost, D.A.; Morris, C.A. Variation in bull beef quality due to ultimate muscle pH is correlated to endopeptidase and small heat shock protein levels. *Meat Sci.* **2009**, *83*, 1–9. [CrossRef]

66. D'Alessandro, A.; Rinalducci, S.; Marrocco, C.; Zolla, V.; Napolitano, F.; Zolla, L. Love me tender: An Omics window on the bovine meat tenderness network. *J. Proteom.* **2012**, *75*, 4360–4380. [CrossRef] [PubMed]

67. Polati, R.; Menini, M.; Robotti, E.; Millioni, R.; Marengo, E.; Novelli, E.; Balzan, S.; Cecconi, D. Proteomic changes involved in tenderization of bovine Longissimus dorsi muscle during prolonged ageing. *Food Chem.* **2012**, *135*, 2052–2069. [CrossRef] [PubMed]

68. Zapata, I.; Zerby, H.N.; Wick, M. Functional proteomic analysis predicts beef tenderness and the tenderness differential. *J. Agric. Food Chem.* **2009**, *57*, 4956–4963. [CrossRef] [PubMed]

69. Morzel, M.; Terlouw, C.; Chambon, C.; Micol, D.; Picard, B. Muscle proteome and meat eating qualities of Longissimus thoracis of "Blonde d'Aquitaine" young bulls: A central role of HSP27 isoforms. *Meat Sci.* **2008**, *78*, 297–304. [CrossRef] [PubMed]

70. Carvalho, M.E.; Gasparin, G.; Poleti, M.D.; Rosa, A.F.; Balieiro, J.C.; Labate, C.A.; Nassu, R.T.; Tullio, R.R.; Regitano, L.C.; Mourao, G.B.; et al. Heat shock and structural proteins associated with meat tenderness in Nellore beef cattle, a Bos indicus breed. *Meat Sci.* **2014**, *96*, 1318–1324. [CrossRef]

71. Thornton, K.J.; Chapalamadugu, K.C.; Eldredge, E.M.; Murdoch, G.K. Analysis of Longissimus thoracis Protein Expression Associated with Variation in Carcass Quality Grade and Marbling of Beef Cattle Raised in the Pacific Northwestern United States. *J. Agric. Food Chem.* **2017**, *65*, 1434–1442. [CrossRef]

72. DeLany, J.P.; Floyd, Z.E.; Zvonic, S.; Smith, A.; Gravois, A.; Reiners, E.; Wu, X.; Kilroy, G.; Lefevre, M.; Gimble, J.M. Proteomic Analysis of Primary Cultures of Human Adipose-derived Stem Cells. *Modul. Adipogenes.* **2005**, *4*, 731–740. [CrossRef]

73. Takeda, K.; Sriram, S.; Chan, X.H.D.; Ong, W.K.; Yeo, C.R.; Tan, B.; Lee, S.-A.; Kong, K.V.; Hoon, S.; Jiang, H.; et al. Retinoic Acid Mediates Visceral-Specific Adipogenic Defects of Human Adipose-Derived Stem Cells. *Diabetes* **2016**, *65*, 1164–1178. [CrossRef]

74. Grabez, V.; Kathri, M.; Phung, V.; Moe, K.M.; Slinde, E.; Skaugen, M.; Saarem, K.; Egelandsdal, B. Protein expression and oxygen consumption rate of early postmortem mitochondria relate to meat tenderness. *J. Anim. Sci.* **2015**, *93*, 1967–1979. [CrossRef]

75. Jia, X.; Ekman, M.; Grove, H.; Faergestad, E.M.; Aass, L.; Hildrum, K.I.; Hollung, K. Proteome changes in bovine longissimus thoracis muscle during the early postmortem storage period. *J. Proteome Res.* **2007**, *6*, 2720–2731. [CrossRef] [PubMed]

76. Wu, W.; Yu, Q.Q.; Fu, Y.; Tian, X.J.; Jia, F.; Li, X.M.; Dai, R.T. Towards muscle-specific meat color stability of Chinese Luxi yellow cattle: A proteomic insight into post-mortem storage. *J. Proteom.* **2016**, *147*, 108–118. [CrossRef] [PubMed]

77. Vasiliiou, V.; Thompson, D.C.; Smith, C.; Fujita, M.; Chen, Y. Aldehyde dehydrogenases: From eye crystallins to metabolic disease and cancer stem cells. *Chem. Biol. Interact.* **2013**, *202*, 2–10. [CrossRef] [PubMed]

78. Dominguez, R.; Pateiro, M.; Gagaoua, M.; Barba, F.J.; Zhang, W.; Lorenzo, J.M. A Comprehensive Review on Lipid Oxidation in Meat and Meat Products. *Antioxidants* **2019**, *8*, 429. [CrossRef] [PubMed]

79. Sasaki, K.; Motoyama, M.; Narita, T.; Hagi, T.; Ojima, K.; Oe, M.; Nakajima, I.; Kitsunai, K.; Saito, Y.; Hatori, H.; et al. Characterization and classification of Japanese consumer perceptions for beef tenderness using descriptive texture characteristics assessed by a trained sensory panel. *Meat Sci.* **2014**, *96*, 994–1002. [CrossRef]

80. Bonnet, M.; Faulconnier, Y.; Leroux, C.; Jurie, C.; Cassar-Malek, I.; Bauchart, D.; Boulesteix, P.; Pethick, D.; Hocquette, J.F.; Chilliard, Y. Glucose-6-phosphate dehydrogenase and leptin are related to marbling differences among Limousin and Angus or Japanese Black × Angus steers1,2. *J. Anim. Sci.* **2007**, *85*, 2882–2894. [CrossRef]

81. Ertbjerg, P.; Puolanne, E. Muscle structure, sarcomere length and influences on meat quality: A review. *Meat Sci.* **2017**, *132*, 139–152. [CrossRef]

82. Lana, A.; Zolla, L. Proteolysis in meat tenderization from the point of view of each single protein: A proteomic perspective. *J. Proteom.* **2016**, *147*, 85–97. [CrossRef]
83. Picard, B.; Gagaoua, M. Muscle Fiber Properties in Cattle and Their Relationships with Meat Qualities: An Overview. *J. Agric. Food Chem.* **2020**, *68*, 6021–6039. [CrossRef]

84. Hocquette, J.-F.; Jurie, C.; Picard, B.; Alberti, P.; Panea, B.; Christensen, M.; Failla, S.; Gigli, S.; Levéziel, H.; Olleta, J.L.; et al. Metabolic and contractile characteristics of Longissimus thoracis muscle of young bulls from 15 European breeds in relationship with body composition. In Proceedings of the International Symposium on Energy and Protein Metabolism and Nutrition (ISEP 2007), Vichy, France, 9–13 September 2007; pp. 111–112.

85. Poleti, M.D.; Regitano, L.C.A.; Souza, G.H.M.F.; Cesar, A.S.M.; Simas, R.C.; Silva-Vignato, B.; Oliveira, G.B.; Andrade, S.C.S.; Cameron, L.C.; Coutinho, L.L. Longissimus dorsi muscle label-free quantitative proteomic reveals biological mechanisms associated with intramuscular fat deposition. *J. Proteom.* **2018**, *179*, 30–41. [CrossRef]

86. Poleti, M.D.; Regitano, L.C.A.; Souza, G.H.M.F.; Cesar, A.S.M.; Simas, R.C.; Silva-Vignato, B.; Montenegro, H.; Pértille, F.; Balieiro, J.C.C.; Cameron, L.C.; et al. Proteome alterations associated with the oleic acid and cis-9, trans-11 conjugated linoleic acid content in bovine skeletal muscle. *J. Proteom.* **2020**, *222*, 103792. [CrossRef]

87. Zhao, C.; Zan, L.; Wang, Y.; Scott Updike, M.; Liu, G.; Bequette, B.J.; Baldwin Vi, R.L.; Song, J. Functional proteomic and interactome analysis of proteins associated with beef tenderness in Angus cattle. *Livest. Sci.* **2014**, *161*, 201–209. [CrossRef]

88. Beldarrain, L.R.; Aldai, N.; Picard, B.; Sentandreu, E.; Navarro, J.L.; Sentandreu, M.A. Use of liquid isoelectric focusing (OFFGEL) on the discovery of meat tenderness biomarkers. *J. Proteom.* **2018**, *183*, 25–33. [CrossRef]

89. Ouali, A. Meat Tenderization: Possible Causes and Mechanisms. A Review. *J. Muscle Foods* **1990**, *1*, 129–165. [CrossRef]

90. Shathasivam, T.; Kislinger, T.; Gramolini, A.O. Genes, proteins and complexes: The multifaceted nature of FHL family proteins in diverse tissues. *J. Cell. Mol. Med.* **2010**, *14*, 2702–2720. [CrossRef]

91. Laville, E.; Sayd, T.; Morzel, M.; Blinet, S.; Chambon, C.; Lepetit, J.; Renand, G.; Hocquette, J.F. Proteome changes during meat aging in tough and tender beef suggest the importance of apoptosis and protein solubility for beef aging and tenderization. *J. Agric. Food Chem.* **2009**, *57*, 10755–10764. [CrossRef]

92. Bjarnadottir, S.G.; Hollung, K.; Hoy, M.; Bendixen, E.; Codrea, M.C.; Veiseth-Kent, E. Changes in protein abundance between tender and tough meat from bovine longissimus thoracis muscle assessed by isobaric Tag for Relative and Absolute Quantitation (iTRAQ) and 2-dimensional gel electrophoresis analysis. *J. Anim. Sci.* **2012**, *90*, 2035–2043. [CrossRef]

93. Boudon, S.; Ounaissi, D.; Viala, D.; Monteils, V.; Picard, B.; Cassar-Malek, I. Label free shotgun proteomics for the identification of protein biomarkers for beef tenderness in muscle and plasma of heifers. *J. Proteom.* **2020**, *217*, 103685. [CrossRef]

94. Wang, Y.H.; Bower, N.I.; Reverter, A.; Tan, S.H.; De Jager, N.; Wang, R.; McWilliam, S.M.; Cafe, L.M.; Greenwood, P.L.; Lehnert, S.A. Gene expression patterns during intramuscular fat development in cattle1. *J. Anim. Sci.* **2009**, *87*, 119–130. [CrossRef]

95. Pillar, N.; Pleniceanu, O.; Fang, M.; Ziv, L.; Lahav, E.; Botchan, S.; Cheng, L.; Dekel, B.; Shomron, N. A rare variant in the FHL1 gene associated with X-linked recessive hypoparathyroidism. *Hum. Genet.* **2017**, *136*, 835–845. [CrossRef]

96. Carafoli, E. Calcium signaling: A tale for all seasons. *Proc. Natl. Acad. Sci. USA* **2002**, *99*, 1115–1122. [CrossRef]

97. Demaurex, N.; Distelhorst, C. Cell biology. Apoptosis—the calcium connection. *Science* **2003**, *300*, 65–67. [CrossRef]

98. Lee, J.Y.; Lori, D.; Wells, D.J.; Kemp, P.R. FHL1 activates myostatin signalling in skeletal muscle and promotes atrophy. *FEBS Open Bio* **2015**, *5*, 753–762. [CrossRef]

99. Malheiros, J.M.; Braga, C.P.; Grove, R.A.; Ribeiro, F.A.; Calkins, C.R.; Adamec, J.; Chardulo, L.A.L. Influence of oxidative damage to proteins on meat tenderness using a proteomics approach. *Meat Sci.* **2019**, *148*, 64–71. [CrossRef]

100. Mato, A.; Rodríguez-Vázquez, R.; López-Pedrouso, M.; Bravo, S.; Franco, D.; Zapata, C. The first evidence of global meat phosphoproteome changes in response to pre-slaughter stress. *BMC Genom.* **2019**, *20*, 590. [CrossRef]

101. Rosa, A.F.; Moncau, C.T.; Poleti, M.D.; Fonseca, L.D.; Balieiro, J.C.C.; Silva, S.L.E.; Eler, J.P. Proteome changes of beef in Nellore cattle with different genotypes for tenderness. *Meat Sci.* **2018**, *138*, 1–9. [CrossRef]
102. Silva, L.H.P.; Rodrigues, R.T.S.; Assis, D.E.F.; Benedeti, P.D.B.; Duarte, M.S.; Chizzotti, M.L. Explaining meat quality of bulls and steers by differential proteome and phosphoproteome analysis of skeletal muscle. *J. Proteom.* **2019**, *199*, 51–66. [CrossRef]

103. Marino, R.; Albenzio, M.; Della Malva, A.; Caroprese, M.; Santillo, A.; Sevi, A. Changes in meat quality traits and sarcoplasmic proteins during aging in three different cattle breeds. *Meat Sci.* **2014**, *98*, 178–186. [CrossRef]

104. Didiasova, M.; Schaefer, L.; Wygrecka, M. When Place Matters: Shuttling of Enolase-1 across Cellular Compartments. *Front. Cell Dev. Biol.* **2019**, *7*, 61. [CrossRef]

105. Sedoris, K.C.; Thomas, S.D.; Miller, D.M. Hypoxia induces differential translation of enolase/MBP-1. *BMC Cancer* **2010**, *10*, 157. [CrossRef]

106. Fisher, A.B. Peroxiredoxin 6 in the repair of peroxidized cell membranes and cell signaling. *Arch. Biochem. Biophys.* **2017**, *617*, 68–83. [CrossRef]

107. Arriga, R.; Pacifici, F.; Capuani, B.; Coppola, A.; Orlandi, A.; Scioli, M.G.; Pastore, D.; Andreadi, A.; Sbraccia, P.; Tesauro, M.; et al. Peroxiredoxin 6 Is a Key Antioxidant Enzyme in Modulating the Link between Glycemic and Lipogenic Metabolism. *Oxidative Med. Cell. Longev.* **2019**, *2019*, 9685607. [CrossRef]

108. Cai, C.; Masumiya, H.; Weisleder, N.; Matsuda, N.; Nishi, M.; Hwang, M.; Ko, J.-K.; Lin, P.; Thornton, A.; Zhao, X.; et al. MG53 nucleates assembly of cell membrane repair machinery. *Nat. Cell Biol.* **2009**, *11*, 56–64. [CrossRef]

109. Jung, S.-Y.; Ko, Y.-G. TRIM72, a novel negative feedback regulator of myogenesis, is transcriptionally activated by the synergism of MyoD (or myogenin) and MEF2. *Biochem. Biophys. Res. Commun.* **2010**, *396*, 238–245. [CrossRef]

<figure><img src="image_32.png" /></figure>

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
<figure><img src="image_33.png" /></figure>
Article

# Proteomic Changes in Sarcoplasmic and Myofibrillar Proteins Associated with Color Stability of Ovine Muscle during Post-Mortem Storage

Xiaoguang Gao <fnref n="1" />/<fnref n="2" />, Dandan Zhao <fnref n="1" />, Lin Wang <fnref n="1" />, Yue Cui <fnref n="1" />, Shijie Wang <fnref n="1" />, Meng Lv <fnref n="1" />, Fangbo Zang <fnref n="1" /> and Ruitong Dai <fnref n="2" />,\*

**Citation:** Gao, X.; Zhao, D.; Wang, L.; Cui, Y.; Wang, S.; Lv, M.; Zang, F.; Dai, R. Proteomic Changes in Sarcoplasmic and Myofibrillar Proteins Associated with Color Stability of Ovine Muscle during Post-Mortem Storage. *Foods* **2021**, *10*, 2989. https://doi.org/10.3390/foods10122989

Academic Editors: Mohammed Gagaoua and Brigitte Picard

Received: 1 November 2021

Accepted: 1 December 2021

Published: 3 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_34.png" /></figure>

Copyright: © 2021 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).

**Abstract:** The objective of this study was to investigate the proteomic characteristics for the sarcoplasmic and myofibrillar proteomes of *M. longissimus lumborum* (LL) and *M. psoasmajor* (PM) from *Small-tailed Han* Sheep. During post-mortem storage periods (1, 3, and 5 days), proteome analysis was applied to elucidate sarcoplasmic and myofibrillar protein changes in skeletal muscles with different color stability. Proteomic results revealed that the identified differentially abundant proteins were glycolytic enzymes, energy metabolism enzymes, chaperone proteins, and structural proteins. Through Pearson's correlation analysis, a few of those identified proteins (Pyruvate kinase, Adenylate kinase isoenzyme 1, Creatine kinase M-type, and Carbonic anhydrase 3) were closely correlated to representative meat color parameters. Besides, bioinformatics analysis of differentially abundant proteins revealed that the proteins mainly participated in glycolysis and energy metabolism pathways. Some of these proteins may have the potential probability to be predictors of meat discoloration during post-mortem storage. Within the insight of proteomics, these results accumulated some basic theoretical understanding of the molecular mechanisms of meat discoloration.

**Keywords:** ovine; color stability; enzyme; proteomics; bioinformatics

## 1. Introduction

Among various sensory characteristics, the color of fresh meat is a critical quality factor, influencing the purchase decision of the consumers [1]. The defects of meat color have always been connected to spoilage and unwholesomeness [2]. Prevention of meat discoloration has always been a challenging task and long-term objective for scientific researchers [3–5]. Consequently, maintaining the “cherry-red” color of fresh meat is very critical for the meat industry.

The relative contents of the reduced deoxymyoglobin (DeoMb, purple), the oxygenated oxymyoglobin (OxyMb, bright red), and the oxidized metmyoglobin (MetMb, brown) determine the color of fresh meat [6]. The reason for meat discoloration is most commonly attributed to the formation of MetMb. Furthermore, the affecting factors of redox forms of myoglobin are formed by many intrinsic and extrinsic ones such as breed, muscle type, temperature, and oxygen partial pressure of storage [7,8]. Because of the distinct metabolic function and biochemical profile, *M. longissimus lumborum* (LL) and *M. psoas major* (PM) were considered to be typical color-stable and color-labile muscles, respectively, by previous researchers [9,10]. Thus, both of the skeletal muscles whose color stabilities (color-stable and color-labile) were opposing can be used as typical experimental subjects to elaborate the underlying biochemical mechanisms of post-mortem

<footnote n="1">¹ College of Food Science and Biology, Hebei University of Science and Technology, No. 26 Yuxiang Street, Yuhua District, Shijiazhuang 050000, China; gaoxiaoguang23@hotmail.com (X.G.); zdd6364@126.com (D.Z.); w592148022@163.com (L.W.); cuiy1223@hotmail.com (Y.C.); mrshjwang@163.com (S.W.); lmeng1126@163.com (M.L.); 15076005611@163.com (F.Z.)</footnote>

<footnote n="2">² College of Food Science and Nutritional Engineering, China Agricultural University, No. 17 Qinghua East Road, Haidian District, Beijing 100083, China</footnote>

<footnote>\* Correspondence: dairuitong@hotmail.com; Tel.: +86-010-62737547</footnote>
color stability [11]. The post-mortem changes in muscles constitute a very complicated biochemical processes influencing meat color stability, as skeletal muscle consist of proteins, lipids, and other biomacromolecules and micro-molecules [11,12]. Proteomics provides an efficient way, in post-genome era, to illustrate the post-mortem changes in meat quality development, including the field of color characteristics [13,14]. In the previous literature, some investigators explored a series of post-mortem proteomic changes in beef and pork [15–18]. However, data of physiological and biochemical mechanisms of ovine color stability during post-mortem storage has rarely been reported [7,19]. The proteomics for the variations of different ovine muscles (color-stable and color-labile) still need to be investigated. Therefore, the potential molecular mechanisms need further research. The objective of this study was to investigate the proteomic characteristics for the sarcoplasmic and myofibrillar proteome of *M. longissimus lumborum* (LL) and *M. psoasmajor* (PM) from *Small-tailed Han* Sheep during post-mortem storage.

## 2. Materials and Methods

### 2.1. Sample Preparation

The skeletal muscle samples were harvested from 6 male Small-tailed Han Sheep, slaughtered at 8 months, with a mean carcass weight of 15.6 ± 0.2 kg. There were 2 typical types of skeletal muscles (*M. longissimus lumborum* (LL) and *M. psoas major* (PM)) that were excised (30 min after slaughter) from both sides of each carcass (n = 6 carcasses). The whole process of slaughter followed the industrial practice. The animals received humanely non-painful manipulations at the termination of the procedure, without regaining consciousness. Muscles samples were cut along the direction perpendicular to the muscle fibers. For each sample (corresponding to 1 carcass and 1 type of muscle), 6 equal-weight (30 g, trimmed free of connective tissue and fat) pieces were randomly assigned to 12 replications [n = 6, (6 muscles per type × 2 sides × 3 sections per muscle)/3 days].

The samples were wrapped individually with polyethylene films (350–400 cm³m⁻²h⁻¹atm⁻¹, Mitsui Chemical, Japan) and put in a refrigerator for a storage time of 5 days (4 ± 1 °C). Muscles were taken out for analyses at day 1, day 3, and day 5. At each time point (day 1, day 3, and day 5), instrumental color parameters and metmyoglobin reducing activity (MRA) of the muscles were measured. Sarcoplasmic and myofibrillar proteins were extracted for comparing the change in the early stage of storage (day 1) and the late stage of storage (day 5). At the storage time of day 3, proteins were not extracted accordingly. At the end of storage period, myofibrillar proteins were extracted for the comparison of proteomic difference between muscle types (LL and PM).

### 2.2. Instrumental Color

All the measurements of color parameters were determined on the surface of each sample. The instrumental color parameters were measured by utilizing a Minolta chromameter (CR-400, Minolta Inc., Osaka, Japan), expressing as CIE Lab lightness (L*), redness (a*), and yellowness (b*) with D₆₅ standard illuminant. Accordingly, the parameters were calibrated through the white and black reference standards. The observer angle and measurement area were 2° and 8 mm, respectively. Hue angle (h°) and chroma (C*) were determined by the following formula: h° = [ATAN(b*/a*) × (180/π)], C* = [(a*² + b*²)⁰.⁵] [20].

LL and PM samples (5 g of each sample) were collected and homogenized for 10 s in a phosphate buffer (25 mL, pH 6.8, 0 °C, 40 mM) by using a homogenizer (F6-10, Fluko, Shanghai, China). The homogenate solution was placed in 4 °C for an hour and centrifuged (4500 × g, 4 °C). Supernatant fluid was collected and filtered (Whatman-No.1 filter paper). The absorbance of the filtrate was tested with a spectrometer (wavelength of 503, 525, 557, 582 nm; TU-1810, PERSEE, Beijing, China). The proportion of myoglobin redox forms was calculated as follows [21]:

$$ 
\begin{align*}
[\text{DeoMb}] &= C_{\text{DeoMb}} \div C_{\text{Mb}} = -0.543R_1 + 1.594R_2 + 0.552R_3 - 1.3290 \\
[\text{OxyMb}] &= C_{\text{OxyMb}} \div C_{\text{Mb}} = 0.722R_1 - 1.432R_2 - 1.659R_3 + 2.599 \\
[\text{MetMb}] &= C_{\text{MetMb}} \div C_{\text{Mb}} = -0.159R_1 - 0.085R_2 + 1.262R_3 - 0.520
\end{align*}
\quad (1) $$
where R₁ = A₅₈₂ ÷ A₅₂₅; R₂ = A₅₅₇ ÷ A₅₂₅; R₃ = A₅₀₃ ÷ A₅₂₅.

## 2.3. Metmyoglobin Reducing Activity

The metmyoglobin reducing activity (MRA) was determined through spectrophotography and based on Mikkelsen et al. [22], including the following experimental equipment: spectrometer (TU-1810, PERSEE, Beijing, China), refrigeration centrifuge (3K15, SIGMA, Hamburg, Germany) and homogenizer (F6-10, Fluko, Shanghai, China). Muscle sample (12 g, connective tissue and fat removed) was homogenized (20 mL, 2.0 mM phosphate buffer, pH 7.0) for 30 s. The homogenate was centrifuged (35,000× g, 4 °C) for 30 min. The supernatant fluid was collected and filtered (Whatman-No.1 filter paper) for further removal of fat.

K₃Fe(CN)₆ was used to oxidize OxyMb into MetMb. The solution was dialyzed against phosphate buffer (2.0 mM, pH 7.0, 4 °C, 14,000 Mw cut-off) and centrifuged for 20 min (15,000× g, 4 °C). The supernatant was adjusted to volume of 20 mL through phosphate buffer (2.0 mM, pH 7.0). The standard assay (pH 6.4, 25 °C) was the mixed solution of NADH (0.1 mL, 2.0 mM, used for initiating the reaction), K₄[Fe(CN)₆] (0.1 mL, 3.0 mM), EDTA (0.1 mL, 5.0 mM), Mb Fe(III) (0.2 mL of 0.75 mM within 2.0 mM phosphate buffer, pH 7.0), phosphate buffer (0.1 mL, 50 mM, pH 7.0), deionized water (0.1 mL) and muscle extract (0.3 mL). The blank contained all chemicals except for NADH (replace NADH with water, no reaction). MRA was followed the changes in absorbance (580 nm) and determined as nmol MetMb reduced (min⁻¹ g⁻¹).

## 2.4. Sarcoplasmic and Myofibrillar Protein Extraction

Sarcoplasmic and myofibrillar protein extraction methods were based on Sayd et al. [23] and Kim et al. [24] with a minor modification, respectively. Each muscle sample (3 g) was homogenized with extraction buffer (30 mL, 2 mM EDTA, 40 mM Tris, pH 7.0, 4 °C) containing a protease inhibitors cocktail by utilizing a homogenizer (F6-10, Fluko, Shanghai, China). The sample was centrifuged (10 min, 10,000× g, 4 °C) and the supernatant (sarcoplasmic proteins) was collected for further experiments (stored in -80 °C). The precipitate of each sample was incubated for 40 min in 2 mL of extraction buffer (1% pH 3–10 bio-lyte ampholytes (Bio-Rad, Hercules, CA, USA), 2 M thio-urea, 8 M urea, 2% CHAPS (w/v), 65 mM DTT, with protease inhibitors cocktail). The sample was centrifuged (40,000× g, 4 °C) for 60 min and the supernatant (myofibrillar proteins) was stored at -80 °C for further experiments.

## 2.5. Two-Dimensional Electrophoresis and Gel Image Analysis

The sarcoplasmic and myofibrillar proteome (800 µg, respectively) was included in the Bio-Rad buffer (8 M urea, 2% CHAPS (w/v), 50 mM DTT, 0.3% carrier ampholyte (v/v), bromophenol blue). The samples were loaded on the IPG strips (immobilized pH gradient, 24 cm, pH 3–10, Bio-Rad). The PROTEAN II XL and PROTEAN IEF (isoelectric focusing) cell system (Bio-Rad, Hercules, CA, USA) performed in the first and second dimension electrophoreses, respectively.

In the first dimension, IEF was subjected to passive rehydration (16 h) and then rapid voltage ramping (80,000 V h) was applied. In the second dimension, proteins were resolved on 12% SDS-PAGE gels, which were stained in Coomassie Brilliant Blue (48 h). Gel images were scanned by an image scanner (ImageScanner III, GE Healthcare, Branford, CT, USA). ImageMaster 2D Platinum software (6.0 version, GE Healthcare, Branford, CT, USA) was used to analyze the scanning images. Through expressing the relative number of each spot as the ratio of the number of single spots to the total number of valid spots, the detected and matched spots were normalized. The mean values of gels for each sample and spot were calculated in triplicate. A spot was considered as differential proteins when it came to 5% statistical significance (p < 0.05) in one-way ANOVA.
## 2.6. Identification of Protein Spots by Mass Spectrometry

The protein spots were carefully excised from gels and then de-stained for 30 min in wash buffer (100 µL, 25 mM NH₄HCO₃/50% acetonitrile (v/v)). Washed-out gel-spots (dehydrated in 100% acetonitrile) were dried completely with a centrifuge (Vacufuge plus, Eppendorf, Hamburg, Germany) and then incubated in 15 ng/µL trypsin and 25 mM NH₄HCO₃ (37 °C, 16 h). The peptides were incubated in trifluoroacetic acid (20 µL, 0.1% (v/v), 37 °C, 40 min) after digestion. The above extraction procedures were repeated by using 50% acetonitrile/0.1% trifluoroacetic acid (v/v). Sediments were washed in trifluoroacetic acid and then vacuum freeze-dried for further analysis. The proteins were identified by using AUTOFLEX II TOF-TOF mass spectrometer (autoflex™ speed, Bruker Daltonik, Bremen, Germany). Samples in 1 µL of buffer (50% acetonitrile and 5 mg/mL α-CHCA in 0.1% trifluoroacetic acid) were loaded on plate (AnchorChip, 384-MPT).

Each crystallized sample was washed by using 0.1% trifluoroacetic acid for removing salt ions. Protein identification was performed by peptide mass fingerprinting (PMF), searching the Mascot (2.2 version, Matrix Science, London, UK), and matched with a sheep (*Ovis aries*) family in the Uniprot database (https://www.uniprot.org/ (accessed on 8 October 2021)), correspondingly.

## 2.7. Statistical Analysis

Proportions of MetMb, OxyMb, DeoMb, and instrumental color attributes (L*, a*, b* and MRA) were analyzed separately. Random terms for all models included the carcass and the processing day. For color models, an additional term for side nested within the carcass was added to the random model to account for the repetition of measures [25]. The correlation between meat color attributes and differentially abundant proteins was analyzed by Pearson's correlation. The mixed models, which were prediction models for muscle color traits, were based on Starkey et al.

Data were expressed as mean ± SE (standard error). The GLM (general linear model) procedure was used to analyze the instrumental color and MRA data with the SAS 8.2 software (Institute Inc., Cary, NC, USA). Means and the correlation data were analyzed by one-way ANOVA by using the SPSS 20.0 software (IBM Inc., Chicago, IL, USA) with the 5% (*p* < 0.05) level of statistical significance difference.

Protein-protein interaction (PPI), gene ontology (GO), and Kyoto encyclopedia of genes and genomes (KEGG) were applied to characterize the functional information of the identified proteins. The protein-protein interaction was analyzed by String 10.0. Bioinformatics analysis was performed through DAVID Bioinformatics Resources 6.7. The species of *Ovis aries* was selected.

# 3. Results and Discussion

## 3.1. Instrumental Color and MRA

The instrumental color parameters of skeletal muscles (*M. longissimus lumborum* (LL) and *M. psoas major* (PM)) are presented in Table 1. There was an insignificant (*p* > 0.05) decrease trend in L*-value. Besides, no significant difference (*p* > 0.05) was found between these two kinds of skeletal muscles within storage time points in L*-value. Although some previous studies [26] showed that the variation in L*-value were very subtle. Besides being a color parameter, the L*-value is also considered an indicator of the water-holding capacity of meat [11]. During ageing in ovine muscle, structural changes of proteins were considered one of the causes of the change in lightness (L*-value) [27].

For the other color parameters, there were decreasing trends in redness (a*-value), yellowness (b*-value), and chroma (C*-value) of both kinds of muscles. Meanwhile, hue angle (h°-value) exhibited increasing trends during 5 days of post-mortem storage. Those results were consistent with studies in the previous literature [26,28,29]. In an associated study [17], the relative content of OxyMb, MetMb, and DeoMb, on the value of the color parameters of minced pork loin, were evaluated. The color parameters were associated with myoglobin redox forms. Our results (Figure 1) were partially in agreement with Karamucki
et al. [17]. For instance, an increase in OxyMb contributes most greatly to an increase in redness (a*-value), yellowness (b*-value), and chroma (C*-value). Chroma (C*-value) was deemed as an indicator for vividness of color [30]. Compared single color coordinate, hue angle (h°-value), showed more realistic perspectives on meat discoloration [8]. In the literature, muscles with unstable color often go along with greater hue angle (h°-value) [10]. In this study, compared with PM muscle, LL showed lower (p < 0.05) hue angle (h°-value) values, indicating a greater color stability.

<figure><img src="image_35.png" /><figcaption>**Figure 1.** Changes in the relative proportions of MetMb, OxyMb, and DeoMb in ovine *M. longissimus* lumborum (LL) and *M. psoas major* (PM) muscles during postmortem storage (4 ± 1 °C).</figcaption></figure>

<table><caption>**Table 1.** Instrumental color attributes of ovine *M. longissimus* lumborum (LL) and *M. psoas major* (PM) muscles at 1, 3, and 5 days of postmortem storage (4 °C).</caption><thead><tr><th rowspan="2">Attribute</th><th rowspan="2">Muscle</th><th colspan="3">Storage Time (Day)</th><th rowspan="2">SEM</th></tr><tr><th>1</th><th>3</th><th>5</th></tr></thead><tbody><tr><td rowspan="3">L*</td><td>LL</td><td>44.25 ax</td><td>43.61 ax</td><td>37.83 by</td><td>0.403</td></tr><tr><td>PM</td><td>44.81 bx</td><td>46.75 ax</td><td>44.70 bx</td><td>0.472</td></tr><tr><td>LL</td><td>21.41 ax</td><td>17.20 bx</td><td>13.62 cx</td><td>0.35</td></tr><tr><td rowspan="2">a*</td><td>PM</td><td>16.96 ay</td><td>12.98 by</td><td>10.59 cy</td><td>0.306</td></tr><tr><td>LL</td><td>15.57 ay</td><td>15.02 ax</td><td>12.72 by</td><td>0.262</td></tr><tr><td rowspan="2">b*</td><td>PM</td><td>17.33 ax</td><td>15.88 bx</td><td>13.28 cx</td><td>0.261</td></tr><tr><td>LL</td><td>26.47 ax</td><td>22.83 bx</td><td>18.64 cx</td><td>0.423</td></tr><tr><td rowspan="2">C*</td><td>PM</td><td>24.25 ay</td><td>20.51 by</td><td>16.99 cy</td><td>0.374</td></tr><tr><td>LL</td><td>36.03 cy</td><td>41.12 by</td><td>43.06 ay</td><td>0.322</td></tr><tr><td rowspan="2">h°</td><td>PM</td><td>45.62 bx</td><td>50.74 ax</td><td>51.43 ax</td><td>0.372</td></tr><tr><td>LL</td><td>0.210 ax</td><td>0.145 bx</td><td>0.067 cx</td><td>0.003</td></tr><tr><td>MRA</td><td>PM</td><td>0.170 ay</td><td>0.121 by</td><td>0.041 cy</td><td>0.004</td></tr></tbody></table>

Means in rows with different superscripts (a-c) are different (p < 0.05) and means in column for LL and PM for each storage time and attribute with different superscripts (x-z) are different (p < 0.05). (a) SEM, standard error of the mean. (b) C* (Chroma) = (a*² + b*²)⁰.⁵. (c) h° (Hue angle) = ATAN(b*/a*) × (180/π).

MRA (metmyoglobin reducing activity) decreased with storage time. LL showed greater metmyoglobin-reducing activity than PM (p < 0.05, Table 1). Muscles with more color stability had greater ability to reduce metmyoglobin into reduction state. In practice, the loss of reducing activity was influenced by many cofactors, such as the decline of
pH, the depletion of substrates, the loss of respiration enzymes or functional properties, and the structural integrity of mitochondria during post-mortem storage of meat [29,31]. MRA gradually decreased, partially due to depletion of the NADH pool and reduction in material consumption [32]. Consequently, LL muscle displayed greater color stability than PM during storage period.

## 3.2. Sarcoplasmic and Myofibrillar Proteome Analysis

The differentially abundant (*p* < 0.05) proteins, which were identified by two-dimensional electrophoresis image analysis of ovine muscles (*M. longissimus lumborum* (LL, color-stable) and *M. psoas major* (PM, color-labile)), were presented in Figure 2. Differential abundance changes of spots were at least twofold (*p* < 0.05). The identification and related information were shown in Table 2. The identified proteins were matched with the sheep (*Ovis aries*) family in the Uniprot database. According to physiological function in metabolism, the differentially abundant sarcoplasmic and myofibrillar proteins were classified into three categories: enzymes, chaperone proteins, and structural proteins.

Please note: The early stage of storage (day 1) of myofibrillar proteins gel images were not shown, as there was little differentially abundant proteins between LL and PM muscles at day 1.

<table><caption>**Table 2.** Differentially abundant proteins in ovine *M. longissimus lumborum* (LL) and *M. psoas major* (PM) muscles at day 1 or day 5 postmortem storage (4 °C).</caption><thead><tr><th>Spot.<br>(a)</th><th>Protein Name <sup>(b)</sup></th><th>Uniprot ID<br>(b)</th><th>Gene<br>Names</th><th>Protein<br>Score</th><th>Mw/pI <sup>(c)</sup></th><th>Matched<br>Peptides<br>(d)</th><th>Sequence<br>Coverage<br>(%) <sup>(d)</sup></th><th>Overabundant<br>in Muscle <sup>(c)</sup></th></tr></thead><tbody><tr><td></td><td colspan="9"><em>Sarcoplasmic proteins</em> (The early stage of storage, day 1, corresponding Figure 2A,B)</td></tr><tr><td>1314</td><td>Serum albumin</td><td>P14639</td><td>ALB</td><td>153</td><td>71,139/5.80</td><td>15</td><td>43</td><td>LD</td></tr><tr><td>1404</td><td>Creatine kinase M-type</td><td>W5PJ69</td><td>CKM</td><td>127</td><td>43,213/6.66</td><td>15</td><td>40</td><td>LD</td></tr><tr><td>1410</td><td>Fructose-bisphosphate aldolase</td><td>W5PCA0</td><td>ALDOB</td><td>96</td><td>41,555/7.57</td><td>12</td><td>34</td><td>LD</td></tr><tr><td>1417</td><td>Glyceraldehyde-3-phosphate dehydrogenase</td><td>W5PDG3</td><td>GAPDH</td><td>111</td><td>36,110/8.51</td><td>11</td><td>40</td><td>LD</td></tr><tr><td>1448</td><td>Carbonic anhydrase 3</td><td>W5PUC1</td><td>CA3</td><td>83</td><td>29,726/7.70</td><td>6</td><td>22</td><td>LD</td></tr><tr><td>1469</td><td>Heat shock protein family B</td><td>W5P9U1</td><td>HSPB11</td><td>75</td><td>15,737/5.08</td><td>5</td><td>34</td><td>LD</td></tr><tr><td>1475</td><td>Adenylate kinase isoenzyme 1</td><td>C5IJA8</td><td>AK1</td><td>172</td><td>21,750/8.40</td><td>15</td><td>62</td><td>LD</td></tr><tr><td></td><td colspan="9"><em>Sarcoplasmic proteins</em> (The late stage of storage, day 5, corresponding Figure 2C,D)</td></tr><tr><td>1058</td><td>Pyruvate kinase</td><td>W5QC41</td><td>PKM</td><td>147</td><td>58,551/7.24</td><td>20</td><td>40</td><td>LD</td></tr><tr><td>1122</td><td>Creatine kinase M-type</td><td>W5PJ69</td><td>CKM</td><td>109</td><td>43,213/6.66</td><td>14</td><td>37</td><td>LD</td></tr><tr><td>1128</td><td>Fructose-bisphosphate aldolase</td><td>W5PCA0</td><td>ALDOB</td><td>87</td><td>41,555/7.57</td><td>12</td><td>36</td><td>PM</td></tr><tr><td>1134</td><td>Glyceraldehyde-3-phosphate dehydrogenase</td><td>Q28554</td><td>GAPDH</td><td>95</td><td>36,110/8.51</td><td>13</td><td>39</td><td>LD</td></tr><tr><td>1172</td><td>Adenylate kinase isoenzyme 1</td><td>C5IJA8</td><td>AK1</td><td>147</td><td>21,750/8.40</td><td>14</td><td>64</td><td>LD</td></tr><tr><td></td><td colspan="9"><em>Myofibrillar proteins</em> (The late stage of storage, day 5, corresponding Figure 2E,F)</td></tr><tr><td>2172</td><td>Actin, alpha 1</td><td>W5NYJ1</td><td>ACTA1</td><td>70</td><td>42,338/5.23</td><td>11</td><td>25</td><td>LD</td></tr><tr><td>2242</td><td>Enolase 2</td><td>W5P5C0</td><td>ENO2</td><td>117</td><td>47,382/7.60</td><td>14</td><td>37</td><td>LD</td></tr><tr><td>2323</td><td>Fructose-bisphosphate aldolase</td><td>W5PCA0</td><td>ALDOB</td><td>62</td><td>41,555/7.57</td><td>7</td><td>26</td><td>LD</td></tr><tr><td>2327</td><td>Creatine kinase M-type</td><td>W5PJ69</td><td>CKM</td><td>104</td><td>43,213/6.66</td><td>16</td><td>32</td><td>LD</td></tr><tr><td>2330</td><td>Myosin light chain 1</td><td>A0A0H3V7A0</td><td>MYL1B</td><td>73</td><td>20,950/4.95</td><td>9</td><td>58</td><td>LD</td></tr></tbody></table>

<footnote>(a) The numbered spots in gel image (Figure 2). (b) Protein names and accession numbers were taken from the Uniprot database (http://www.uniprot.org (accessed on 8 October 2021)). (c) Theoretical protein mass (Mw; kDa) and isoelectric pH (pI). (d) Number of peptides that matched the protein sequence and total percentage of sequence coverage. (e) Muscle with greater abundance of the protein significance level was indicated (*p* < 0.05).</footnote>
<figure><img src="image_36.png" /><figcaption>**Figure 2.** Two-dimensional gel images (pH range between 3 and 10 and molecular weight from about 20 to 205 kDa) of sarcoplasmic (A,B), the early of storage, day 1; (C,D), the late stage of storage, day 5) and myofibrillar (E,F), the late stage of storage, day 5) proteome from ovine *M. longissimus lumborum* (LL) and *M. psoas major* (PM) muscles.</figcaption></figure>

### 3.3. Correlation of Differentially Abundant Proteins with Meat Color Attributes

The correlation data between meat color attributes and differentially abundant sarcoplasmic proteins are presented in Table 3. In this research, the proteomic changes of proteins may affect the color traits of ovine muscle. Therefore, the method of Pearson's Spearman's correlation was used to infer the possible correlation between meat color attributes and differentially abundant sarcoplasmic proteins.
<table><caption>**Table 3.** Correlation coefficients (Pearson and Spearman) of differentially abundant proteins in ovine *M. longissimus lumborum* (LL) with color attributes (n = 6).</caption><thead><tr><th>Protein Name</th><th>Trait <sup>(a)</sup></th><th>Pearson's Correlation Coefficient <sup>(b)</sup></th><th>Spearman's Correlation Coefficient <sup>(b)</sup></th></tr></thead><tbody><tr><td rowspan="4">Fructose-bisphosphate aldolase (ALDOA)</td><td>a*</td><td>0.126</td><td>0.086</td></tr><tr><td>C*</td><td>0.179</td><td>0.257</td></tr><tr><td>h°</td><td>-0.061</td><td>0.086</td></tr><tr><td>MRA</td><td>0.145</td><td>0.143</td></tr><tr><td rowspan="4">Gyceraldehyde-3-phosphate dehydrogenase (GAPDH)</td><td>a*</td><td>0.127</td><td>0.314</td></tr><tr><td>C*</td><td>0.072</td><td>0.371</td></tr><tr><td>h°</td><td>-0.139</td><td>-0.143</td></tr><tr><td>MRA</td><td>0.232</td><td>0.486</td></tr><tr><td rowspan="4">Carbonic anhydrase 3 (CA3)</td><td>a*</td><td>0.702</td><td>0.829*</td></tr><tr><td>C*</td><td>0.857*</td><td>0.771</td></tr><tr><td>h°</td><td>-0.543</td><td>-0.771</td></tr><tr><td>MRA</td><td>0.346</td><td>0.200</td></tr><tr><td rowspan="4">Adenylate kinase isoenzyme 1 (AK1)</td><td>a*</td><td>0.950**</td><td>0.771</td></tr><tr><td>C*</td><td>0.876*</td><td>0.829*</td></tr><tr><td>h°</td><td>-0.956**</td><td>-0.714</td></tr><tr><td>MRA</td><td>0.876*</td><td>0.657</td></tr><tr><td rowspan="4">Pyruvate kinase (PKM)</td><td>a*</td><td>0.575</td><td>0.543</td></tr><tr><td>C*</td><td>0.416</td><td>0.543</td></tr><tr><td>h°</td><td>-0.644</td><td>-0.429</td></tr><tr><td>MRA</td><td>0.653</td><td>0.657</td></tr><tr><td rowspan="4">Creatine kinase M-type (CKM)</td><td>a*</td><td>0.955**</td><td>0.886*</td></tr><tr><td>C*</td><td>0.877*</td><td>0.886*</td></tr><tr><td>h°</td><td>-0.957**</td><td>-0.878*</td></tr><tr><td>MRA</td><td>0.889*</td><td>0.771</td></tr></tbody></table>

During post-mortem storage, six proteins were correlated with MRA and instrumental color data in *M. longissimus lumborum* muscle (Table 3). There was a significant correlation between redness (a*-value) and adenylate kinase isoenzyme 1 or creatine kinase M-type (r = 0.950, 0.955; *p* < 0.01). Contrarily, a negative relationship was found between hue angle (h°-value) and adenylate kinase isoenzyme 1 or creatine kinase M-type (r = -0.956, 0.957; *p* < 0.01). Chroma (C*-value) also has been found a positive correlation with carbonic anhydrase 3, adenylate kinase isoenzyme 1 or creatine kinase M-type (r = 0.857, 0.876, 0.877; *p* < 0.05). MRA of LL and PM muscles were positively correlated with adenylate kinase isoenzyme 1 and creatine kinase M-type (r = 0.876, 0.889; *p* < 0.05).

Those proteins (overabundant in color-stable muscles), such as fructose-bisphosphate aldolase [16], glyceraldehyde-3-phosphate dehydrogenase [33], adenylate kinase isoenzyme 1, pyruvate kinase [34], and creatine kinase M-type [16,29], exhibited positive correlations with several meat color parameters (a*-value, b*-value, and C*-value), associating them with stability of meat color. In addition, Joseph et al. [12] observed that some sarcoplasmic proteins (antioxidant) showed positive correlations with a*-values and MRA. Those results were partially in agreement with present study.

### 3.4. Functional Roles of Differentially Abundant Proteins and Their Relevance to Color Stability

#### 3.4.1. Glycolytic Enzymes

Four sarcoplasmic proteins (glyceraldehyde-3-phosphate dehydrogenase, fructose-bisphosphate aldolase, pyruvate kinase, enolase 2), which belong to glycolytic enzymes, were over-abundant in LL group (color-stable muscle, Table 2).

Fructose-bisphosphate aldolase catalyzes the reaction that splits the fructose 1,6-bisphosphate into glyceraldehyde 3-phosphate, and dihydroxyacetone phosphate [35]. During the early stage of storage, fructose-bisphosphate aldolase was over-abundant in LL muscle. Nevertheless, an opposite result was observed during late stage of storage (Table 2). Fructose-bisphosphate aldolase was closely related to fast-twitch fiber displaying a higher

<footnote>(a) C* means chroma; h° means hue angle. (b) * means significant at *p* < 0.05, ** means significant at *p* < 0.01.</footnote>
glycolysis metabolic activity [2]. The down-regulated trend of fructose-bisphosphate aldolase involved in glycolytic metabolism indicated that glycolytic activity was highly active in PM muscle which exhibited lower redness (a*-value) and color stability [34].

The glycolysis pathway continues with the conversion of glyceraldehyde 3-phosphate into dihydroxyacetone phosphate, the reaction of which is catalyzed by glyceraldehyde 3-phosphate dehydrogenase. The reaction is NAD-dependent, resulting in the production of pyruvate and NADH [34,35]. Pyruvate kinase catalyzes conversion of phosphoenolpyruvate (PEP) and ADP into ATP and pyruvate, which is one of the main rate-limiting enzymes in glycolysis [35,36]. Enolase (also called muscle-specific enolase) catalyzes the conversion of 2-phosphoglycerate to PEP [37]. The three glycolytic enzymes mentioned above were over-abundant in LL group (Table 2). Glyceraldehyde 3-phosphate dehydrogenase and pyruvate kinase exhibited a positive correlation with MRA (Table 3). Some researchers have discussed the differences between color-stable and color-labile muscles in the literature. Several glycolytic enzymes (β-enolase, pyruvate kinase M2, and glyceraldehyde-3-phosphate dehydrogenase), which were over abundant in *Longissimus lumborum* (color-stable muscle), were found to be positively correlated with a*-values [12,34]. The over-abundance of the enzymes mentioned above could result in higher metabolic activity of glycolytic pathway [34]. This stimulated biological process could accelerate the NADH and pyruvate production. Furthermore, pyruvate is one of the mitochondrial substrates that promote the regeneration of NADH [38]. Besides, NADH supplied by glycolysis is necessary for the NADH–cytochrome b~5~ reductase system of mitochondria [39]. The NADH–cytochrome b~5~ reductase system was deemed as an important electron transport-mediated reduction system for metmyoglobin, whose mechanism has been expounded by some researchers [7]. Researchers [30] considered NADH a critical component of MRA, whose metabolic function plays a key role in reducing metmyoglobin accumulation [19]. In fact, a direct relationship was found between color stability and MRA, in that muscles with enzyme activity (MRA) in higher had greater color stability [7]. Consequently, NADH may play a key role in MetMb enzymatic reduction, which further influences color stability of meat [2].

In this study, Enolase 2 was overabundant in the LL group (color-stable muscle, Table 2). Enolase, one of the key enzymes in glycolysis, catalyzes the formation of phosphoenolpyruvate from 2-phosphoglycerate, leading to an increased rate of glycolysis [40]. Previous literature findings showed that enolase 1 was related to meat color development and quality variation [23,41,42]. In some proteome research, the profile of muscles of different color stability was compared, revealing the correlation between enolase and attributes of meat color [12,32]. The findings in proteomic analysis indicated that enolase was a key enzyme associated with the color of meat.

### 3.4.2. Energy Metabolism Enzymes

Adenylate kinase isoenzyme 1 was overabundant and positively correlated with MRA and chroma (C*-value) in the color-stable LL muscle (Tables 2 and 3). Adenylate kinase (also known as myokinase), a reversible enzyme, plays a role in adenine synthesis and energy metabolism [43]. It provides a way to catalyze the conversion of ADP to ATP and convert AMP to ADP in the cytoplasm [44]. There were different findings within adenylate kinase in the previous associated muscle proteomic studies. Canto et al. [34] noted that the exact mechanism of how adenylate kinase affects color stability remains unclear. Wu et al. [2] found a positive relationship between this enzyme and MRA. It was considered to be a possible potential predictor for meat color stability of *M. longissimus lumborum*. Thus, the exact function mechanism of adenylate kinase in proteomics of muscle color stability is waiting for further investigation.

Creatine kinase M-type was over-abundant in LL muscle (Table 2), demonstrating positive relationships ($r = 0.702, 0.955; p < 0.05$) with redness (a*-value) and MRA (Table 3). Creatine kinase M-type was critical in maintenance of ATP-ADP level, reversibly catalyzing the mutual transformation of phosphate between ATP and phosphocreatine [16]. Creatine kinase was degraded during the early post-mortem stages due to the depletion of ATP [42].
Findings from previous literature showed that the denaturation of creatine kinase and several other related proteins were considered to be related to the paler colored (lower *a*\*-value) *longissimus lumborum* muscle [45]. The present study was in agreement with a previous study, which found that *longissimus lumborum* muscle that creatine kinase M-type was positively correlated (p < 0.05) with *a*\*-values and MRA [16]. Creatine can be used as a natural selective antioxidant because of its ability to scavenge oxygen free radicals [46,47].

Creatine kinase M-type, which is relatively overabundant, can increase creatine levels and minimize myoglobin oxidation, which leads to a lower formation of metmyoglobin and enhanced color stability in color-stable muscles [34].

The functions of carbonic anhydrase 3 were considered as an oxyradical scavenger. In skeletal muscle, this enzyme prevents cells for oxidative damage [48]. In agreement with our results, Yu et al. found that carbonic anhydrase 3 was overabundant in *M. longissimus lumborum* (color-stable muscle) compared with *M. psoas major* (color-labile muscle) [11]. Moreover, we also found this enzyme was correlated with redness (*a*\*-value) and hue angle (h°-value) (r =0.702, -0.543; *p* < 0.05). In practice, an increased glycolysis and oxidative stress can be indicated by carbonic anhydrase 3 in skeletal muscle [49]. In color-stable muscles, the overabundance of this enzyme can increase glycolysis levels and minimize myoglobin oxidation, possibly resulting in lower redness (*a*\*-value) and stabilized color attributes [5].

### 3.4.3. Chaperone Proteins and Structural Proteins

One of the physiological functions of chaperone protein is preventing protein from denaturation [23]. Myoglobin and meat color stability can be compromised by protein denaturation in the process of transformation from muscle to meat [12]. The major chaperone proteins and structural proteins identified were heat shock protein family B, actin (alpha 1), myosin light chain 1, serum albumin, and carbonic anhydrase 3.

The heat shock protein family B (HSP beta-1, also known as HSP27) belongs to chaperone proteins. HSP27 is involved in small heat shock protein family which can prevent structure damage or degradation of proteins in muscle cells and promotes the survival of cell under physiological stress [50,51]. In pig muscle, lower abundance of HSP27 was related to the development of lighter color and PSE zones in semimembranosus [23,52]. Higher abundance of HSP27 was found in beef which was lower in *a*\* and *L*\* values compared with their counterparts. [24]

The physiological function of myosin and actin are mainly related to contractile properties of muscle [43]. During muscle contraction, the function of alpha actin is associated with ATP binding. [53]. Actin alpha has been identified as related to the color of meat in some previous proteomic analyses [4,54]. In a similar study, alpha actin and myosin light chain 1, identified in porcine longissimus muscle, were related to *L*\* value. In the study of Gagaoua et al. [15], actin alpha was correlated with all color coordinates. This protein and several structural proteins were deemed as biomarkers for meat color and of great interest in discriminating between beef color classes. Thus, the functional role of structural proteins in meat color proteomics was important.

### 3.5. *Bioinformatics Analysis of Differentially Abundant Proteins*

The analytical method of bioinformatics has been applied to the related research area of meat color investigation to explore the biological processes and molecular functions of proteins [2,41,55]. Bioinformatics analysis of differentially abundant proteins between muscle types (color-stable and color-labile) provides an efficient tool for deeper understanding of the molecular mechanisms of meat discoloration, based on the proteomic results. In this study, Protein-protein interaction (PPI), gene ontology (GO), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were applied to characterize the functional information of the identified proteins (DAVID Bioinformatics Resources 6.7).
Protein-Protein Interactions

The protein–protein interactions were illustrated in Figure 3 and assessed by String 10.0. There are nodes and edges networks, representing the identified proteins and functional annotation of protein–protein interactions with different colors. The local clustering coefficient and PPI enrichment *p*-value were 0.49 and 1.81 × 10<sup>-8</sup> respectively.

<figure><img src="image_37.png" /><figcaption>**Figure 3.** Protein–protein interaction (PPI) network of differentially expressed proteins from ovine *M. longissimus lumborum* (**A**) and *M. psoas major* (**B**) muscles.</figcaption></figure>