Muscle transcriptome analysis reveal candidate genes and pathways related to fat and lipid metabolism in Yunling cattle

Abstract Yunling cattle (YL) is a recently developed beef breed harboring a quarter of Yunnan ancestral cattle genome, spanning over past 30 years. Compared with Diqing cattle (DQ), a Yunnan native cattle breed, YL presents various advantages, including rapid growth and exquisite meat quality. However, the molecular mechanisms underlying these phenotypic differences are not clearly understood. To further identify the candidate genes responsible for the quality of the meat in the muscle, longissimus dorsi (LD) muscle was used for RNA-Seq analysis. A total of 508 differentially expressed genes (DEGs) were identified in YL (adjusted p-value <0.01 and log2FoldChange >1), of which 243 were up-regulated and 265 were down-regulated. Functional association analysis showed that the identified DEGs mainly enriched the lipid and fat metabolism pathways. Moreover, it was also observed that several fat-related genes were differentially expressed in both cattle breeds, including three up-regulated genes (MOGAT1, ACSM3, PLPP2) and two down-regulated genes (ADIG, GPAT3). In addition, alternative splice analysis was also performed revealing an important 9–11 exon skipping variation of GPAM gene (crucial for beef marbling) in YL, which is three times higher than that in DQ, suggesting that this variation might have played the central role in the ‘snow beef’ effect in YL. We believe that our results will help in understanding the mechanism of muscle development and promote the further breeding programs in YL cattle.


Introduction
Yunling cattle (YL) is a newly developed beef breed harboring well-balanced body stature and stable genetic performance, established by the Yunnan Academy of Grassland and Animal Science in China for more than 30 years. 1 The breeding process of YL includes 1/ 2 Brahman cattle, 1/4 Murray Gray and 1/4 Yunnan native cattle, harboring three-way crossbreeding and inter se breeding. 2 It is characterized by its enhanced tolerance against warm-wet climate, rapid growth, and improved fat deposition. In addition, the 'snow beef' trait related to the marbling of the meat in YL is popular with people as a high-end consumer product. YL is reportedly better than Wenshan cattle (Yunnan native cattle) in various traits, such as body weight and body size index. In addition to the phenotype differences, the transcriptomic analysis also revealed differences in two important signaling pathways related to meat quality. 3 Previously, the studies about YL have mainly focused on its ancestral 1 and maternal lineage, 4 however some studies have also explained DNA methylation of the meat quality related genes 5 and growth-related pituitary differentially expressed genes. 6 Nevertheless, there is still a lack of transcriptomic data on the regulatory mechanism associated with the longissimus dorsi (LD) muscle in YL, which affects the further selection and breeding process in YL.
Beef is an important dietary source of protein and fat in human diet. 7,8 With the improvement of consumption levels, high-quality meat (including flavor, tenderness and juiciness) is also sought after by consumers. 9,10 Previously, various genes related to meat quality have been identified, such as, FABP4 and ADIPOQ have been found to be associated with fat deposition in Shandong black cattle, 11 SIRT1 gene related intramuscular fat content in Qinchuan cattle. 12 The content of intermuscular fat (IMF) directly affects the quality of meat; therefore, an appropriate amount of IMF is a crucial indicator in beef breeding. 13 IMF deposition or marbling in beef meat, especially in YL, is one of the most important traits, affecting the color and palatability attributes of beef. 14 The intense and well-marbled beef is also named 'snow beef' for its tenderness and juicy beef. Recently, a transcriptomic study in YL revealed a lower IMF and monounsaturated fatty acids, while higher proportions of saturated fatty acids, polyunsaturated fatty acids, and short-chain fatty acids were observed in its LD muscle. 15 However, there is limited to insufficient research on muscle fat deposition in YL, which needs to be addressed.
RNA-Seq was used for the transcriptomic study, in order to explain the differentially expressed genes (DEGs) and differentially alternative splicing (DAS) affecting the quality of meat in the LD muscles of YL. The subsequent analysis divulged some candidate genes related to IMF deposition in YL. Our research provides a powerful transcriptomic reference for the further genetic improvement of fat deposition in YL.

Ethics statement
All animal based experimental protocols were approved by the biological studies animal care and use committee at institute of animal science, Northwest A&F University, Shaanxi, China (Permit number: NWAFAC1019).

Animals and sample collection
In the current experiment, 6 steers (3 Yunling cattle and 3 Diqing cattle, Table S1) from same feeding conditions for the whole experimental period (from birth to 30-Month-old) till slaughter were used in the Xiaoshao Farm of Yunnan Academy of Grassland and Animal Science. The complete sample set was made sure of any direct or collateral family relationship for at least three consecutive generations. After slaughter, the specimens were collected from the longissimus dorsi (LD) muscle between the 12th and 13th ribs (left half carcass), immediately stored in liquid nitrogen, and transferred to À80 C before RNA extraction.

RNA isolation and illumina sequencing
The total RNA of each muscle sample was extracted using Trizol reagent (Sangon, Shanghai, China), according to the manufacturer's instructions. The RNA integrity was verified by measuring the Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA). Following Illumina RNA-Seq library preparation, paired end sequencing (150 bp) of three muscles from each group was performed on Illumina HiSeq-1000 Platform.

Analysis of sequence data
After sequencing, raw data was cleaned using Trimmomatic software 16 (LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36). The clean data was aligned to the latest Bos taurus reference genome assembly (ARS-UCD1.2) using HISAT2 (version 2.1.0). 17 To obtain the Fragments Per Kilobase of transcript per Million mapped reads (FPKM), StringTie software 18,19 was used for the quantification of transcripts and genes. A python script (prepDE.py) was used to convert StringTie 18,19 result into DESeq2 20 analyzable data. Finally, the DESeq2 20 package in R was used for data normalization and analysis of DEGs. The adjusted p-value was used to correct the p-value with the false discovery rate (FDR). The adjusted p-value <0.01 and absolute value of log2FoldChange >1 were used as the threshold to statistically define the significant DEGs. We used DESeq2 20 to perform principal component analysis (PCA) and visualize the gene expression profiles of each individual.

Enrichment analysis
The DEGs and DAS gene set enrichment analysis was carried out by gene ontology (GO) analysis, and Kyoto Encyclopedia of Genes and Genomes (KEGG) using the online KOBAS annotation tool (v.3.0; Beijing, China; http://kobas.cbi.pku.edu.cn). The significant pathways and genes were selected stringently with an adjusted probability (Fisher's exact test, p < 0.05).

Analysis of DAS genes
The rMATS, 21 a computational tool to detect the differential alternative splicing (DAS) events from RNAseq data, was employed to detect the DAS events between YL and DQ. The DAS events are classified as alternative 5 0 splice site (A5SS), exon skipping (ES), mutually exclusive exons (MXE), alternative 3 0 splice site (A3SS), and retained intron (RI). Significant DAS genes were screened based on the threshold of jDWj>0.1 and false discovery rate (FDR) 0.05. 22

Overview of RNA sequencing data
In total, $431.3 million clean reads were generated from Illumina HiSeq-1000 Platform. The Q30 value of an average 150 bp clean read ranged around 94% whereas, the GC content of each sample lied between 49 and 51%, indicating the reliability and high quality of the sequenced data ( Table 1). The clean reads were mapped to the latest Bos taurus reference assembly ARS-UCD1.2 using the HISAT2 software. 17 Approximately 98% (94 $ 99%) of each individual was mapped to the reference genome while 70% (60 $ 80%) of reads were uniquely mapped to the genic regions of the reference genome (Table 1).

Differentially expressed gene analysis
The box plot distribution of Fragments Per Kilobase of exon model per Million (FPKM) intuitively compared the overall gene expression level differences between different samples, depicting similar expression patterns of YL and DQ (Fig. 1A). Principal component analysis showed that there was a significant difference between both groups, suggesting that the transcriptomic genes in the LD muscle tissue had different expression pattern in YL and DQ (Fig. 1B). A total of 508 candidate genes were differentially expressed in comparison to YL and DQ in the LD muscle (Fig. 1C). Among these DEGs, 243 were up-regulated, while 265 were downregulated genes in YL (adjusted p-value <0.01 and log2FoldChange >1) (Fig. 1D, E).
To further investigate the biological function of these DEGs, GO and KEGG analysis (v.3.0, Beijing, China; http://kobas.cbi.pku.edu.cn) were conducted separately for both the up and down-regulated genes (Fig. 2). The enrichment analysis represented 20 significantly enriched KEGG pathways and GO terms, each. The enrichment was filtered out with corrected p-values <0.05 ( Fig. 2; Table S2, 3). Some of the significant KEGG pathways were associated with the lipid metabolism (TGF-beta signaling pathway, corrected p-value ¼ 0.001; glycerolipid metabolism, corrected p-value ¼ 0.031) and metabolic pathways (corrected p-value ¼ 0.003). Moreover, three significant GO terms were found to be associated with fat metabolism (fatty acid biosynthetic process, corrected p-value ¼ 0.002; positive regulation of fat cell differentiation, corrected p-value ¼ 0.008; cellular lipid metabolic process, corrected p-value ¼ 0.009).

Identification and characterization of DAS events
Alternate splicing (AS) events produces protein subtypes with different functions, increasing the diversity in the gene expression regulation. rMATS software 21 was used to understand the alternative splicing functions of the muscle tissue in YL and DQ from the transcriptomic data. A total of 42,831 AS events were curated from 15,585 genes. Exon skipping was the most frequent AS event while retained introns were the rarest AS events. Amongst all AS events, 1,431 AS events and 1,148 genes showed significant differences between YL and DQ ( Table 2).
To further define the DAS genes biological functions, KEGG and GO functional enrichment analysis (v.3.0, Beijing, China; http://kobas.cbi.pku.edu.cn) was performed. The results of KEGG pathways analysis showed significant enrichment of glycerophospholipid metabolic pathways in DAS genes (Table S4). Interestingly, we also found muscle-related GO biological processes including 'microtubule cytoskeleton organization', 'calmodulin binding' and 'actin filament' (Table S5).
Moreover, an ES event was identified to the GPAM gene, on chromosome 26. Differential splicing skipped exon 10, resulting in a direct connection between exon 9 and exon 11 (Fig. 4A). Although, the splicing mutations are 3Â in YL, however, there was no significant differential expression of the GPAM gene between YL and DQ (Fig. 4B). The existence of the ES event was confirmed through the Ensemble database (Fig. 4C). Also, GPAM gene is related to fat metabolism, suggesting that this gene may be potentially related to the meat quality in YL.

Discussion
Balanced intermuscular fat (IMF) is crucial to increase the marbling beef production for beef breeders. As a commercial beef cattle, the excellent meat quality of Yunling cattle (YL) is praised by peoples, but there are only a few reports on the potential molecular mechanism of the meat quality of YL. In this regard, the longissimus dorsi (LD) muscles of YL and DQ were processed for transcriptomic data generation. The RNA-Seq analysis identified candidate differentially expressed genes (DEGs) and alternate splicing (AS) events, which provide references for meat-quality trait improvement in YL.
To further understand the mechanisms underlying regulation of 508 DEGs, we performed GO and KEGG pathway analyses. As expected, several wellknown GO terms and signaling pathways relating to lipid and fat metabolism were identified, including fatty acid biosynthetic process (6 DEGs), positive regulation of fat cell differentiation (6 DEGs), cellular lipid metabolic process (4 DEGs), TGF-beta signaling pathway (10 DEGs) and glycerolipid metabolism pathway (6 DEGs). 15,[23][24][25][26] Large numbers of DEGs involved in those signaling pathways have been proven to be functional in lipid metabolism, such as BMP5, THBS1, MOGAT1, ACSM3, PLPP2, ADIG and GPAT3. [27][28][29][30][31][32] In this study, multiple GO terms and signaling pathways related to lipid and fat metabolism were screened, suggesting that the longissimus dorsi muscle of YL and DQ has significant difference in fat deposition. Furthermore, DAS genes are significantly   enriched in glycerophospholipid metabolic pathways (14 DASs), suggesting that the differential alternative splicing have important roles in fat biosynthesis and metabolism. 33,34 In the present study, we identified three significantly up-regulated genes (MOGAT1, ACSM3, PLPP2), all of which were involved in adipogenesis. Among them, the MOGAT1 (Monoacylglycerol O-Acyltransferase 1, MGAT1) gene plays an important role in the process of triglyceride synthesis. MOGAT1 À/À mice lead to increased lipogenesis suggesting that this gene plays an important role in fat decomposition and absorption. 30 A recent study pointed out that the polymorphism of this gene is related to the growth and development of beef cattle, and the homozygous mutation of g.25940 T > C leads to a significant increase in chest depth in YL. 31 Furthermore, ACSM3 (Acyl-CoA Synthetase Medium Chain Family Member 3) is an acyl-CoA synthetase which takes part in the first step of fatty acid metabolism. 32 In addition, PLPP2 is a protein-coding gene and belongs to a phospholipid phosphatase family. In our results, these three genes were significantly upregulated in the muscle of YL, suggesting that these genes are one of the potential reasons for the higher fat content in YL than that in DQ. Moreover, two significantly down-regulated genes (ADIG, GPAT3) were also identified. ADIG (Small Adipocyte Factor 1, SMAF1) gene specifically expressed in the plasma membrane of adipocyte, 35 which is a protein coding gene that regulate the differentiation and proliferation of adipose tissues. 36 Overexpression of ADIG in vitro, promotes fat accumulation. 37 In contrast, the absence of ADIG impairs adipogenesis. 38 In addition, GPAT3 is a glycerol-3-phosphate acyltransferases gene in the white adipose tissue and plays an important role in regulating lipid homeostasis. 39 These results suggest that ADIG and GPAT3 may be associated with a high quality 'snow beef' of YL.
Alternative splicing is a crucial stage of RNA processing and increases protein diversity through complex splicing mechanisms. 40,41 In cattle, at least one-third of the genes are alternatively spliced. 34 In a study conducted by Kim et al., 42 four miRNA of the GPAM gene are significantly related to the quality of Korean beef. GPAM gene as a candidate gene related to the marbling score of cattle. 43,44 In this study, we report a vital AS event of GPAM gene as a potential gene related to the marbling score of YL 'snow beef'. Our research adds the way that GPAM increases the intramuscular fat for cattle.
YL is a developed beef breed with rapid growth, excellent meat quality, and enhanced resistance to humid and hot weather. The current study revealed the differential genes and alternative splicing events between YL and DQ muscle tissue at the transcriptome level. We have discovered important genes related to intramuscular adipogenesis and ES event related to YL 'Snow Beef'. Our research provides some important information for the further genetic breeding in YL.

Conclusion
In summary, we compared DEGs and DAS genes in the muscle transcriptome of YL and DQ. We identified 5 DEGs (MOGAT1, ACSM3, PLPP2, ADIG, GPAT3) associated with intramuscular fat of YL. In addition, an exon skipping of GPAM is significantly related to the meat quality in YL. The results will provide significant assistance in the future investigation of the molecular mechanisms of YL meat quality traits.

Acknowledgement
Finally, we thank High-Performance Computing (HPC) of Northwest A&F University (NWAFU) for providing computing resources.

Disclosure statement
No potential conflict of interest was reported by the author(s). The authors declare that they do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.