Transcriptomic Study on Variations in Tenderness of Tibetan Lamb Meat

Tibetan sheep are an important part of the pastoral livestock on the Tibetan Plateau, significantly influencing the ecosystem and the lifestyle of the local population. These sheep are adapted to the harsh conditions of the plateau, such as extreme cold and low oxygen levels, thanks to their unique physiological and behavioral adaptations. For instance, Tibetan sheep can maintain nitrogen balance even with low energy and protein intake, indicating their high adaptability to diets low in energy and protein (Zhou et al., 2019). Moreover, Tibetan sheep exhibit good growth performance and meat quality characteristics in response to diets with different protein levels, which is crucial for enhancing the production and consumption of Tibetan sheep meat in the Tibetan Plateau region (Wang et al., 2020). However, the tenderness of Tibetan sheep meat is poorer than other cattle, which has become a great challenge to promote Tibetan sheep meat products. Tibetan sheep meat is favored by consumers for its unique flavor and nutritional value. The rearing management and nutritional supplementation strategies for Tibetan sheep under extremely cold conditions significantly affect their meat quality and growth performance. These strategies enhance the absorptive capacity of the rumen by increasing the abundance of ruminal microbes and improving the development of the ruminal epithelium (Jing et al., 2018). Additionally, the feeding mode of Tibetan sheep has a significant impact on their blood mineral status and appropriate mineral supplementation can improve the nutritional status of Tibetan sheep in winter, thereby affecting meat quality (Xin et al., 2011). The significance of Tibetan sheep on the Tibetan Plateau is not only reflected in their adaptability to harsh environments but also in the advantages of their meat quality, which is important for the development of the local economy and consumer dietary preferences. The intricacies of meat tenderness, particularly in Tibetan sheep, are influenced by a multitude of factors ranging from genetic predispositions (Pelmuş et al., 2020) to environmental interactions (Bai et al., 2022) and post-mortem processing conditions (Devatkal et al., 2022). Recent advances in transcriptomic analyses have provided a novel lens through which the molecular underpinnings of meat tenderness can be explored. For instance, a comparative transcriptome analysis by Wen et al.  (2022) unveiled the dynamic changes in meat quality attributes of the Longissimus Thoracis muscle in Tibetan sheep across different growth stages. This study identified key genes and signaling pathways, such as the AMPK signaling pathway, that are instrumental in regulating muscle fiber transformation and intramuscular fat content, both of which are critical determinants of meat tenderness. These findings underscore the potential of transcriptomic approaches in unraveling the complex biological networks that govern meat quality traits. Furthermore, the interplay between muscle biochemistry, including protein degradation, sarcomere length and collagen characteristics, significantly contributes to variations in meat tenderness. Previous research has elucidated the role of proteolytic enzymes and their inhibitors, such as calpains and calpastatin, in modulating post-mortem muscle proteolysis, thereby affecting the tenderness of meat (Warner et al., 2010). Additionally, the structural and metabolic attributes of muscle fibers, influenced by genetic factors such as myostatin and calpain gene polymorphisms, have been shown to have a considerable impact on meat tenderness.

Given the multifaceted nature of meat tenderness and the sophisticated regulatory mechanisms involved , this paper aims to delve deeper into the transcriptomic landscape of Tibetan sheep muscle tissues. By leveraging high-throughput RNA sequencing technologies, we seek to identify novel gene expression profiles and molecular pathways that are pivotal in determining the tenderness of Tibetan sheep meat. This endeavor not only contributes to the scientific understanding of meat science but also holds practical implications for the sheep industry, aiming to enhance meat quality through genetic and management strategies.

Tenderness is the main evaluation index of meat quality and the main factor affecting consumers’ purchase (Bhat and Pathak, 2012). It is pertinent to start with the Warner-Bratzler Shear Force (WBSF) as a critical indicator for assessing meat tenderness. The WBSF test, widely used in the field of meat science, measures the force required to cut through a meat sample, serving as an objective evaluation of meat tenderness. Generally, lower WBSF values indicate more tender meat, while higher values suggest tougher meat. Meat tenderness significantly impacts consumer acceptance of meat products, making the improvement of tenderness a vital goal within the meat production and processing industry.

In this study, WBSF is employed as the objective standard to determine the tenderness of Tibetan sheep meat, categorizing the experimental animals into high-tenderness and low-tenderness groups based on the WBSF values. This classification provides a foundation for a deeper understanding of the molecular basis underlying the muscle characteristics of Tibetan sheep at varying levels of tenderness.

The purpose of this research is to conduct a transcriptomic analysis of the muscle tissues from both high-tenderness and low-tenderness groups of Tibetan sheep, aiming to uncover the molecular mechanisms behind the variations in meat tenderness. We seek to identify gene expression patterns and regulatory networks closely associated with meat tenderness, in search of potential molecular markers for breeding purposes. These molecular markers not only aid in elucidating the genetic regulatory mechanisms of meat tenderness but also offer theoretical and practical guidance for improving the quality of Tibetan sheep meat through molecular breeding techniques.

Experimental animals and sample collection
 
The experimental subjects in this study comprised 40 healthy, disease-free male Tibetan sheep, with an average body weight of 55.45±4.12 kg. These experimental animals came from Haibei Tibetan Autonomous Prefecture, China’s Plateau Ecological Animal Husbandry Science and Technology Demonstration Park. These animals were allowed to graze freely until they reached the age of 2 months before being slaughtered. During the 12-24 hours preceding slaughter, only water was provided to the animals. The animal slaughter experiment was carried out in Xiahua Food Company (located in Haiyan County, Qinghai Province, China). Muscle tissues were harvested immediately post-slaughter and were promptly submerged in liquid nitrogen, followed by storage at -80^oC until the time for RNA extraction arrived. Additionally, a portion of the muscle samples was reserved for the determination of shear force (Warner-Bratzler Shear Force, WBSF). Shear force analysis was carried out in the sample analysis laboratory of Xihua Company. All experimental procedures were conducted in strict accordance with the ethical guidelines set forth by the Ethics Committee of Qinghai University. The study period began in April 2023 and ended in December 2023.
 
Determination of shear force and experimental grouping
 
To measure the Warner-Bratzler Shear Force (WBSF) (Akanno et al., 2014) and group the experimental animals based on tenderness, the standardized procedure involves cooking one-inch thick steaks from the longissimus lumborum muscle to an internal temperature of 71^oC, then cooling and coring the steaks to obtain samples. These samples are sheared using a Warner-Bratzler shear device at a crosshead speed of 20 cm/min. The force required to shear through the muscle fibers is recorded, with lower forces indicating higher tenderness. Based on the shear force values, the animals were categorized into high-tenderness (low shear force: 45.51 N d≤ WBSF d≤ 53.05 N; n = 15) and low-tenderness groups (high shear force: 101.08 N d≤ WBSF d≤ 109.52 N; n = 15).
 
RNA extraction, library construction and sequencing
 
For RNA extraction, library construction and sequencing, six samples (three from each high-tenderness and low-tenderness group) were selected. Total RNA was isolated using TRIzol Reagent ((Invitrogen), followed by integrity and purity assessments using agarose gel electrophoresis and NanoPhotometer spectrophotometer and Agilent 2100. Libraries were constructed through fragmentation, cDNA synthesis, end repair, adapter ligation using TruseqTM RNA sample preparation Kit (Illumina) and size selection for 250-300 bp cDNA fragments using AMPure XP beads. PCR amplification and a second AMPure XP bead purification were performed to finalize the library. Sequencing was conducted on the HiSeq2000 platform (Illumina).
 
Quality control and reference gene mapping
 
Following sequencing, raw reads were subjected to quality control using FastQC software (version 0.11.8 ) to filter out low-quality reads, reads containing ploy-N.and adapter sequences, yielding clean reads. These clean reads were then aligned to the Tibetan sheep reference genome version 3.1 using HISAT2 software (version 2.1.0) (Kim et al., 2015) . Transcriptome assembly, gene annotation and prediction of novel transcripts were conducted with StringTie software (Pertea et al., 2016)(version 1.3.6). Gene expression levels were quantified using Cufflinks software, utilizing FPKM (fragment mapped per million fragment per thousand base exon) values for normalization.
 
Identification of differentially expressed genes
 
After calculating gene expression levels, differential gene expression analysis was performed using DESeq2 software. The criteria for identifying differentially expressed genes were set at a fold change of e≥2 and a P-value < 0.05.
 
Enrichment analysis of DEGs
 
Following the identification of differentially expressed genes, Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were conducted to understand the biological functions and pathways associated with these genes.
 
qPCR validation
 
For qRT-PCR validation, six differentially expressed genes were selected for confirmation. Bactin was used as an internal reference gene and the Ct (2^-DDCt) method was employed for data analysis. The statistical significance between the high-tenderness and low-tenderness groups was assessed using a t-test in R language. Primers used for qPCR in this study are listed in Table 1.

Transcriptome sequencing quality analysis
 
In the pursuit of elucidating the molecular distinctions between high and low tenderness muscle tissues in Tibetan sheep, transcriptome sequencing was employed to generate comprehensive statistics. In this study, we obtained a total of 45.361 Gb of clean data (sequencing data after quality control). The average clean data per sample was 7.560 Gb. As shown in Table 2, the sequencing output yielded an abundance of raw reads, with subsequent quality filtration processes ensuring a high yield of clean reads. The table reflects a detailed account of the sequencing effort, which is bifurcated into high tenderness (H) and low tenderness (L) groups, each consisting of three replicates. Quality metrics of the sequencing data are notably high, with Q30 values exceeding 89% across all samples, indicating a high accuracy level in the base callings. The GC content hovers around 52%, suggesting a genomic representation without AT or GC bias. The total mapped reads, constituting over 88% on average, reflect the alignment’s efficiency to the reference genome, providing a reliable foundation for downstream differential gene expression analysis. The unique mapping rate stands out, with an excess of 84%, underscoring the specificity of the sequence reads to distinct genomic locations. This specificity is critical for the accurate quantification of gene expression levels and the identification of differentially expressed genes, which are pivotal in understanding the molecular underpinnings of muscle tenderness. The minimal multi-mapped read percentages, averaging around 3.8%, further validate the sequencing data’s precision. These transcriptome sequencing statistics form the bedrock for comparative analyses between high and low tenderness muscle tissues in Tibetan sheep. Such comparative analyses are essential for discovering genetic markers associated with meat quality, thus aiding in the selection and breeding of sheep with desired traits. The data encapsulated in the table are indispensable for subsequent bioinformatic analyses and biological interpretations that may lead to actionable insights into the genetic factors that regulate muscle tenderness in Tibetan sheep.


 
Identification of differentially expressed genes
 
In this experiment, 824 differentially expressed genes were identified across six samples, comprising 241 up-regulated and 583 down-regulated genes (Fig 1). The identified genes were designated as the differential gene set, on which clustering analysis was conducted. The clustering results across different samples are presented in Fig 2.
 




GO and GSEA enrichment analysis of differentially expressed genes
 
The GO enrichment analysis revealed a significant overrepresentation of terms across three primary categories: Biological Process (BP), Cellular Component (CC) and Molecular Function (MF), suggesting a diverse but interconnected range of biological implications (Fig 3). The BP category was particularly enriched for terms related to ‘cellular response to stimulus’, ‘regulation of biological process’ and ‘signal transduction’, indicative of the dynamic and responsive nature of the cellular functions involved. Within the CC category, enrichment was observed for ‘cell part’, ‘organelle’ and ‘membrane’, highlighting the structural aspects of cellular architecture that contribute to the biological roles under investigation. The MF category showed significant representation of terms such as ‘binding’, ‘catalytic activity’ and ‘transporter activity’, underscoring the functional capabilities of the gene products analyzed. These enriched GO terms point towards a comprehensive set of biological themes that may be pivotal in the context of the studied biological system or condition. The scatterplot visualization utilized in the analysis provides a succinct representation of the enriched terms, allowing for an immediate grasp of the most prominent functional groups. The terms are clustered based on semantic similarity, reinforcing the interconnectedness of the biological processes, cellular components and molecular functions represented. The use of semantic similarity measures ensures that closely related terms are plotted proximally, thus facilitating a clearer interpretation of the overarching biological themes.



Gene Set Enrichment Analysis (GSEA) was utilized to ascertain the predominant KEGG pathways that were significantly enriched in our dataset (Fig 4). The analysis revealed differential pathway activation profiles, as depicted in the enrichment plots. The left panel illustrates the top six pathways with the most significant enrichment scores. These pathways include ‘Cysteine and methionine metabolism’, ‘Fat digestion and absorption’, ‘Glycine, serine and threonine metabolism’, ‘Maturity onset diabetes of the young’, ‘Pancreatic secretion’ and ‘Valine, leucine and isoleucine degradation’, which are indicative of the metabolic shifts associated with the condition under investigation. The enrichment scores, represented by the peak heights of the running sum, suggest a concerted biological response or adaptation. The right panel presents the subsequent six pathways, which include ‘Acute myeloid leukemia’, ‘Basal cell carcinoma’, ‘Cell cycle’, ‘Galactose metabolism’, ‘IL-17 signaling pathway’ and ‘p53 signaling pathway’. These pathways underscore the perturbations in cellular processes such as proliferation, differentiation and stress response mechanisms. Each plot portrays a running-sum statistic across a ranked list of genes, where the presence of genes in the indicated pathway contributes positively and the absence negatively, to the cumulative score. The peak of the running enrichment score (ES) for each pathway reflects the degree to which the genes are over-represented at the top or bottom of the ranked list. The ranked list metric is a product of the expression change direction and the statistical significance, thus ensuring that genes with small p-values and consistent regulation patterns contribute more heavily to the ES.



Validation of differentially expressed genes by qRT-PCR
 
The eight candidate genes screened in this study were subjected to qRT-PCR verification of HS / LS muscle tissue (Fig 5). The qPCR results were consistent with the sequencing results, thus verifying the accuracy of the sequencing results. Primers are shown in Table 1.



The Tibetan sheep, native to the high-altitude Tibetan Plateau, thrives in harsh environments, enduring extreme cold, winds and low oxygen levels. prized for its high-quality wool, long, dense and naturally crimped, ideal for luxurious textiles. Additionally, they provide lean, flavorful meat and nutrient-rich milk for local communities. Well-suited to extensive grazing, Tibetan sheep are raised in semi-nomadic or nomadic traditions, integral to Tibetan livelihoods and cultural heritage. However, the differential molecular mechanisms underlying the tenderness of Tibetan sheep meat remain to be fully elucidated to date.

Transcriptome sequencing of muscle tissues from Tibetan sheep with different tenderness levels revealed a series of differentially expressed genes (DEGs), including key genes related to muscle development, growth and metabolism, such as MyoD, PAK6, PDGFRA, Myf5, MRF4, DCN, IGF, GH, FGF and GDF8. Firstly, we can initiate the discussion from the impact of fat deposition on Tibetan sheep meat tenderness. Adipose tissue, as an endocrine organ, plays a crucial role in regulating lipid metabolism (Wheeler et al., 1994). Nishimura (2010) observed histologically that connective tissue was damaged and mechanical strength decreased during intramuscular fat deposition, resulting in improved Tibetan sheep meat tenderness Additionally, gene families like growth hormone (GH) (Liu et al., 1992; Ma et al., 2010), growth differentiation factor 8 (GDF8) (Fu et al., 2019; Li et al., 2013; Sakuma et al., 2000) and insulin-like growth factor (IGF) (An et al., 2014) are involved in regulating muscle cell growth, differentiation and metabolism. In our study, we observed some genes closely associated with muscle development and growth, such as MyoD, Myf5 and MRF4 (Ren et al., 2014). These genes belong to the myogenic regulatory factor (MRFs) family, which influences muscle development by regulating the proliferation, differentiation and fusion of muscle cells (Ren et al., 2014). Furthermore, our research also identified genes related to cellular signaling pathways, such as PDGFRA and FGF. These genes may affect muscle cell proliferation and metabolism by activating specific signaling pathways like the PI3K/Akt pathway.

Subsequently, our focus shifted to the impact of energy metabolism on meat tenderness. Carbohydrates primarily serve the physiological function of providing energy for vital activities (Jian et al., 2021). The oxidative breakdown of sugars involves three main pathways: anaerobic glycolysis, aerobic oxidation and the pentose phosphate pathway (Granot and Kelly, 2019). Typically, postmortem glycogen in lamb is fermented into lactate, while ATP is also degraded, yielding phosphate. Accumulation of lactate and phosphate in muscle tissue increases H+ concentration, lowering pH, releasing calcium ions, catalyzing fiber degradation and enhancing meat tenderness (Obrenovitch et al., 1988). This aligns with our findings, indicating that energy metabolism pathways may be a critical factor influencing meat tenderness across different tenderness groups. Our research has identified numerous differentially expressed genes related to energy metabolism, including MyoD, Myf5, MRF4, among others. These genes are involved in the regulation of muscle development, growth and metabolism. Specifically, MyoD (Nicolas et al., 1996; Tan et al., 2014), Myf5 (Indriulyte and Miceikiene, 2010; Maak et al., 2006) and MRF4 belong to the myogenic regulatory factor family (MRFs family), which influences muscle development by controlling the proliferation, differentiation and fusion of muscle cells (Odle et al., 2020; Zhang et al., 2018). MyoD is a muscle-specific transcription factor crucial for the proliferation, differentiation and fusion of muscle cells (Nicolas et al., 1996; Tan et al., 2014). Surprisingly, these genes have been reported in other mammals; for instance, research has found that the Myf5 gene can influence the meat quality traits of cattle (Yang and Chen, 2010) and pigs (Indriulyte and Miceikiene, 2010), suggesting its candidacy as a molecular marker for meat quality traits (Indriulyte and Miceikiene, 2010; Yang and Chen, 2010; Zhang et al., 2014). Genes such as IGF (An et al., 2014), GH (Liu et al., 1992; Ma et al., 2010) and GDF8 (Fu et al., 2019; Li et al., 2013; Sakuma et al., 2000) also participate in the growth and differentiation of muscle cells, playing crucial roles in muscle development. Indeed, MSTN (Myostatin) serves as a suppressor of muscle growth, with mutations linked to elevated muscle mass traits observed across various animal species, including other sheep breeds (Boman et al., 2010), cattle (Luo et al., 2014), etc. Therefore, in Tibetan sheep muscle tissue, the expression level of MSTN may influence muscle growth and development, consequently affecting meat tenderness. In summary, our research further supports the influence of energy metabolism on meat tenderness and reveals a series of differentially expressed genes related to muscle development, growth and metabolism, providing a new theoretical basis for improving meat quality.

Furthermore, we investigated the impact of muscle structural proteins on meat tenderness. Muscle structural proteins, particularly myofibrillar proteins, have been demonstrated to be closely associated with meat tenderness (Li et al., 2013a). Our research results indicate differential expression of certain genes such as PAK2, PAK6, PDGFRA, among others, across high and low tenderness groups, potentially reflecting the importance of muscle structural proteins in different meat tenderness qualities. Specifically, the PAK6 gene exhibits higher expression in the high tenderness group compared to the low tenderness group, possibly associated with its facilitation of muscle actin fiber skeleton degradation. Moreover, the relevant studies have found that the upregulation of the PAK6 gene may facilitate the degradation of the actin cytoskeleton, thereby affecting muscle structure and tenderness (Ribeiro et al., 2014). PAK6 is a protein kinase closely associated with processes like cytoskeletal remodeling and cell motility. In muscle tissue, elevated PAK6 expression may enhance muscle cell activity, affecting muscle structure and function and consequently, meat tenderness. PDGFRA encodes a receptor tyrosine kinase involved in regulating cell proliferation, growth and development (Gao et al., 2018). Lastly, Fibroblast Growth Factor (FGF) participates in regulating cell growth, differentiation and muscle repair processes (Chen and Robert, 2009). In other animals muscle tissue, FGF may influence muscle tissue structure and meat tenderness by promoting cell proliferation and repair. These results align with previous studies on yak and goat muscle tissues, indicating potential commonalities in meat tenderness regulatory mechanisms across different species. For instance, key genes in yak muscle tissue play crucial roles in meat tenderness formation and in our study, similar genes in Tibetan sheep muscle tissue exhibit analogous expression patterns, providing important clues for further exploration of the molecular mechanisms underlying meat quality formation.

According to the results of Gene Ontology (GO) enrichment analysis, significant enrichment phenomena were observed in three major categories: biological processes, cellular components and molecular functions. Enrichment in biological processes primarily involves cellular responses to stimuli, regulation of biological processes and signal transduction, indicating the dynamic nature and responsiveness of the studied cellular functions. Enrichment in cellular components includes cellular parts, organelles and cell membranes, emphasizing attention to cellular structural aspects. In terms of molecular functions, enrichment encompasses binding, catalytic activity and transporter proteins, highlighting the functional capabilities of gene products. Based on the results of gene set enrichment analysis, several significantly enriched pathways were identified, involving metabolic shifts and condition-related biological responses. These pathways include pathways related to asparagine and methionine metabolism, fat digestion and absorption, glycine, serine and threonine metabolism, maturity onset diabetes of the young, pancreatic secretion and valine, leucine and isoleucine degradation, among others. These enriched pathways reflect metabolic changes associated with the studied conditions. Additionally, pathways involving cellular processes such as proliferation, differentiation and stress response mechanisms were observed, such as acute myeloid leukemia, basal cell carcinoma, cell cycle, galactose metabolism, IL-17 signaling pathway and p53 signaling pathway. Through comparison with results from other studies on yak and goat, we can better understand the enrichment of KEGG pathways observed in our sequencing data. Consistent findings can further confirm the reliability of our results and enhance our understanding of the regulation of biological processes and pathways.

Through transcriptomic analysis, we elucidated the molecular basis of meat tenderness in Tibetan sheep, identifying 824 differentially expressed genes, including 241 upregulated and 583 downregulated genes. Genes such as MyoD, PAK6, PDGFRA, Myf5, IGF, GH, FGF and GDF8, involved in crucial processes like muscle development, growth and metabolism, serve as candidate molecular markers. GO and KEGG enrichment analysis revealed significant enrichment of these genes in biological processes, cellular components and molecular functions. These findings provide valuable insights into enhancing the quality of Tibetan sheep meat through gene regulation.

This work was supported by China Agriculture Research System of MOF and MARA (CARS-39-35), the Open Project of State Key Laboratory of Plateau Ecology and Agriculture, Qinghai University (2023-ZZ-11) and Qinghai Provincial Middle-aged and Young Scientific and Technological Talents Support Project (2023QHSKXRCTJ17).

Zhina Du and Yang Xiang designed the experiments and drafted the manuscript. Sequencing data processing and visualization was done by Wu Sun and Xiayang Jin. Shike Ma and Lizhuang Hao carried out quantitative experiments.

Informed consent has been obtained from all individuals included in this study.

The datasets generated during and analyzed during the current study are available from the corresponding author on reasonable request.

Authors state no conflict of interest.