Variations in the diacylglycerol acyltransferase-1 (DGAT1) and its association with meat tenderness in Gannan yaks (Bos grunniens)

Abstract Diacylglycerol acyltransferase-1 (DGAT1) has been recognised as one of the functional genes during fat deposition and closely related to meat quality. It is well known that yak meat has low intramuscular fat (IMF) and poor tenderness, influencing the purchase decision for consumers. In this study, the variations of DGAT1 gene were detected using polymerase chain reaction single-stranded conformational polymorphism (PCR-SSCP) analysis so as to ascertain the association of DGAT1 variation with yak carcase and meat quality traits. Two variants (A1 and B1 ) in intron 1 – exon 2 and another three variants (A2 to C2 ) in intron 15 – exon 17 were identified in the yaks investigated. Three sequence variations (c.192-117 G > A, c.1252-23 C > T and c.1339 C > T) were detected among these variants, with c.192-117 G > A located in intron 1, c.1252-23 C > T located in intron 15 and c.1339 C > T located in exon 17. The mutation c.1339 C > T resulted in an amino acid change (p. Arg 447 Cys). Variants A1 and B2 were associated with a decrease and increase in Warner–Bratzler shear force (WBSF) (p = 0.027 and p < 0.001), respectively. Variant C2 was associated with a decrease in WBSF and hot carcase weights (HCW) (p = 0.001 and p = 0.037, respectively). The presence of haplotype H3 and absence of H5 were highly associated with a decrease in WBSF (p < 0.05). These finding suggested that the variations in DGAT1 could be potential targets for gene-assisted selection to improve meat tenderness in Gannan yak. Highlights Three SNPs were identified in the DGAT1 gene in Gannan yak. Variants (A1 , B2 and C2 ) and haplotypes (H3 and H5) have a significant influence on Warner–Bratzler shear force (WBSF) in Gannan yak. DGAT1 variations may have potential in selection for improving meat tenderness in Gannan yaks.


Introduction
Carcase and meat quality traits play an economically important role in the value assessment of the beef cattle industry. Among these traits considered, tenderness is a key factor in determining meat quality in the beef industry (Fiems et al. 2000). Gannan yak is a unique indigenous Bos grunniens breed well adapted to high altitudes and hypoxic conditions, which is located in the Southwest of the Tibetan Plateau in China. Yak meat has a high protein and minerals content that is beneficial for human health and would be a new resource for meat industry (Tian et al. 2013). However, yak breeds are known to produce meat with low intramuscular fat (IMF) and poor tenderness.
Meat tenderness is positively correlated with fat content and composition (Hocquette et al. 2010;Nishimura 2015). Fat deposition is a complex biological process that is regulated by controlling the activity of enzymes in lipid synthesis and metabolism. Triglycerides (triacylglycerols, TGs) are the major storage form of lipid droplets in adipocytes (Chitraju et al. 2019). Diacylglycerol acyltransferases (DGATs), including DGAT1 and DGAT2, are key enzymes for catalysing the terminal step in the formation of TGs through the acylation of diacylglycerol (DAG) with a fatty acyl-CoA and thus regulating lipid digestion, absorption, and glycerol lipid metabolism pathways (Bhatt-Wessel et al. 2018). DGAT1 is nearly ubiquitous in eukaryotes and widely expressed in tissues (Turchetto-Zolet et al. 2011). Studies have shown that DGAT1 is highly expressed in tissues with TGs synthesis and metabolism, such as subcutaneous fat, longissimus dorsi muscle, breast and small intestine (Cases et al. 1998;Chen 2006). It is well documented that mice lacking the DGAT1 enzyme may have a 50% reduction in TGs storage in adipose tissue (Smith et al. 2000) and a slower rate of fat absorption in the intestines (Hung and Buhman 2019). Thus, DGAT1 plays a pivotal role in fat deposition.
DGAT1 is located within a quantitative trait locus (QTL) for fat deposition mapped to bovine chromosome 14 (Casas et al. 2003). In previous studies, the expression of DGAT1 was positively correlated to IMF content in Angus cattle and pigs (Anton et al. 2011;Cui et al. 2019). Variation in DGAT1 has been reported to closely relate to mutton tenderness (Xu et al. 2009) and level of beef marbling , as well as carcase weight in sheep (Armstrong et al. 2018).
Currently, relatively little is known about the effect of DGAT1 on fat-related traits of muscles in Gannan yak. Therefore, the objective of this study was to identify genetic variations of the DGAT1 gene and to evaluate any associations with carcase and meat quality traits to assist with breeding and selection to improve meat quality in Gannan yak.

Experimental animals and genome DNA extraction
A total of 616 Gannan yaks, obtained from twelve populations that grazed on the plateau pasture under the same feeding and management conditions all year round in Gannan Tibetan Autonomous Prefecture (Gansu, China), were investigated in this study. All yaks ranged in age from three to seven years. When yaks investigated have been slaughtered in slaughterhouse in October, gender and age were recorded, and blood samples were collected from each yak onto Flinders Technology Associates ( (Zhou et al. 2006) from FTA cards.
All yaks were humanely sacrificed to alleviate suffering and handled in accordance with an animal use protocol approved by the animal use committee of the Chinese Ministry of Agriculture (Approval number 2006-398).

Carcase and meat quality traits measurements
Hot carcase weights (HCW; kg) from 301 yaks were measured immediately after slaughter. A cross-section of the longissimus dorsi muscle between the twelfth and thirteenth rib of the right carcase side was used to measure the ribeye area (REA; cm 2 ), and meat samples were taken from longissimus dorsi muscle to assess the meat quality, at 48 h post-mortem. Warner-Bratzler shear force (WBSF) was determined according to Fortes et al. (2009) as follows. Meat samples were cooked until reaching an internal temperature of 71 C, and then removed from the water bath and cooled for 24 h to 5 C-6 C. A total of six cores (1.27 cm diameter) were obtained from each sample parallel to the muscle fibre orientation using a sampler. The analysis using the Digital Muscle Tenderometer (Model C-LM3, Northeast Agricultural University, Harbin, China) that had crosshead speed force set at 200 mm/min to measure the WBSF (kg). The 6 core values were used to calculate mean WBSF value for each sample. Meat samples used for drip loss rate (DLR; %) measurement were cut from the longissimus dorsi muscle and immediately weighed (initial weight). Then, samples were suspended in an inflated bag and were again weighted after a storage period at chill temperature. At the time of cooked meat rate (CMR; %) measurement, the initial weight of samples were firstly weighted. Then, samples were heated in a waterbath. When the endpoint temperature has been attained, the samples were removed and cooled in an ice slurry and used to weight ( Honikel.1998).

Polymerase chain reaction (PCR) amplification and Single-Stranded conformational polymorphism (SSCP) analysis
To rapidly detecting sequence variation in DGAT1 gene, genomic DNA from 50 yak individuals were randomly selected and used as a template for polymerase chain reaction (PCR) amplification. Eight sets of primers (Table 1) for PCR were designed based on the bovine sequence (GenBank no. AY065621.1) to amplify the entire CDS region and some introns of DGAT1. All primers were designed by Primer-BLAST (NCBI https:// www.ncbi.nlm.nih.gccov) and synthesised by Huada Biological Engineering Co., Ltd. PCR amplification was performed in a 20 lL reaction mixture, containing 100 ng of genomic DNA (one 1.2-mm punch of blood spot on FTA card), 0.25 lM of each primer, 150 lM of each dNTP (Takara, Japan), 2.5 mM Mg 2þ , 0.5 U of Taq DNA polymerase (Takara, Japan), 2.0 lL of 10 Â PCR buffer and ddH 2 O to make up the volume. The cycling protocol was initial denaturation at 94 C for 5 min, followed by 35 cycles of denaturation at 94 C for 30 s, annealing at 57-64 C (Table 1) for 30 s, and extension at 72 C for 30 s and a final extension at 72 C for 5 min. PCR products were checked using 1% agarose gel electrophoresis and then were sequenced in both directions at Sangon Biotech (Shanghai, China).

Variant sequencing and analysis
PCR amplicons that were confirmed to be homozygous by SSCP were directly sequenced in both directions at Sangon Biotech (Shanghai, China). Variants that were only found in heterozygous yaks were sequenced using an approach described by Hu et al. (2010). Briefly, a band corresponding to the variant was cut as a gel slice from the polyacrylamide gel, macerated, and then used as a template for re-amplification with the original primers. This second amplicon was then sequenced directly as described above for the homozygous patterns. Nucleotide sequence alignments were analysed by MEGA 5.0 software (CEMI, Tempe, AZ, USA) to determine the single-nucleotide polymorphisms (SNPs) location.

Haplotype determination
Haplotypes were constructed for intron 1exon 2 and intron 15exon 17 of DGAT1. Progeny that typed as homozygous in either of the regions could directly infer their haplotypes. If progeny typed as heterozygous in both regions, the haplotypes could not be analysed. This reduced the number of yaks studied in the haplotype association analyses.

Statistical analysis
Association analyses of DGAT1 variants and haplotypes with carcase and meat quality traits were performed using the General Linear Mixed Models (GLMMs) of IBM SPSS 22.0 software (IBM Corp., Armonk, NY, USA). The model was calculated as follows: where Y ijk is the trait of interest, l is the overall mean, G i (H i ) is the fixed effect of the ith variant (i ¼ 0 or 1) or the fixed effect of the ith haplotype (i ¼ 0 or 1) or the fixed effect of the ith number of copies of the haplotype ((i ¼ 0, 1, 2), P j is the fixed effect of the jth population (j ¼ 1, … ,12), S k is the fixed effect of the sex (k ¼ 1(male), 2 (female)), A l is the covariate effect of age (l ¼ 3, … ,7), and e ijkl is the random error. Each variant/haplotype was coded as either presence (1) or absence (0) for each yak investigated. Any variant and haplotype that had associations in the single-variant/haplotype models with a p-value of less than 0.20, and thus could potentially impact on the carcase and meat quality traits, were factored into the subsequent multi-variant/haplotype models. For these haplotypes with sufficiently common homozygous forms (>1% of all genotypes), a second set of analyses was performed with the number of haplotype copies present included (in place of presence or absence), followed by planned orthogonal contrasts to ascertain whether additive, dominant or recessive effects were present. These models were conducted in an identical manner to the GLMMs used for testing the presence/absence of each haplotype. Unless otherwise indicated, all p-value less than 0.05 were considered to be significantly different.

Identification of variation in yak DGAT1
In 616 Gannan yaks investigated, three distinct genotypes A 1 A 1 , A 1 B 1 and B 1 B 1 were detected in intron 1exon 2 and a further five distinct genotypes A 2 A 2 , A 2 B 2 , A 2 C 2 , B 2 B 2 and C 2 C 2 were detected in intron 15exon 17 by PCR-SSCP (Figure 1(A)). Three SNPs were detected among those sequences, with c.192-117 G > A in intron 1, c.1252-23 C > T in intron 15 and c.1339 C > T in exon 17 (Figure 1(B,C)). Moreover, the c.1339 C > T was a non-synony mous SNP and would result in an arginine to cysteine substitution (p. Arg 447 Cys). All sequences were submitted to GenBank (accession number MN012932 to MN012936).

Polymorphisms in yak DGAT1
The genotype and variant frequencies of yak DGAT1 are shown in Table 2. In intron 1exon 2, A 1 A 1 and A 1 were the most common genotype and variant, with a frequency of 50.65% and 70.13%, respectively. In intron 15exon 17, A 2 A 2 and A 2 were the most common genotype and variant, with a frequency of 37.34% and 57.55%, respectively.

Construction and frequencies of haplotype in yak DGAT1 gene
Six haplotypes were constructed for intron 1exon 2 and intron 15exon 17 of DGAT1 in a subset of 519 individuals out of 616 Gannan yak (Table 3). Of these haplotypes, H1, H2, H3 and H4 were the most common, occurring at frequencies of 42.20%, 14.45%, 17.24% and 16.28%, respectively. H5 and H6 haplotypes were rare, with frequencies of 5.68% and 4.14%, respectively.

Association between DGAT1 variants and carcase and meat quality traits
In the single-variant (presence/absence) model (Table 4), the presence of variant A 1 (in intron 1) was found to be associated with decrease WBSF (p ¼ 0.027). The presence of variant B 2 (in exon 17) was associated with an increase in WBSF (p < 0.001), whereas the presence of C 2 (in intron 15) was related to decrease WBSF and HCW (p ¼ 0.001 and p ¼ 0.037, respectively). These associations remained significant when the other variants (p < 0.20) were factored into the models. No other significant associations were detected for any variant in both regions of yak DGAT1 (p > 0.05).

Association between haplotypes in DGAT1 and carcase and meat quality traits
Haplotype frequencies lower than 5% were excluded from association analyses. In the single-haplotype (presence/absence) model (Table 5), the presence of haplotype H3 and the absence of haplotype H5 were associated with a decrease in WBSF (p ¼ 0.002, p < 0.001).These associations persisted in the multihaplotype models. No associations were found between the other haplotypes and DLR, CMR, REA and HCW (p > 0.05). A second set of analyses was performed with the number of haplotype copies present (in place of presence/absence). For WBSF traits, the presence of two copies of H3 and the absence of H5 were associated with decrease WBSF in both singlehaplotype and multi-haplotype models (Table 6).

Discussion
Most of the meat quality traits, especially meat tenderness, have been presenting genetic differences among cattle populations (Xie et al. 2012;Boudon et al. 2020). Numer ous studies have been shown that calpain and calpastatin genotypes were associated with meat tenderness in different cattle populations (Casas et al. 2006;Curi et al. 2009;Allais et al. 2011;Tait et al. 2014) because of their key role in proteolysis. However, meat tenderness is under polygenic control. Identifying more gene markers related to meat quality traits would improve meat tenderness prediction. DGAT1 is an important functional enzyme involved in TGs biosynthesis and metabolism, and the gene encoding this protein has been considered to be a genetic factor influencing IMF deposition and meat tenderness in cattle (Thaller et al. 2003;Harris et al. 2011). This is the first study to report the associations between DGAT1 variants or haplotypes and carcase or meat quality traits in yaks and the findings suggest that these genetic variations may play an important role in yak meat tenderness trait.    Studies have shown that 17 and 5 SNPs in DGAT1 were detected in bovine and sheep, respectively (Scat a et al. 2009;Yuan et al. 2013). Also, Yuan et al. (2007) reported that 7 SNPs were observed in the buffalo DGAT1 gene, and the c. 1481 C > T (originally g.11,785 C > T) was located in exon 17, potentially resulting in an alanine to valine substitution (p. Ala 484 Val). In this study, a total of 3 SNPs c.192-117 G > A, c.1252-23 C > T and c.1339 C > T were detected in intron 1, intron 15 and exon 17 of DGAT1 in Gannan yak, respectively. Substitution c.1339 C > T would potentially result in the 447th amino acid in the mature peptide from Arg to Cys.
DGAT1 is the key enzyme to catalyse TGs synthesis and the DGAT1 variations have been reported to be associated with meat quality traits (Thaller et al. 2003;Pannier et al. 2010). The K232A substitution at exon 8 of the DGAT1 was associated with IMF content and beef marbling level . Yuan et al. (2013) investigated c.1416 T > G (p. Gly 468 Val) in exon 17 of bovine DGAT1 gene and reported that the sequence variation had an effect on fatness deposited, WBSF and marbling score. Xu et al. (2009) also reported that the c.1461 T > C in exon 17 of DGAT1 in sheep was positively correlated with WBSF, DLR, and negative with IMF and muscle marbling score.
Similarly, Mohammadi et al. (2013) found that the variations in exon 16 À 17 of DGAT1 were significantly associated with back fat thickness and other fatrelated traits in Iranian sheep. In this study, the presence of variant B 2 in exon 17 was associated with an increase in WBSF, and then had a negative effect on meat tenderness in Gannan yak. The significance of the substitution of Arg to Cys is unknown, but given that cysteine is essential for the formation of disulphide bonds (Wang et al. 2017), it is conceivable that this amino acid substitution may affect the DGAT1 protein folding. It is reported that the p. Arg 92 Cys substitution in colipase impaired its function and secretion by increasing protein misfolding (Xiao et al. 2013). Therefore, we speculate that the c.1339 C > T SNP may affect the structure and function of the DGAT1 protein, nevertheless further investigation is required.The association between intron variations of DGAT1 and the fat-related traits has been reported in numerous studies (Angiolillo et al. 2007;Armstrong et al. 2018). For instance, the SNP in intron 16 of DGAT1 may influence milk fat content and other economical traits in goat (Angiolillo et al. 2007). Armstrong et al. (2018) showed that c. 512 þ 411 C > T (originally rs411875883) in intron 6 of DGAT1 in sheep was highly associated with REA, HCW, and fat thickness. Consistent with these previous studies, the presence of variant A 1 (in intron 1) and C 2 (in intron 15) were associated with a decrease in WBSF and thus positively correlated with yak meat tenderness. In addition, variant C 2 was also associated with HCW in Gannan yak. Although these variations were not located in the coding regions of DGAT1, they may impact gene function in other ways. There is some evidence that variations in introns may not directly affect the structure of gene, but could influence the transcription efficiency by affecting regulatory elements, such as enhancers, silencers, or other DNA structures (Chorev and Carmel 2012). Moreover, variation in introns may be linked to other variation in coding regions of the gene or regions that may affect the processing of the primary transcript or its stability . Scheike et al. (2010) reported that molecular markers could be accurately identified when integrating the SNPs and haplotypes into association analysis. Remarkably, haplotypes are particular combinations of variant sequences on a single chromosome (Hu et al. 2019) and may be considerably more powerful than single SNP analysis in association studies (Hagenblad et al. 2004). Aslan et al. (2010) reported that the HAP2 was association with increased IMF and juiciness and HAP4 was correlated with tenderness in bovine, and these two haplotype in Ankyrin 1 gene could be used to select individuals with improved IMF, juiciness or tenderness. In current investigation, haplotype were constructed with 3 SNPs, and the presence of haplotype H3 and absence of H5 were highly related to a decrease in WBSF, suggesting that H3 and H5 may have potential in selection for yak meat tenderness. Meat tenderness is a complex quantitative trait regulated by multiple genes. Therefore, it is necessary to increase the sample size or study the effects of other related genes on tenderness to determine the common effect of DGAT1 gene and other related genes on meat tenderness.

Conclusions
In summary, three SNPs were identified in the DGAT1 gene in Gannan yak. Variants and haplotypes of intron 1, intron 15 and exon 17 in DGAT1 were significantly associated with WBSF of the longissimus dorsi muscle and had an important effect on yak meat tenderness. These results suggest that the DGAT1 variations may have potential in selection for improving meat tenderness in Gannan yaks.