A comprehensive study on the longissius dorsi muscle of Ashdan yaks under different feeding regimes based on transcriptomic and metabolomic analyses

Abstract

Yak is an important dominant livestock species at high altitude, and the growth performance of yak has obvious differences under different feeding methods. This experiment was conducted to compare the effects of different feeding practices on growth performance and meat quality of yaks through combined transcriptomic and metabolomic analyses. In terms of yak growth performance, compared with traditional grazing, in-house feeding can significantly improve the average daily weight gain, carcass weight and net meat weight of yaks; in terms of yak meat quality, in-house feeding can effectively improve the quality of yak meat. A combined transcriptomic and metabolomic analysis revealed 31 co-enriched pathways, among which arginine metabolism, proline metabolism and glycerophospholipid metabolism may be involved in the development of the longissimus dorsi muscle of yak and the regulation of meat quality-related traits. The experimental results increased our understanding of yak meat quality and provided data materials for subsequent deep excavation of the mechanism of yak meat quality.

Keywords:

• Yak
• longissimus dorsi muscle
• meat quality
• transcriptomics
• metabolomics

Introduction

The yak is an important large mammal living in highland areas,¹ which has exceeded 20 million worldwide, and is mainly distributed in high-altitude areas between 3000 and 5000 m.² More than 95% of the world’s yaks are mainly in the alpine and low-oxygen areas cantered on the Qinghai-Tibetan Plateau in China,³ which is an important dominant livestock species in plateau animal husbandry.⁴ The distribution area of yaks is characterized by low temperature, large temperature difference between day and night, strong radiation and low oxygen partial pressure, etc. Under the influence of long-term special ecological environment and natural selection, yaks have the advantages of high resistance to disease, strong resistance to adversity and strong vitality, etc.,⁵ and the unique genetic mechanism of yaks has been formed in the process of natural evolution. Due to the yak living in the cold region of plateau for a long time, coupled with the lack of technology of grassroots personnel and the backward form of feeding management, these reasons limit the growth and development of the yak.⁶ As a result, the economic value of yaks cannot be properly realized, and although the relevant grass-roots research departments have taken corresponding countermeasures, the situation of ‘relying on the sky to feed the animals’ has not been fundamentally resolved.

The ‘Ashdan Yak’ is a new breed of yak bred by the Lanzhou Institute of Animal Husbandry and Veterinary Medicine of the Chinese Academy of Agricultural Sciences (CAAS) and the Datong Breeding Farm of Qinghai Province after more than 20 years of efforts.⁷ The Ashdan yak is characterized by high meat yield, high reproductive capacity and mild temperament.⁸ Compared with other horned yaks, Ashdan yaks are easier to be managed and reared by herdsmen breeders. As the first hornless yak breed in the world, the Ashdan yak is conducive to the improvement of the ecological environment as well as the economic development of Qinghai and other highland areas,⁹ and has been well promoted with a wide range of applications. Therefore, it is of great significance to study the complex quantitative traits and important economic traits of the Ashdan yak.¹⁰

Transcriptomics is a discipline that reveals complex biological pathways and molecular mechanisms through the study of gene expression and regulation,^11,12 which allows for the mining of candidate genes related to specific functions,¹³ and is able to reflect the molecular response to specific physiological processes at the transcriptional level.^14,15 The transcriptome has spatial and temporal properties,¹⁶ which allows for a more comprehensive interpretation of gene expression in specific tissues,¹⁷ which in turn leads to a deeper excavation of the biological functions of the relevant genes.^18,19 Metabolomics is the study of metabolites in tissues and changes in metabolites due to different factors (genetic, environmental or disease, etc.).^20,21 Because metabolites are the end products of the genome,²² the corresponding changes in metabolite species and content can form a more direct relationship with changes in animal phenotype, so metabolomics is also the closest to phenotypic changes.^23,24 Untargeted metabolomics allows for a comprehensive analysis of endogenous small molecule metabolites in organisms^25,26 and is also more widely used than targeted metabolomics. A single transcriptomics analysis of results acting on multiple substances cannot directly determine the extent of changes in metabolites that lead to changes in metabolites, whereas metabolomics can unearth differences in metabolite species in the metabolic pathway, which is important for a comprehensive analysis of biosynthesis and regulatory mechanisms. Carrillo et al.,²⁷ using transcriptome and metabolome, identified biological pathways associated with differences in grass-fed and grain-fed Angus cattle. Grass-fed cattle produce lower total fat and less blood stress. This will provide greater insight into the merits and mechanisms of muscle development. Li et al.,²⁸ applied multi-omics technology to study the longest dorsal muscle tissues of Guangling donkeys with different meat tenderness and screened differentially expressed genes and differential metabolites. The joint analysis showed that the significantly differentially expressed genes and metabolites played key regulatory roles in the metabolic pathway, and could be used as potential candidate genes and metabolites for meat tenderness, which could provide a theoretical basis for the investigation of the molecular regulatory mechanism of meat tenderness and molecular breeding. Zhan et al.,²⁹ conducted comprehensive transcriptomics and metabolomics analyses of intramuscular fat (IMF) and meat colour traits in Enshi black pigs, revealing multiple pathways, such as adipogenesis and catabolism, and screened several functional genes, emphasizing nicotinamide metabolism and glutathione metabolism as important trait-related pathways. The results provide a basis for further research on the biological mechanism of IMF deposition. Liu et al.³⁰ analysed the effects of CRTC3 overexpression on the metabolomics and transcriptome profiles of porcine subcutaneous white adipocytes through a multi-omics association, revealing the metabolic effects of CRTC3 on porcine adipocytes.CRTC3 overexpression altered the expression of genes and pathways related to adipogenesis, fatty acid metabolism and glycerol phospholipid metabolism in vitro. Sun et al.³¹ utilized transcriptomic and metabolomic approaches to analyse the transcriptomic data and fatty acid composition of four vital metabolic tissues (rumen, liver, muscle and dorsal fat) in cattle, exploring the molecular regulatory mechanisms within these tissues that promote the formation of fatty acids in cattle. This research holds significant promise for future cattle breeding initiatives. Xi et al.³² studied the differences in the transcriptome-level fatty acid metabolism pathway of Bamei pork. The results showed that the metabolism levels of DX3 and DX4 were significantly higher than those of other fatty acids. Additionally, the study found that the differences in the G1 and G2 genes played an important regulatory role in fatty acid metabolism, which affects the quality of pork. These findings are of great significance for understanding the species’ traits and gene functions.

Multi-omics analysis can better illustrate the cellular life process from different dimensions, obtain the dynamic changes in the growth and development of the animal body, and then deeply study the regulatory mechanisms related to complex traits.³³ Yak as an important genetic resource pool, analysing the genetic mechanisms of economically important traits using multi-omics is an important task in this field of research.^34,35 In this experiment, we selected the longissimus dorsi muscle of 12 Ashan yaks from the Datong Breeding Farm in Qinghai Province for transcriptomic and metabolomic studies to screen out differentially expressed genes and metabolites, to reveal the biogenetic mechanism of yak meat quality, and to provide a scientific basis for the selection of yak meat for breeding and efficient production.

Materials and methods

Ethics statement

All animal studies conducted in this study were approved by the Animal Administration and Ethics Committee of Lanzhou Institute of Husbandry and Pharmaceutical Sciences, CAAS (SYXK-2014–0002).

Experimental design and sample collection

In this experiment, 12 male Ashdan yaks (4 years old) were selected from the Datong Breeding Yak Farm in Qinghai Province, and were divided into fattening group (Group F) and grazing group (Group Y) according to the feeding method in the farm, of which 6 were in each group. The experimental period was from April 5 2019 to October 2 2019, and the experimental group was reared in full housing, with a preparatory trial of 10 d and a formal trial of 170 d. Group F was prepared with total mixed ration (TMR) feed in strict accordance with the Beef Cattle Feeding Standard (NY/T 815-2004), as shown in Table 1. During the preparatory test, all the yaks were dewormed internally and externally, and the experimental group was provided with basic food once a day in the morning and once in the evening, and water was freely available for drinking, and the food intake of yaks in the experimental group was kept and recorded at any time. The control group (grazing group) was treated by grazing alone in natural pasture without anthropogenic supplementation.

Table 1. Composition of the basal diet in the experimental group.

Item	Content (%)
Corn	48.53
Bran (esp. of wheat)	6.22
Rapeseed cake	7.25
Baking soda	0.25
Table salt	0.25
4% Premixed compound	1.50
Oat hay	22.00
Alfalfa hay	14.00
Total	100.00

All test yaks were fasted and dehydrated for 24 h before sampling and slaughtering, and yak slaughtering was carried out in accordance with the National Standard of the People’s Republic of China, ‘Cattle Slaughtering Operation Procedures’ (GB/T 19477-2004). Carcass and live weights of pre-slaughter yaks were recorded and eye muscle area and slaughter rate were measured. Meat colour, shear, cooking loss and pH of the longissimus dorsi muscle between the 12th and 13th pairs of ribs of yaks were measured at 45 min and 24 h. Finally, the longissimus dorsi muscle samples were cut into 1–3 cm³ pieces, transported back to the laboratory using dry ice, and stored in a −80 °C refrigerator for long-term storage.

RNA-Seq analysis

Total RNA was extracted from yak longissimus dorsi muscle tissue using Trizol reagent, and the whole process was carried out in an ultra-clean bench. Total RNA concentration and purity were examined using an Agilent 2100 Bioanalyzer, and finally total RNA integrity was examined by 1% agarose gel electrophoresis. The mRNA was then enriched with magnetic beads with Oligo (dT), and the library was constructed using TruSeq stranded total RNA with Ribo-Zero Gold. The constructed libraries were qualified with an Agilent 2100 Bioanalyzer for QC and then sequenced using the Illumina HiSeq^TM 2500 Sequencer to produce 125 or 150 bp bipartite data. After passing the quality control (QC), sequencing was performed using an Illumina sequencer.

The raw data generated in high-throughput sequencing are fastq format sequences. In order to get high quality data for subsequent analysis, quality filtering of raw reads is required. The data were first quality controlled, joints removed, and low quality bases and N bases filtered using Trimmomatic³⁶ software to obtain high quality clean reads. The use of Hisat³⁷ was compared to the BosGru_v2.0 reference genome. Gene FPKM³⁸expression was quantified using Cufflinks³⁹ software, HTseq-Count⁴⁰ software was applied to obtain the number of reads falling on the gene in each sample, and the data were normalized using the estimateSizeFactors function that comes with DESeq (2012) R package⁴¹ software. Differential genes with p values less than 0.05 and multiplicity of difference greater than 2 were selected for GO and KEGG functional enrichment analysis. The reads were assembled by StringTie⁴² software to splice new transcripts, which should be more than 200 bp away from the existing annotation gene and not less than 180 bp in length. Variable shear events were detected using ASprofile⁴³software. SNP and INDEL analyses were performed using Samtools⁴⁴ and Bcftools⁴⁵ software.

LC/MS untargeted metabolomic analysis

The metabolome samples were first pretreated by weighing 30 mg of the samples, adding 20μL of each of the internal standards (L-2-chlorophenylalanine, 0.3 mg/mL; Lyso PC17:0, 0.01 mg/mL, both configured from methanol), followed by the addition 400 μL of aqueous formaldehyde solution (CH3OH:H2O V:V = 4:1) for metabolite extraction. The processed samples were transferred to LC injection vials and stored at −80 °C in a refrigerator until LC/MS analysis. QC samples were prepared by mixing equal volumes of extracts from all samples, ensuring that the QC samples had the same volume as the rest of the samples. A liquid-mass spectrometry (MS) system consisting of AB ExionLC ultra-high performance liquid chromatography and AB Triple TOF 6600 high-resolution mass spectrometer was used for metabolomics analysis in this experiment. Chromatographic conditions included ACQUITY UPLC BEH C18 (100 mm × 2.1 mm, 1.7 um), column temperature (45 °C), mobile phase (A-water, B-acetonitrile), flow rate (0.4 mL/min) and injection volume (5 μL). MS conditions included ion source (ESI), and the samples’ MS signals were acquired using positive and negative ions respectively.

Joint transcriptomics and metabolomics analysis

This is because differential genes obtained from transcriptomics analysis alone will correspond to substances in multiple metabolic pathways, affecting the direct determination of the amount of metabolites accumulated. Metabolomics analysis can better identify differences in the types and ion abundance of metabolites in the pathways. Therefore, the joint analysis of transcriptomics and metabolomics on the same sample will be helpful to dig out the biological mechanisms related to yak meat quality traits and screen out the related genes and metabolites. Based on the transcriptomic and metabolomic data, expression correlation analysis was performed to construct KGML networks, etc.

Statistical analysis of data

Results such as production data were analysed using the t-test in GraphPad Prism 8 (GraphPad Software, La Jolla, CA) and expressed as mean standard deviation. Differences were judged to be significant when p < 0.05.

Results

Measurement of growth performance, slaughtering performance and meat quality

The experiment was conducted to determine the growth performance and slaughter performance of Ashdan yaks within groups Y and F, respectively (Figure 1(A,B)). It was found (Table 2) that final body weight, average daily weight gain, carcass weight, net meat weight, slaughter rate and eye muscle area of Ashdan yaks in group F were significantly increased (p < 0.01) compared to group Y. The final body weight, carcass weight and net meat weight of the yaks in group F were 377.17, 208.04 and 176.23 kg, respectively, which were 70.34, 93.76 and 81.5 kg higher than that of the yaks in group Y, respectively. It indicates that the growth performance indexes of Ashdan yaks can be effectively improved by housed feeding in this experiment. Meat quality of the longissimus dorsi muscle of Ashdan yaks was determined at the time of sample collection and the results are shown in Table 3. The differences in pH (45mIn, 24h), cooking loss and shear were highly significant (p < 0.01) in yaks of group F compared to group Y. The differences in pH (45 min, 24 h), cooking loss and shear were highly significant (p < 0.01).

Figure 1. Radar plots of growth performance and longissimus dorsi muscle quality of Ashdan yaks in group Y and group F. (A) Slaughtering performance of Ashdan yak, (B) Determination of meat quality of Ashdan yak.

Table 2. Determination of growth performance and slaughtering performance of Ashdan yaks.

Item	Group Y	Group F
Starting weight/kg	208.33 ± 12.53	215.67 ± 21.83
Final weight/kg	306.83 ± 24.34^B	377.17 ± 15.77^A
Average daily weight gain/kg	0.52 ± 0.12^B	0.86 ± 0.15^A
ketone weight/kg	124.28 ± 8.77^B	208.04 ± 8.75^A
net meat weight/kg	94.73 ± 4.55^B	176.23 ± 10.61^A
bone weight/kg	31.08 ± 1.96^b	27.17 ± 2.46^a
slaughtering rate/%	42.49 ± 10.24^B	55.26 ± 14.25^A
eye muscle area/cm^2	57.89 ± 8.83^A	76.64 ± 10.45^B

¹Different upper case letters in the same row of the table indicate highly significant differences i.e., (p < 0.01) (same for the following table).

Table 3. Determination of the longissimus dorsi muscle quality of Ashdan yaks.

Item	Group Y	Group F
L*45 min	6.32 ± 0.80	6.98 ± 1.27
a*45 min	28.12 ± 2.32	29.29 ± 2.49
b*45 min	5.67 ± 1.04	6.42 ± 1.38
pH*45 min	6.89 ± 0.32^A	6.34 ± 0.24^B
L*24 h	7.61 ± 0.99	8.86 ± 1.20
a*24 h	32.68 ± 0.91	33.23 ± 1.19
b*24 h	9.94 ± 1.18	8.62 ± 0.81
pH*24 h	5.49 ± 0.25^B	6.03 ± 0.28^A
Steaming losses/%	17.82 ± 3.84^B	27.21 ± 1.70^A
water damage from dripping/%	15.06 ± 2.30	17.27 ± 3.12
shearing force/N	17.70 ± 1.51^A	14.69 ± 2.35^B

Transcriptome results and analysis

Since one yak died during the trial in group Y, we performed reference transcriptome sequencing on samples from a total of 11 yaks from both groups. A total of 77.79 G of CleanData was obtained, and the effective data volume of each sample was distributed in the range of 6.79–7.76 G, the Q30 bases were distributed in the range of 91.65–92.65%, and the average GC content was 51.82%; and the comparison rate with the reference genome was 86.63–95.05% (Schedule 1). We applied the FPKM method to the gene expression and the results are shown in Figure 2(A–C) (Schedule 2). As seen in Figure 2(A), the results of the test grouping were more reasonable, with more homogeneous results within groups F and Y, and differences between the two groups, with the mean FPKM values of the genes in groups F and Y being 17.89 and 16.85, respectively. The FPKM densities of the genes in the two groups were more evenly distributed (Figure 2(B)), and intervalization of the FPKM revealed that the highest number of genes expressed in the samples was located at 0.5 < FPKM < 1 and the lowest at 1 < FPKM < 10 (Figure 2(C)). We tested the correlation between the samples, and the results were all greater than 0.8, indicating that the biological replication was good enough to continue with subsequent analyses (Figure 2(D)). Figure 2(E,F) shows a good distribution of samples between the two groups with more dispersed samples between groups and relatively clustered samples within groups. The differentially expressed genes were analysed and a total of 1265 differentially expressed genes were obtained from the two groups, including 566 up-regulated genes and 699 down-regulated genes (Figure 3(A)). Hierarchical clustering results showed significant differences in the expression profiles of the two groups of genes (Figure 3(B,C), Schedule 3). GO enrichment analysis was performed to identify significantly enriched entries in differentially expressed genes, and at the transcriptome level, the differentially expressed genes compared between the two groups were significantly enriched in 1489 GO terms (p < 0.05). The main distributions were in the regulation of steroid metabolism, collagen fibre organization, denervation response, muscle adaptation regulation, negative regulation of phosphorylation of pathway-restricted SMAD proteins, positive regulation of myotube differentiation, positive regulation of the mitotic cell cycle, cell migration and regulation of myofibrils (Figure 3(D), Schedule 4). KEGG enrichment analysis was performed used to identify significantly enriched entries in differentially expressed genes, and at the transcriptome level, the differential genes compared between the two groups were significantly enriched in 34 KEGG terms (p < 0.05). The main distribution was in the FoxO signalling pathway, HIF-1 signalling pathway, glutathione metabolism and tyrosine metabolism (Figure 3(E), Schedule 5).

Figure 2. Expression level analysis of protein-coding genes of Ashtan yak in group Y and group F. (A) Box plots of gene expression levels in the two groups. Box plots for each region for five statistics (maximum, upper quartile, median, lower quartile and minimum, from top to bottom), (B) Density distribution plot of FPKM values for each sample. The horizontal coordinates of the points on the curve indicate the logarithmic value of the corresponding sample FPKM, and the vertical coordinates of the points indicate the probability density, (C) Distribution of gene expression of samples, (D) Plot of correlation coefficients between samples based on gene expression was obtained, (E) Sample 3D distribution, (F) Obtaining sample clustering based on gene expression. The horizontal coordinate indicates the sample name, the vertical coordinate indicates the corresponding sample name, and the colour represents the magnitude of the correlation coefficient.

Figure 3. Number of differentially expressed mRNAs in Ashdan yaks in groups Y and F. (A) Histogram of differential gene expression statistics. The horizontal axis is the individual comparison groups; the vertical axis is the number of differentiated genes in the comparison groups, where up is the number of up-regulated genes that are significantly different and down is the number of down-regulated genes that are significantly different, (B) Volcano map of differentially expressed genes. Gary is non-significantly different genes, red and green are significantly different genes, (C) Cluster analysis of group differences in subgroups. The red colour in the figure indicates highly expressed genes and the blue colour indicates lowly expressed protein-coding genes, (D) GO enrichment analysis, (E) KEGG enrichment analysis.

Metabolome results and analysis

The test was prepared by mixing F1-F5 into F6 for metabolomics analysis. The stability and reproducibility of the analytical methods were evaluated using QC. Ion peaks with RSD > 0.4 in the QC group were removed. Unsupervised PCA identification model and QC metabolite intensity distribution were applied to assess inter-sample variability and QC variability, and the samples in the QC group were more clustered, indicating that there were differences in metabolite species and quantities between the two groups, and that the analytical method was stable and reproducible (Figure 4(A,B)). A total of 2752 substance peaks were extracted by total ion chromatography (924 in the negative ion mode and 1828 in the positive ion mode); 1500 metabolites were obtained by combined metabolite annotation using primary MS and secondary MS (468 in the negative ion mode and 1032 in the positive ion mode) (Figure 4(C)). Multivariate statistical analysis using unsupervised PCA and supervised OPLS-DA was seen to be effectively divided into two categories, and it was also found that within-group samples were closer and clustered together, with better data reproducibility; between-group samples were farther apart, suggesting that there were differences in metabolite phenotypes between the two groups (Figure 4(D,E)). The OPLS-DA model was built to prevent overfitting, and it was found that all Q2 points were lower than the rightmost original Q2 point and the regression line of Q2 was < 0 at the intersection of the vertical coordinates, and the intercepts of R2 and Q2 were 0.925 and −0.265, respectively, which indicated that the model was more stable and did not suffer from overfitting (Figure 4(F), Schedule 6).

Figure 4. Principal component analysis and orthogonal partial least squares discriminant analysis of metabolites. (A) Principal component analysis of QC samples, (B) QC metabolite intensity distribution assessment, (C) Qualitative quantitative results. Substance peaks and metabolite statistics graphs, (D,E) Unsupervised PCA and supervised OPLS-DA for multivariate statistical analysis, (F) Permutation diagram. Loads can range from −1 to 1. Loads close to −1 or 1 indicate that the variable has a very strong effect on the component. A load close to 0 indicates that the variable has a very weak effect on the component.

A total of 1500 metabolites were obtained in this study, and the differential metabolites were screened by the VIP value (VIP > 1) of the OLPS-DA model and the p value of the t-test (p < 0.05). From the results, it can be seen that 105 differential metabolites were initially identified, of which 56 were up-regulated and 49 were down-regulated, accounting for 53.33% and 46.67% of the total number of differential metabolites, respectively (Figure 5(A), Schedule 7). The 105 differential metabolites screened were subjected to cluster analysis (Figure 5(B)), from which two groups can be clustered into two categories, suggesting that these clustered metabolites may have similar functional annotations or be in the same metabolic pathway, and that the differences in expression between the samples are obvious. The results of correlation analysis of the differential metabolites are shown in Figure 5(C). KEGG enrichment analysis at the metabolome level revealed that a total of 32 pathways were significantly enriched (p < 0.05). The main distribution was in lipid metabolism, glycosylphosphatidylinositol (GPI), pyrimidine metabolism and oxidative phosphorylation (Figure 5(D,E), Schedule 8).

Figure 5. Differential metabolite analysis. (A) Volcano diagram for VIP and P screening. The red dots represent differential metabolites that were significantly up-regulated in the experimental group, the blue dots represent differential metabolites that were significantly down-regulated, and the grey dots represent differential metabolites that were not significant, (B) Heat map of differential metabolites, (C) Correlation analysis plot of TOP-50 differential metabolites, (D) TOP-20 metabolic pathway enrichment map, (E) TOP-20 bubble chart.

Results of combined transcriptome and metabolome analysis

In order to understand the meat quality as well as the deeper regulatory mechanisms in Ashdan yak, the transcriptome and metabolome were comprehensively analysed by expression correlation analysis, pathway mapping, KEGG enrichment analysis and KGML network map in this study. Correlation analysis of the differentially expressed genes and metabolites of the two groups showed that the correlation between differentially expressed genes and metabolites was high, which indicated that the changes in differentially metabolites regulating the longissimus dorsi muscle of the Ashdan yak were related to the changes in differentially expressed genes (Figure 6(A,B)). Differential gene and differential metabolite association network maps were plotted according to p value < 0.05 (Figure 6(C)). Joint KEGG pathway coanalysis was performed, resulting in 31 co-enriched pathways, among which arginine metabolism, proline metabolism and glycerophospholipid metabolism may be involved in the regulation of the development of the longissimus dorsi muscle and meat quality traits in Ashdan yaks. The KGML network analysis was then plotted (Figure 6(D)) (Schedule 9).

Figure 6. Joint transcriptomic and metabolomic correlation analysis. (A) Correlation matrix of gene expression and metabolites, (B) Correlation heat map, (C) Relevance network diagram, (D) KGML network diagram.

Discussion

As one of the important local livestock breeds in China, yak plays an important role in people’s diet, life and economic development. Yaks mainly feed on natural pasture and live in a pollution-free environment, so yak products are considered a natural green food.⁴⁶ Meanwhile, in order to adapt to the high altitude and low oxygen environment, the yak’s body will make corresponding adjustments, such as the metabolism of various types of amino acids, oxidative phosphorylation, the expression of genes related to glycolysis and the increase of myoglobin. In this experiment, we investigated the biological mechanisms regulating yak meat quality and longissimus dorsi muscle development through different feeding methods and in combination with transcriptome and metabolome.^47,48 In this experiment, compared with group Y, there were significant differences in final body weight, average daily gain and net meat weight of yaks in group F, which increased by 70.34, 0.34 and 81.5 kg, respectively, which indicated that housed feeding could significantly improve yak performance. Shear force is a direct indicator of meat tenderness,⁴⁹ and is one of the most important factors influencing consumers’ choice of meat, which directly affects the edible and commercial value of meat. Meat tenderness is in turn affected by a variety of factors, including the content of intramuscular connective tissue, IMF content and muscle fibre diameter.⁵⁰ Yak meat is generally tougher and chewy because of its thicker muscle fibres. In this experiment there was a significant difference in the shear force of the housed F group compared to the Y group, which indicated that the shear force could be effectively reduced and the tenderness of the meat could be improved by housed feeding, and it has also been shown in numerous studies that housed feeding can effectively reduce the shear force value.^51–53 The pH in muscle is also one of the most important indicators of meat quality, with corresponding effects on tenderness and flavour, among others.⁵⁴ There is a strong positive correlation between pH and tenderness, the higher the pH the higher the tenderness.⁵⁵ The pH decreases gradually due to the production of lactic acid from glycolysis and the catabolism of ATP to produce phosphoric acid in the muscle after slaughter,⁵⁶ which is also consistent with the result that the pH at 24 h was lower than that at 45 min in this experiment. At the same time, this reduction in pH puts the cells in an acidic environment, where the ionic balance is disrupted, promoting the release of Ca²⁺ which in turn activates calpain, which then degrades fibrous proteins, resulting in an increase in the tenderness of the meat.⁵⁷ However, the decrease in pH also causes autolysis of μ-calcium-activating enzymes and a decrease in enzyme activity, which in turn leads to a decrease in the rate of myofibrillar protein hydrolysis.⁵⁸ This suggests that changes in pH will have an effect on the catalytic activity of endogenous enzymes in the organism, which in turn will have an effect on the degradation rate of myofibrillar proteins, and ultimately affect the tenderness of meat. In conclusion, reasonable housing can improve the production performance of yaks and yak meat quality.

In the transcriptome analysis, a total of 1265 differential genes (566 up-regulated and 699 down-regulated genes) were screened in both groups and were significantly enriched (p < 0.05) in 1489 GO entries, mainly in the collagen fibre organization (ADAMTS2, COL1A1, COL1A2, COL5A2, DPT, and P3H4), positive regulation of myotube differentiation (MAML1, MAMSTR, MYOG and SMYD1), myofibrillar ganglia (ARHGEF25, CALM3, FHOD3 and HABP4), cytoplasmic (CRYAB, ACTA1, HSPB1, SEPW1 and APOBEC2), and denervation and muscle adaptive regulation (FBXO32, HDAC4, MYOG and TRIM63), among others. A total of 34 KEGG entries were significantly enriched (p < 0.05), mainly in the FOXO signalling pathway (CCNB3, FOXO6, INSR, PIK3R1, PIK3CD, TGFB2, CCNG2, FBXO32 and FOXO3), tyrosine metabolism (FAH, AOX1, PNMT, TH and DDC), oxytocin signalling pathway (CACNG7, CACNA2D2, MYLK3, ROCK2 and MYLK4) and glutathione metabolism (MGST3, GPX4, GPX7 and GPX1). This suggests that the regulation of yak meat quality is not determined by a single factor, but is co-regulated by multiple factors. Actin alpha 1 (ACTA1) is a class of heterotrimeric forms of actin, which constitutes the main component of skeletal muscle contractile proteins, and is involved in the development of fibres and the movement of cells and organelles in skeletal muscle.⁵⁹ Sibut et al.⁶⁰ found that the expression of ACTA1 in the longissimus dorsi muscle of low-moon-old Korean cattle was significantly higher than that of high-moon-old Korean cattle, suggesting that this gene may be involved in the regulation of beef quality traits in Korean cattle. Venturini et al.,⁶¹ in their study of ACTA1 and meat production traits in broiler chickens found that ACTA1 is involved in the production performance of broiler chickens and can be used as an important molecular marker to increase economically important traits in broiler chickens. Shin and Chung⁶²found that ACTA1 was differentially expressed in different grades of marbled beef, with low-grade marbled beef being significantly higher than high-grade marbled beef, which may indicate that ACTA1 is involved in the regulation of meat quality traits. The results of numerous studies have shown that ACTA1 is a major candidate gene that affects the quality of livestock and poultry meat. F-box protein 32 (FBXO32) is an E3 ubiquitin-protein ligase that plays an important role in the ubiquitin-proteinase pathway and has been shown to be one of the major pathways for degrading myosin.⁶³ Wang et al.⁶⁴ found that FBXO32 plays an important role in regulating muscle development and muscle proteolysis. Chen et al. found that FBXO32 was significantly associated with carcass and meat quality traits, but not with growth traits.⁶⁵ However, it has also been shown that this gene may be a candidate locus for muscle and fat deposition, and its mutation may affect the production performance of livestock.⁶⁶

In the metabolome level analysis, 105 differential metabolites were identified (56 up-regulated and 49 down-regulated). One of them, L-Carnitine (L-Carnitine), is a vitamin-like nutrient involved in fatty acid metabolism in the body and promotes beta-oxidation of fatty acids. Ying et al. found⁶⁷ that the addition of L-carnitine to the feed of fattening pigs significantly reduced the shear force of their longest dorsal muscles and increased their tenderness. The expression of L-carnitine in group F was significantly higher than that in group Y in this experiment. Also in the study of James et al.,⁶⁸ it was found that L-carnitine could improve meat quality by changing the meat colour. However, the specific regulatory mechanism of L-carnitine on yak meat quality needs to be further investigated and deeply explored. Carnosine is a natural water-soluble dipeptide composed of β-alanine and L-histidine, which is involved in animal growth, has been the oxidation of lipids and improves the meat quality of livestock and poultry. Ma et al. found⁶⁹ that the colour redness value of the longissimus dorsi muscle of the pig’s back increased and the drip loss decreased after proper supplementation of carnosine, indicating that the quality of pork was effectively improved. The same results were also reflected in the study of Ma et al. The addition of carnosine significantly reduced drip loss, improved muscle tethering force and enhanced muscle quality and tenderness in muscles stored for 72 and 96 h.⁷⁰ Muscle tethering force measures the water retention capacity of meat, which is important for both flavour and tenderness.⁷¹ D'Astous-Pagé⁷² found that pork with high expression of myostatin had better water holding capacity and meat colour value. The results of numerous studies have shown that myostatin can improve meat quality by reducing drip loss and increasing meat colour value.⁷³ By discovering the transient state of gene expression through transcriptome data, and then observing the final result of gene expression through transcriptome data, joint analysis can be performed to study the intrinsic formation mechanism of traits from different stages of transcription to metabolism. A total of 31 co-enriched channels were screened by co-analysing the functional enrichment of differential genes and differential metabolites, and a high correlation between differential genes and differential metabolites was also found. The regulation of the longissimus dorsi muscle of the yak is not determined by a single factor, but by a combination of multiple factors. The results revealed that the differential gene FBXO32 and the differential metabolite L-carnitine showed a positive and strong correlation, suggesting that they are involved together in the regulation of the longissimus dorsi muscle in yaks. Similarly FBXO32 and DL-acetylcarnitine also showed strong positive correlation. However, the current research on the key regulatory genes for yak meat quality is not deep enough, especially the correlation between the key genes and traits of yak meat quality still needs to be further explored, so as to improve the yak meat industry and yak meat quality.

Conclusion

In this study, 1265 differential genes and 105 differential metabolites were obtained from the two histological studies by combining the transcriptome and metabolome to compare the longissimus dorsi muscle characteristics of grazing and housed yaks, respectively. Genes (ACTA1 and FBXO32, etc.) and metabolites (L-carnitine and myostatin, etc.) potentially involved in yak meat quality regulation were screened. And these related genes were found to have significant correlations with metabolites, which may be jointly involved in the regulatory role of yak meat quality. Yak meat is rich in nutrients, unique in flavour, and rich in a variety of essential amino acids. However, due to the special environment of yaks and the fact that most yaks are raised in grazing mode, the tenderness of yak meat is poorer than that of other cattle breeds, and screening for related regulatory genes or metabolites to improve the quality of yak meat is particularly critical. Meat tenderness, as an important indicator of meat quality, needs to be further investigated so as to reveal the molecular mechanisms regulating yak meat tenderness and provide a theoretical basis for improving yak meat quality marker-assisted breeding.

Authors’ contribution

Conceptualization, T.W., X.M. and P.Y.; methodology, T.W.; software, T.W. and X.M.; formal analysis, Q.Z. and C.M.; investigation, Z.Z.; resources, X.W., M.C., and H.P.; data curation, H.P., C.L. and P.Y.; writing – original draft preparation, T.W.; writing – review and editing, T.W., X.M. and C.M.; visualization, Z.Z. and X.G.; supervision, P.Y.; project administration, P.Y.; funding acquisition, P.Y. All authors contributed to the article and approved the submitted version.

Institutional review board Statement

All animal studies conducted in this study were approved by the Animal Administration and Ethics Committee of Lanzhou Institute of Husbandry and Pharmaceutical Sciences, Chinese Academy of Agricultural Sciences (SYXK-2014–0002).

Acknowledgements

The authors would like to thank the Key Laboratory of Yak Breeding Engineering in Gansu Province for providing the experimental conditions and the teachers and brothers for their help.

Disclosure statement

The authors declare no conflict of interest.

Data availability statement

Data used in this study are available from the corresponding author on reasonable request.

Additional information

Funding

This research was funded by the National Natural Science Foundation of China (32102524), the Major Science and Technology Special Projects in Gansu Province (21ZD10NA001, GZGG-2021-1), the National Beef Cattle Industry Technology & System (CARS-37), the Agricultural Science and Technology Innovation Program (25-LZIHPS-01) and the Science and Technology program of Gansu Province (20JR5RA580) the Shenzhen Science and Technology Program (KCXFZ20201221173205015).