A Proteomic Study for the Discovery of Beef Tenderness Biomarkers and Prediction of Warner–Bratzler Shear Force Measured on Longissimus thoracis Muscles of Young Limousin-Sired Bulls

Beef tenderness is of central importance in determining consumers’ overall liking. To better understand the underlying mechanisms of tenderness and be able to predict it, this study aimed to apply a proteomics approach on the Longissimus thoracis (LT) muscle of young Limousin-sired bulls to identify candidate protein biomarkers. A total of 34 proteins showed differential abundance between the tender and tough groups. These proteins belong to biological pathways related to muscle structure, energy metabolism, heat shock proteins, response to oxidative stress, and apoptosis. Twenty-three putative protein biomarkers or their isoforms had previously been identified as beef tenderness biomarkers, while eleven were novel. Using regression analysis to predict shear force values, MYOZ3 (Myozenin 3), BIN1 (Bridging Integrator-1), and OGN (Mimecan) were the major proteins retained in the regression model, together explaining 79% of the variability. The results of this study confirmed the existing knowledge but also offered new insights enriching the previous biomarkers of tenderness proposed for Longissimus muscle.


Introduction
Meat-eating quality consists of a complex set of sensory traits including tenderness, flavour, and juiciness, each of which plays an important role in defining the appeal of beef to consumers [1,2]. Amongst these quality attributes, however, tenderness is considered to be one of the most important factors in purchase decisions regarding beef, with negative experience on toughness contributing to a lower likelihood of repeat purchase [3]. To meet the expectations of consumers, beef producers must pursue the provision of consistent highquality beef. The underlying mechanisms involved in dictating the final meat tenderness are intricate, with muscle biochemistry interacting with processing, influenced by several factors including breed [4,5], gender [6], age at slaughter [7], muscle type [8,9], cooking temperature [5], stress at slaughter [10], and post-slaughter management and many other factors from farm-to-fork [9,11].
There have been a number of studies using omics tools to, firstly, enhance our understanding of the pathways and processes contributing to beef tenderness variation [12,13] and secondly, to propose prediction equations to explain the observed variability in this The hydrolysis of the protein bands was carried out with 48 µL of a 25 mM ammonium bicarbonate buffer-12.5 ng/µL trypsin solution (Promega) per band for 5 h in an oven at 37 • C. Then, 30 µL buffer was added periodically during hydrolysis so that the bands were always covered with liquid. The extraction of the peptides was carried out under ultrasound (15 min) with acetonitrile and trifluoroacetic acid. Then, the supernatant was transferred into 500 µL Eppendorf tubes and dry concentrated using a Speedvac for 2 h. The volume was adjusted exactly to 20 µL with a solution of isotopologic peptides (50 pmol/µL) that was diluted 18 times in a 0.05% Trifluoroacetic acid (TFA) solution. After passing through the ultrasonic bath (10 min), the entire supernatant was transferred to the High performance liquid chromatography (HPLC) vial before LC-MS/MS analysis.
For the separation, the hydrolysate was injected into the nano-LC-MS/MS (Thermo Fisher Scientific) using an Ultimate 3000 system coupled to a QExactive HF-X mass spectrometer (MS) with a nanoelectrospray ion source. Briefly, 1 µL of hydrolysate was first preconcentrated and desalted at a flow rate of 30 µL/min on a C18 pre-column 5 cm length × 100 µm (Acclaim PepMap 100 C18, 5 µm, 100 Å nanoViper) equilibrated with trifluoroacetic acid 0.05% in water to remove contaminants that could potentially disrupt the efficiency of the mass spectrometry analysis. After 6 min, the concentration column was put in line with a nano debit analytical column operating at 400 nL/min. The peptides were then separated according to their hydrophobicity (column C18, length 25 cm, diameter 75 µm, SN 10711310), using a gradient of a solution of acetonitrile (ACN/FA-99.9/0.1) of 4 to 25% in 50 min.

LC-MS/MS Data Processing and Protein Identification
The raw files from the LC-MS/MS were aligned against the Bos taurus database (i.e., ref_bos_taurus, 23,970 sequences) with Mascot V.2.5.1 (http://www.matrixscience.com, accessed on 30 August 2020). The precursor and fragment mass tolerance were set up at 10 ppm and 0.02 Da, respectively. The variable modifications included carbamidomethy-Foods 2021, 10, 952 4 of 20 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).

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.

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.proteinside.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].

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.

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). 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).  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 3. Best regression equation of WBSF based on the list of the significant differential proteins from Table 2. 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). The R-square of the correlation is given.

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

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.

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

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

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.

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 myozenin 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 alphaactinin 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.

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

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.

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

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 H 2 O 2 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].

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

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.