Lipidomics and Transcriptome Reveal the Effects of Feeding Systems on Fatty Acids in Yak’s Meat

The differences of fatty acids in yak’s meat under graze feeding (GF) and stall feeding (SF) regimes and the regulation mechanism of the feeding system on the fatty acids content in yak ’s meat was explored in this study. First, the fatty acids in yak’s longissimus dorsi (LD) muscle were detected by gas liquid chromatography (GLC). Compared with GF yaks, the absolute content of ΣSFAs, ΣMUFAs, ΣUFAs, ΣPUFAs and Σn-6PUFAs in SF yak’s LD were higher, whereas Σn-3PUFAs was lower; the relative content of ΣMUFAs, ΣPUFAs, Σn-3PUFAs and ΣUFAs in SF yak’s LD were lower, whereas ΣSFAs was higher. The GF yak’s meat is healthier for consumers. Further, the transcriptomic and lipidomics profiles in yak’s LD were detected by mRNA-Sequencing (mRNA-Seq) and ultra-high performance liquid chromatography-mass spectrometry (UHPLC-MS), respectively. The integrated transcriptomic and lipidomics analysis showed the differences in fatty acids were caused by the metabolism of fatty acids, amino acids, carbohydrates and phospholipids, and were mainly regulated by the FASN, FABP3, PLIN1, SLC16A13, FASD6 and SCD genes in the PPAR signaling pathway. Moreover, the SCD gene was the candidate gene for the high content of ΣMUFA, and FADS6 was the candidate gene for the high content of Σn-3PUFAs and the healthier ratio of Σn-6/Σn-3PUFAs in yak meat. This study provides a guidance to consumers in the choice of yak’s meat, and also established a theoretical basis for improving yak’s meat quality.


Introduction
With living standards improving, consumers are more and more focused on meat quality [1], and request to get more superior meat. Meanwhile, the meat industry focuses on developing valuable nutritional properties in meat. Moreover, meat safety incidents frequently happen in the world, and meat safety becomes a focus too [2]. The content and type of fatty acids in livestock meat are among the key factors which affect the sensory qualities of meat such as tenderness, color, and cooking loss [3,4]; on the other hand, the composition of fatty acids in livestock meat is a critical index in evaluating the meat nutrition. Saturated fatty acids (SFAs) in livestock meat can increase the risk of cardiovascular disease for consumers [5], while more polyunsaturated fatty acids (PUFAs), especially n-3 PUFAs like eicosapentaenoic acid (EPA) and docosahexaenoic acid (DHA), can protect blood vessels and prevent cardiovascular and cerebrovascular diseases [6,7]. The livestock meat with a lower ratio of Σn-6/Σn-3 PUFAs is a benefit to human health [8], so is very popular with consumers.
Researchers in livestock husbandry are attempting to alter the composition of fatty acids in livestock meat by increasing the content of unsaturated fatty acids (UFA), especially n-3 PUFAs in recent years. There are significant differences in digestion and absorption gene functions and lipid metabolism pathways in yak muscle under different feeding system can be revealed by the approach combined lipidomics and transcriptomics.
In this study, we hypothesized that the feeding system could affect the fatty acids metabolism in yak muscle by the differential expression of genes associated with fatty acids, and ultimately cause a change to the content and composition of fatty acids in it. The absolute and relative content of fatty acids in yak's longissimus dorsi (LD) muscle were determined, and the differences in the characteristics of fatty acids in yak muscle between GF and SF were observed. Further, the differentially expressed genes (DEGs) and significantly different lipids (SDLs) in GF and SF yak's LD were identified by RNA-Seq transcriptomics and non-targeted lipidomics, respectively. The DEGs and SDLs were annotated by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genome (KEGG), then the molecular mechanism of the feeding system affecting fatty acids in yak muscle was explored. This study investigated the differences of fatty acids in yak's LD under SF and GF and established a theoretical basis for the evaluation of yak's meat under different feeding system and provided a reference for improving intramuscular fat (IMF) content in yaks.

Animals and Sample Collection
The animal experiment was approved by the Ethics Committee of the Lanzhou Institute of Husbandry and Pharmaceutical Sciences, Chinese Academy of Agricultural Sciences (Permit No. SYXK-2020-0166). A total of twelve healthy male yaks (the same genetic background, two-year old, 210.33 ± 10.23 kg) were selected in May and were divided into two groups in the study. The GF group: six yaks (n = 6) were only grazed in natural pasture with no supplements; the SF group: six yaks (n = 6) were fed with total mixed ration (TMR) food in a stall. The GF and SF yaks were fed in Qinghai province in China, and the yaks in each group can freely eat grass or MTR and drink water by themselves, respectively. All experimental animals were dewormed before the test. By September, all yaks were humanely slaughtered by professional technicians at a commercial abattoir. The slaughter procedure was conducted in accordance with European Commission Regulation. The LD samples from between the 12th and 13th ribs of the left side of the carcass were immediately collected after the yaks were slaughtered. Part of the LD samples were put in liquid nitrogen immediately, and the rest were frozen at −20 • C. The grass sample was collected in September too. The composition of TMR and the content of common nutrition and major fatty acids in grass and TMR are shown in Table 1.

Determination of Fatty Acids
The absolute content of fatty acids in the yak's LD was detected according to the method described in references [42][43][44] with some modifications. First, the lipid in yak's LD was extracted with a mixed solvent of chloroform and methanol (v:v, 2:1) by Soxhlet extraction for three times, then the combined extract solution was dried under blown nitrogen in a water bath. Second, the extracted lipid in the yak's LD was decomposed into free fatty acids with a sodium hydroxide methanol solution by a basic hydrolysis in a water bath. Third, the fatty acids were transferred into the fatty acid methyl esters (FAMEs) by an esterification reaction with boron fluoride-methanol solution. The FAMEs were extracted with normal heptane and were dried under nitrogen. The residue was redissolved with a mobile phase, then was determined on an gas chromatography system (7890A, Agilent Corp., Santa Clara, CA, USA). The analytes were confirmed by comparing their retention times with the values of standard compounds. The absolute content of fatty acids was calculated by the absolute content of FAMEs, and the relative content was calculated as follows: the relative content of fatty acid = (the absolute content of a certain fatty acid)/(the absolute content of total fatty acids) × 100. Table 1. The composition of the total mixed ration (TMR) and the content of common nutrition and major fatty acid in natural grass and TMR (air-dry basis).

RNA Extraction and Sequencing
Three samples of yak's LD were randomly selected from six samples in the SF and GF groups, respectively. According to the protocol of the manufacturer, the total RNA in yak's LD was extracted with the mirVanaTM miRNA Isolation Kit (Ambion Inc., Foster City, CA, USA). The purity of the extracted RNA was determined with a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Santa Clara, CA, USA). Moreover, the RNA integrity was evaluated with the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The mRNA libraries for sequencing were prepared using TruSeq Stranded mRNA LTSample Prep Kit (Illumina, San Diego, CA, USA), then these libraries were sequenced on the Illumina sequencing platform (HiSeqTM 2500) and 125 bp pairedend reads were generated. Raw data containing ploy-N and the low-quality reads were removed, then the clean reads were obtained. The clean reads were mapped to a reference genome using hisat2 and the Q30 and GC-content of the clean reads were calculated. The numbers of reads were normalized against the reads per kilo base of transcripts per million (RPKM) to compute the gene expression levels. The RPKM value of each gene was calculated using cufflinks, and the read counts of each gene were obtained by htseq-count.

Lipids Extraction and MS Data
Fifty mg of the LD sample was put in an eppendorf (EP) tube, and Lyso PC17:0 solution in methanol, internal standard 2-chloro-l-phenylalanine solution in methanol and pure water were successively added to the tube. The LD sample was broken by ultrasonication, then the mixture of solutions and LD sample was transferred into a centrifuge tube. The solution of chloroform-methanol was added into the centrifuge tube, and the mixture was vortexed for 1 min. After centrifuging, the subnatant was transferred into a sample bottle and dried. The residue was dissolved in the solution of isopropanol-methanol and vortexed for 30 s. The solution was transferred into an EP tube and centrifuged. Two hundred µL supernatant was collected and filtered into a LC vial. A Nexera UPLC (Shimadzu, Kyoto, Japan) with Waters ACQUITY UPLC BEH C 18 (100 mm × 2.1 mm, 1.7 µm) was used to separate the extract. The moving phase was composed of (A) acetonitrile and water (v:v, 60:40), containing 0.1% formic acid and 10 mmol/L ammonium formate and (B) acetonitrile and isopropyl alcohol (v:v, 10:90), containing 0.1% formic acid and 10 mmol/L ammonium formate. The column temperature, the flow rate of moving phase and injection volume were 45 • C, 0.3 mL/min and 5 µL, respectively. The Q Exactive MS system (Thermo Scientific™, Waltham, MA, USA) was used to collect the MS data and the parameters of the MS system were as follows: Positive ion mode-heater temp: 300 • C; the flow rate of sheath gas: 45 arb; the flow rate of aux gas: 15 arb; the flow rate of sweep gas: 1 arb; spray voltage: 3.5 KV; capillary temp: 320 • C; S-Lens RF level: 50%; MS1 scan ranges: 120-1800. Negative ion mode-heater temp 300 • C; the flow rate of heath gas: 45 arb; the flow rate of aux gas: 15 arb; the flow rate of sweep gas: 1 arb; spray voltage: 3.1 KV; capillary temp: 320 • C; S-Lens RF level 50%; MS1 scan ranges: 120-1800.

Quantitative Real-Time PCR (qPCR)
Six DEGs, including PLIN2, FABP3, LEP, SCD, ACAA1 and ACOT7, were selected to confirm the mRNA-Seq results. The information of primers is listed in Table S1. The reference gene for qPCR in the yak muscle was the β-actin gene. The qPCR analysis for the target genes was carried out with ABI Prism 7500 instrument (Applied Biosystems, Carlsbad, CA, USA). A volume of 1.5 µg RNA, 0.5 µL gDNA remover, 5 µL 5 × TransScript All-in-one SuperMix, 0.5 µL of each primer and nuclease-free water were used in reverse transcription (RT) reaction. The reactions were carried out on a GeneAmp ® PCR System 9700 thermal cycler (Applied Biosystems, Foster City, CA, USA) for 15 min at 42 • C and 5 s at 85 • C. Next, the mixture for the RT reaction was diluted 10 times with nuclease-free water, then stored at 20 • C. The Real-time PCR was carried out on a LightCycler ® 480 Real-time PCR Instrument (Roche, Basel, Switzerland), and 0.2 µL cDNA, 5 µL 2 × PerfectStartTM Green qPCR SuperMix, 0.2 µL of each primer and 3.6 µL nuclease-free water were used in a 10 µL PCR. The parameters for the Real-time PCR were 94 • C for 30 s, followed by 45 cycles of 94 • C for 5 s and 60 • C for 30 s. Relative gene expression levels were determined using the 2 −∆∆Ct method [45].

Statistical Analysis
The data of the fatty acids was analyzed with the independent-sample T test in the software SPSS 16.0 (SPSS Inc., Chicago, IL, USA), and the results are shown by means ± standard error of the mean (SEM). Correlations among the fatty acids, DEGs and SDLs in yak's LD were conducted by Pearson correlation analysis in SPSS 16.0 too. Significant differences and correlation were considered at p < 0.05 and correlation coefficient > 0.8 or <−0.8. These genes with fold change > 2 or < 0.5 and p < 0.05 were deemed as DEGs. The original Q Exactive mass data in raw format were processed with the software Lipid Search. A data matrix was combined from the positive and negative ion data. The matrix was imported into R to carry out principal component analysis (PCA). Further, orthogonal partial least squares discriminant analysis (OPLS-DA) was used to distinguish the lipids which were different in yak LD between the GF and SF groups. The SDLs were selected with VIP > 1.0 and p < 0.05. The GO analysis for DEGs was carried out with R based by running queries for each DEGs against the GO database. The KEGG analysis for DEGs and SDLs was performed to identify those KEGG pathways with p < 0.05.

Lipidomics Results
The PCA score plots for the samples of GF yak's LD, SF yak's LD and quality control (QC) is shown in Figure 1a. It can be found that the QC group samples were congregated tightly in a small area, which indicated that the lipidomics platform UHPLC-MS possessed high stability and the analysis possessed excellent reliability. Further, the OPLS-DA analysis showed the variables responsible for differentiation in yak's LD between GF and SF ( Figure 1b). As observed herein, R2Y(cum) = (0, 0.638) and Q2(cum) = (0, −0.734) indicated the model was robust and without overfitting (Figure 1c). It was found that there was an obvious difference between GF and SF yak's LD, which demonstrated the feeding system can significantly affect the lipids in yak's LD. Lipids with significant differences between the GF and the SF group are shown in volcano plots (Figure 1d). A total of 75 SDLs in yak's LD under GF and SF were obtained (Table S2). Of these, 26 SDLs were up-regulated in LD of SF yaks, whereas 49 SDLs were down-regulated (Figure 2a). A total of 11 SDLs were enriched in 17 KEGG pathways (Table S3), and the bubble chart of KEGG pathways is shown in Figure 2b. The KEGG pathway being relevant to fat acids metabolism included glycerophospholipid metabolism (bom00564), fat digestion and absorption (bom04975), regulation of lipolysis in adipocytes (bom04923), insulin resistance (bom04931), linoleic acid metabolism (bom00591), alpha-linolenic acid metabolism (bom00592) and arachidonic acid metabolism (bom00590). The above KEGG pathways mainly focus on UFAs metabolism and glycerophospholipid metabolism.

Transcriptome Results
The results of the PCA displaying the distinct biological variation among GF and SF yak's LD are shown in Figure 3a, and the effect of feeding system to the yak's muscle was obvious. The DEGs in LD of GF and SF yaks are shown in Table S4. Compared with LD of GF yaks, there were 1682 DEGs in LD of SF yaks, and 429 DEGs were up-regulated, while 1253 DEGs were down-regulated. The heatmap of the DEGs in GF and SF yak's LD is shown in Figure 3b. After GO enrichment of DEGs, the biological processes involving to fatty acids contained fatty acids metabolism process, glucose metabolic process, regulation of protein processing (Figure 4a). The DEGs were enriched in 45 major KEGG pathways (Table S5), and the KEGG pathways involving to fatty acids included PPAR signaling pathway (bom03320), glycolysis/gluconeogenesis (bom00010), glycine, serine and threonine metabolism (bom00260), carbohydrate digestion and absorption (bom04973), fatty acid biosynthesis (bom00061), protein digestion and absorption (bom04974) (Figure 4b). The enrichment scores of PPAR signaling pathway (3.70) and fatty acid biosynthesis (3.52) were the highest in all KEGG pathways, and contained 19 DEGs and 4 DEGs, respectively. The heatmap of the crucial DEGs (PLIN2, DBI, CPT2, FASN, FABP3, FADS6, SLC16A13, SCD, LEP, ACOT7, ACAA1 and ACSL6) involving the regulation of fatty acids content is shown in Figure 4c. in LD of SF yaks by contrast with GF yaks. Red dots represent the significant differences lipids (SDLs) whose concentration was up-regulated in SF yak LD, while the blue dots represent the SDLs whose concentration was down-regulated and grey dots represent non-SDLs. The volcano plot of total lipids in LD of SF yaks by contrast with GF yaks. Red dots represent the significant differences lipids (SDLs) whose concentration was up-regulated in SF yak LD, while the blue dots represent the SDLs whose concentration was down-regulated and grey dots represent non-SDLs.

Verification of mRNA Sequencing Using qPCR
As shown in Figure 4d, by contrast with GF yaks, the expressions of PLIN2, FABP3, LEP and SCD genes in the LD of SF yaks were up-regulated, whereas the expressions of ACAA1 and ACOT7 genes were down-regulated. All six DEGs in the LD of SF and GF yaks possessed the similar expression patterns in qPCR and mRNA-Seq data, which indicated the reliability of mRNA-Seq data for yak's LD in this study.

Transcriptome Results
The results of the PCA displaying the distinct biological variation among GF and SF yak's LD are shown in Figure 3a, and the effect of feeding system to the yak's muscle was obvious. The DEGs in LD of GF and SF yaks are shown in Table S4. Compared with LD of GF yaks, there were 1682 DEGs in LD of SF yaks, and 429 DEGs were up-regulated, while 1253 DEGs were down-regulated. The heatmap of the DEGs in GF and SF yak's LD is shown in Figure 3b. After GO enrichment of DEGs, the biological processes involving to fatty acids contained fatty acids metabolism process, glucose metabolic process, regulation of protein processing (Figure 4a). The DEGs were enriched in 45 major KEGG pathways (Table S5), and the KEGG pathways involving to fatty acids included PPAR signaling pathway (bom03320), glycolysis/gluconeogenesis (bom00010), glycine, serine and threonine metabolism (bom00260), carbohydrate digestion and absorption (bom04973), fatty acid biosynthesis (bom00061), protein digestion and absorption (bom04974) ( Figure  4b). The enrichment scores of PPAR signaling pathway (3.70) and fatty acid biosynthesis

Joint Analysis of the Lipidomics and Transcriptome
The joint analysis of the transcriptome and lipidomics data for yak's LD under SF and GF was carried out to explore the molecular mechanism for the differences of lipid metabolism in LD of SF and GF yaks. The regulatory processes of crucial DEGs to the lipid metabolism was shown in Figure 5. The biggest differences in lipid metabolism between SF and GF yak's LD relates to the phospholipid metabolism and triglyceride metabolism by the regulation of DEGs in PPAR signaling pathway.

Correlation of Fatty Acids, DEGs and SDLs
Results of Pearson correlation analysis among the fatty acids, crucial DEGs in fatty acids regulation, crucial SDLs in fatty acids metabolism are shown in Figure 6.

Verification of mRNA Sequencing Using qPCR
As shown in Figure 4d, by contrast with GF yaks, the expressions of PLIN2, FABP3, LEP and SCD genes in the LD of SF yaks were up-regulated, whereas the expressions of ACAA1 and ACOT7 genes were down-regulated. All six DEGs in the LD of SF and GF yaks possessed the similar expression patterns in qPCR and mRNA-Seq data, which indicated the reliability of mRNA-Seq data for yak's LD in this study.

Joint Analysis of the Lipidomics and Transcriptome
The joint analysis of the transcriptome and lipidomics data for yak's LD under SF and GF was carried out to explore the molecular mechanism for the differences of lipid metabolism in LD of SF and GF yaks. The regulatory processes of crucial DEGs to the lipid metabolism was shown in Figure 5. The biggest differences in lipid metabolism between SF and GF yak's LD relates to the phospholipid metabolism and triglyceride metabolism by the regulation of DEGs in PPAR signaling pathway.

Correlation of Fatty Acids, DEGs and SDLs
Results of Pearson correlation analysis among the fatty acids, crucial DEGs in fatty acids regulation, crucial SDLs in fatty acids metabolism are shown in Figure 6.

Discussion
The deposition velocity of SFAs and MUFAs is faster than the value of PUFAs, which results in the decrease in relative content of PUFAs and the ratio of PUFA/SFA [46]. Compared with the grass-fed beef, the beef fed on TMR possesses higher SFAs and n-6 PUFAs, Figure 6. The results of the Pearson correlation analysis for the fatty acids, SDLs and DEGs in yak's LD. Red and blue circles show the significant positive and negative correlation, respectively. When the color is darker the correlation is higher. * shows p < 0.05; ** shows p < 0.01.

Discussion
The deposition velocity of SFAs and MUFAs is faster than the value of PUFAs, which results in the decrease in relative content of PUFAs and the ratio of PUFA/SFA [46]. Compared with the grass-fed beef, the beef fed on TMR possesses higher SFAs and n-6 PUFAs, and lower n-3 PUFA [47]. The content of fatty acids in yak's LD under different feeding systems showed the SF can also decrease the relative content of Σn-3 PUFAs, ΣPUFAs, and the ratio of ΣPUFAs/ΣSFAs, whereas it increases the relative content of ΣSFAs. The World Health Organization and Food and Agriculture Organization of United Nations made a suggestion that the intake of SFAs and trans fatty acids (TFAs) should be decreased, while the intake of n-3 PUFAs should be increased. The ratio of ΣPUFAs/ΣSFAs in meat should be more than 0.4, while the ratio of Σn-6/Σn-3 PUFAs should be less than 4 [48,49]. Three TFAs (t-C18:1, t11-C18:1 and t-C18:2n6) were simultaneously determined in LD of yaks under GF and SF, and the absolutely content of t-C18:1 and t11-C18:1 in LD of SF yaks were higher than these values in the GF group. The ratio of PUFAs/SFAs and n-6/n-3 PUFAs in GF yak's LD were more than 0.4 and less than 4, respectively; whereas these values in LD of SF yaks were less than 0.4 and more than 4, respectively. Therefore, the meat of GF yaks is healthier than the meat of SF yaks from the viewpoint of nutritional value. Fat metabolism in ruminants is closely related to rumen microorganisms and the feeding system directly affects the rumen microorganisms. The biohydrogenation by rumen microorganisms in ruminants is the important way of UFAs transferring into MUFAs or SFAs [50]. The differences of rumen microorganisms in yaks under different feeding systems can cause the differences of biohydrogenation for UFAs in the diet, which is one of major factors that led to the differences in the fatty acids composition in yak muscle.
The lipidomics showed the fatty acids in yak's LD mainly existed in the form of triglyceride and phosphatidylcholine, and the major fatty acids in SDLs were C16:0, C16:1, C18:0, C18:1 and C18:2. Moreover, the determination of the fatty acids revealed the major ones in yak's LD under both SF and GF were C16:0, C16:1, C18:0, C18:1 and C18:2, and the absolute contents of C16:0, C18:0, C18:1 and C18:2 in GF and SF yak's LD were different. So the differences of fatty acids in SF and GF yak's LD were mainly due to C16:0, C18:0, C18:1 and C18:2. The DEGs in GY and SY yak's LD were mainly involved in the regulation of the fatty acids metabolism including fatty acids biosynthesis (bom00061), fat digestion and absorption (bom04975), amino acid metabolism including glycine, serine and threonine metabolism, protein digestion and absorption (bom04974) and carbohydrate metabolism including glycolysis/gluconeogenesis (bom00010), and carbohydrate digestion and absorption (bom04973). The fatty acids, amino acids and carbohydrate could transform to each other by the tricarboxylic acid (TCA) cycle [51,52]. Therefore, it can be inferred that the effect of the feeding system on the fatty acids in yak muscle was realized by the reciprocal transformation of three major macronutrients. The carbohydrate and amino acid intake by yaks from the grass or forage can transfer into short chain fatty acids, then the short chain fatty acids were further transformed to long chain SFAs. Moreover, the long chain SFAs can translate into MUFAs. In a word, the GF and SF yaks can intake different carbohydrates, amino acids, and fatty acids from their diet, which lead to the differences in fatty acids in yak muscle in the larger sense.
The regulation of fatty acid transport, fatty acid catabolism and fatty acid anabolism in bovines are driven by these genes in the PPAR signaling pathway [53]. The transcriptome analysis showed the differences of fatty acids in GS and SF yaks were mainly realized by the regulation of the PPAR signaling pathway too, and these downstream genes (FABP, FASN, ME1, SCD, ACBP, LPL, CPT1, ACSLs, ACAA1 and PLINs, which are also involved in the regulation of fatty acid metabolism in beef cattle [54,55]) in the PPAR signaling pathway were essential in the regulation of fatty acids metabolism in yak muscle. On the other hand, the SDLs were mainly enriched in phospholipid metabolism like glycerophospholipid metabolism (bom00564), glycerolipid metabolism (bom00564), triglyceride metabolism like fat digestion and absorption (bom04975), and UFAs metabolism like linoleic acid metabolism (bom00591), alpha-linolenic acid metabolism (bom00592), arachidonic acid metabolism (bom00590). Therefore, the fatty acids in yak muscle under different feeding systems were directly affected by the metabolism of UFAs, phospholipids and triglycerides regulated by the related genes in the PPAR signaling pathway.
As can be seen from the angle of molecular biology, the content of fatty acids in yak muscle is largely determined by some crucial enzymes in fatty acids metabolism. Identification of the genes encoding these enzymes can help understand the genetic variation underlying the content of fatty acids in yak meat. The PLIN1 gene negatively regulates the lipolysis in dairy cows [56]. The expression of the PLIN1 gene in yak muscle was positively correlated with the absolute content of ∑SFAs, ∑MUFA, Σn-6PUFAs, ΣUFAs and the ratio of Σn-6/Σn-3 PUFAs, so the PLIN1 gene in yak muscle can positively regulate all types of fatty acids content except Σn-3PUFAs, which is due to preventing fatty acid degradation. Biosynthesis of SFAs can be realized by the de novo synthesis of fatty acids and elongation of fatty acids, and biosynthesis of MUFAs can be realized by elongation of fatty acids and dehydrogenation of SFAs. The expressions of FASN, FABP3, SLC16A13 and DBI genes in yak muscle were positively correlated with the absolute content of ∑SFAs and ∑MUFAs, respectively. Moreover, the expression of the SCD gene in yak muscle was positively correlated with the absolute content of ∑MUFAs. The expression of the FASN gene was significantly associated with the concentrations of C14:0, C16:0 and C18:1 in cattle muscle [57]; the FABP3 plays an important role in the oxidation, esterification and metabolism of SFAs and MUFAs [58]; the DBI gene regulates the synthesis and degradation of ketone bodies, which is a key step in glucose and amino acid translating to SFAs [59]; the SCD gene is responsible for the dehydrogenation in the biosynthesis of MUFAs especially c-C18:1 in Holstein cattle [60]; the SLC16A13 gene positively regulates the long chain fatty acids transport which is necessary for the biosynthesis of long chain SFAs and MUFAs [61]. The absolute content of c-C18:1, ∑SFA and ∑MUFAs in SF yak's muscle were higher than values in GF yak's muscle. Therefore, FASN, FABP3, SLC16A13 and DBI genes in yak's muscle positively regulate the biosynthesis of SFAs and MUFAs by de novo synthesis of fatty acids, elongation of fatty acids, whereas the SCD gene in yak muscle was the only controlling gene in the process of SFAs translating to MUFAs. Some complicated PUFAs cannot directly be synthesized in liver and adipose tissue of livestock and must be transferred from the precursor compounds PUFAs from grass and feed [62]. The n-6 PUFAs like arachidonic acid are derived from linoleic acid by the action of desaturases and elongases, and n-3 PUFAs like DHA and EPA are derived from linolenic acid. The FADS6 gene acts in the dehydrogenation of linoleic acid and linolenic acid, and plays a crucial role in PUFA biosynthesis [63,64]. The concentrations of C18:3n3, c-C20:5n3 and Σn-3PUFAs in SF yak muscle were lower than the values in GF yak muscle. Meanwhile, the expression of the FADS6 gene in yak muscle was positively correlated with the absolute content of Σn-3PUFAs, and negatively correlated with the value of Σn-6/Σn-3 PUFAs. Therefore, the FADS6 gene is a crucial candidate gene for regulating the Σn-3 PUFAs content and the ratio of Σn-6/Σn-3 PUFAs in yak muscle.

Conclusions
The characteristics and regulation of fatty acids in yak muscle under GF and SF were investigated in this study. The absolute content of ΣSFAs, ΣUFAs, ΣPUFAs, ΣMUFAs and Σn-6 MUFAs in SF yak's meat were higher than these values in GF yak's meat, whereas the relative content of Σn-3PUFAs, ΣPUFAs, ΣMUFA and ΣUFAs in SF yak's meat were lower. The ratio of Σn-6/Σn-3 PUFAs, ΣPUFA/ΣSFA in GF yak's meat was less than 4, and more than 0.4, respectively, thus the composition of fatty acids in GF yak's meat is healthier for consumers from the point of nutritional value. The metabolism of fatty acids, amino acids, carbohydrates and phospholipids in GF and SG yak's meat were different, and the differences of fatty acids in yak meat were mainly regulated by the FASN, FABP3, PLIN1, SLC16A13, FASD6, DBI and SCD genes in the PPAR signaling pathway.

Supplementary Materials:
The following supporting information can be downloaded at https:// www.mdpi.com/article/10.3390/foods11172582/s1. Table S1: The information of primers for qPCR determination. Table S2: The information of total 75 significant differences lipids (SDLs) in yaks longissimus dorsi (LD) muscle between graze feeding (GF) and stall feeding (SF) regimes. Table S3: The information of 17 KEGG pathways enriched by SDLs. Table S4: The information of differentially expressed genes (DEGs) in yak's LD between GF and SF. Table S5: The information of KEGG pathway enriched by DEGs.