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 (A₁ and B₁) in intron 1 – exon 2 and another three variants (A₂ to C₂) 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 A₁ and B₂ were associated with a decrease and increase in Warner–Bratzler shear force (WBSF) (p = 0.027 and p < 0.001), respectively. Variant C₂ 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 (A₁, B₂ and C₂) 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.

Keywords:

• DGAT1
• variations
• meat tenderness
• Gannan yak

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 (Li et al. 2013), 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.

Materials and methods

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 (FTA) cards (Whatman BioScience, Middlesex, UK). Genomic DNA was .............................................................................................. .......................................................................................................................... ............................................................................................... obtained using a two-step procedure (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²), 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 water-bath. 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 μL reaction mixture, containing 100 ng of genomic DNA (one 1.2-mm punch of blood spot on FTA card), 0.25 μΜ of each primer, 150 μΜ of each dNTP (Takara, Japan), 2.5 mM Mg²⁺, 0.5 U of Taq DNA polymerase (Takara, Japan), 2.0 μL of 10 × PCR buffer and ddH₂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).

Table 1. PCR primers for sequencing in yak DGAT1.

Primer	Amplified region	Sequence of primer (5'–3')	Length of amplicon (bp)	PCR annealing temperature (°C)
P1	Exon 1	ATGGGCGACCGCGGCG	303	64
		CCAGACTCTCGGACTC
P2	Intron 1 – exon 2	GTCTTGCATCACCAGCTCCT	248	62
		GTACCAGCATCACCACACAC
P3	Exon 3 – exon 5	GATCTTAAGCAACGCACGGT	444	60
		TGAAGCACAGTGGGACCAG
P4	Intron 5 – exon 7	CTGCACCTCAGGCCTGGT	446	63
		GAGGCAAGTGGCAGCCC
P5	Exon 8 – intron 8	GTAGCTTTGGCAGGTAAG	201	57
		GGGAAGTTGAGCTCGTAG
P6	Exon 9 – intron 12	TACTACTTCCTCTTCGCC	400	63
		ACTCACCGCCAGCTTCAG
P7	Exon 13 – exon 15	GTCCCCAACCACCTCATCTG	600	63
		TCGTGGAAGAAGGCGGAG
P8	Intron 15 – exon 17	CCTCCGCCTTCTTCCACGA	426	59
		AAAGCAGTCCAACACCCACGA

According to the sequencing results, variations were detected in intron 1-exon 2 and intron 15-exon 17 of DGAT1. Therefore, these two fragments were selected for the following SSCP analysis. The mixture of 2.0 μL PCR products with 8 μL denaturing solution (98% formamide, 10 mM ethylenediaminetetraacetic acid (EDTA), 0.025% bromophenol blue, 0.025% xylene-cyanol) was denatured at 95 °C for 5 min. Subsequently, this mixture was rapidly chilled on ice and then loaded on 16 cm × 18 cm 12% acrylamide: bisacrylamide (37.5:1) gels containing 3.5% v/v glycerol. Samples were electrophoresed in 0.5 × Tris-boric-acid-EDTA (TBE) buffer under electrophoretic conditions of 15 °C, 240 V, 11 h (intron 1 – exon 2) or 23 °C, 250 V, 15 h (intron 15 – exon 17). Gels were stained with silver nitrate using the method of Byun et al. (2009).

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 1 – exon 2 and intron 15 – exon 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: $${\mathrm{Y}}_{\text{ijkl}}=\mu +{\mathrm{ }\mathrm{G}}_{\mathrm{i}}\left({\mathrm{H}}_{\mathrm{i}}\right)\mathrm{ }+{\mathrm{ }\mathrm{P}}_{\mathrm{j}}+{\mathrm{ }\mathrm{S}}_{\mathrm{k}}+{\mathrm{ }\mathrm{A}}_{\mathrm{l}}+{\mathrm{ }\mathrm{e}}_{\text{ijkl}}$$ where Y_ijk is the trait of interest, μ 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.

Results

Identification of variation in yak DGAT1

In 616 Gannan yaks investigated, three distinct genotypes A₁A₁, A₁B₁ and B₁B₁ were detected in intron 1 – exon 2 and a further five distinct genotypes A₂A₂, A₂B₂, A₂C₂, B₂B₂ and C₂C₂ were detected in intron 15 – exon 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).

Figure 1. PCR-SSCP analysis of the DGAT1 gene in Gannan yaks. (A) SSCP banding patterns for two regions of yak DGAT1 gene. (B) Single-nucleotide polymorphism (SNP) detected in intron 1 – exon 2 and intron 15 – exon 17 of DGAT1 gene. (C) Confirmation of single nucleotide variations in two regions using sequence.

Polymorphisms in yak DGAT1

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

Table 2. Genotype and variant frequencies of DGAT1 gene in Gannan yak.

Locus	Genotype frequencies (%)					Variant frequencies (%)
Intron 1 – exon 2	A1A1 (n = 312)		A1B1 (n = 241)		B1B1 (n = 63)	A1	B1
	50.65		39.12		10.23	70.13	29.87
Intron 15 – exon 17	A2A2 (n = 230)	A2B2 (n = 105)	A2C2 (n = 144)	B2B2 (n = 74)	C2C2 (n = 63)	A2	B2	C2
	37.34	17.04	23.38	12.01	10.23	57.55	20.53	21.92

Construction and frequencies of haplotype in yak DGAT1 gene

Six haplotypes were constructed for intron 1 – exon 2 and intron 15 – exon 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.

Table 3. Haplotype construction and frequencies of the DGAT1 gene in Gannan yak.

Haplotypes	Intron 1 – exon 2	Intron 15 – exon 17		Haplotype frequencies (%)
	c.192-117 G > A	c.1252-23 C > T	c.1339 C > T
H1	G	C	C	42.20%
H2	G	C	T	14.45%
H3	G	T	C	17.24%
H4	A	C	C	16.28%
H5	A	C	T	5.68%
H6	A	T	C	4.14%

Association between DGAT1 variants and carcase and meat quality traits

In the single-variant (presence/absence) model (Table 4), the presence of variant A₁ (in intron 1) was found to be associated with decrease WBSF (p = 0.027). The presence of variant B₂ (in exon 17) was associated with an increase in WBSF (p < 0.001), whereas the presence of C₂ (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).

Table 4. Association of DGAT1 variants with carcase and meat quality traits in Gannan yak.

Locus	Traits	Variants	Other variants in model	Numbers		Mean ± SE^a		p Value^b
				Presence	Absence	Presence	Absence
Intron 1 – exon 2	WBSF (kg)	A1	None	553	63	5.22 ± 0.05	5.56 ± 0.14	0.027
		B1	None	304	312	5.32 ± 0.07	5.19 ± 0.06	0.167
		A1	B1	553	63	5.23 ± 0.05	5.63 ± 0.15	0.019
		B1	A1	304	312	5.41 ± 0.08	5.24 ± 0.10	0.135
	DLR (%)	A1	None	553	63	21.59 ± 0.24	22.55 ± 0.70	0.196
		B1	None	304	312	21.80 ± 0.32	21.59 ± 0.32	0.647
	CMR (%)	A1	None	553	63	65.70 ± 0.24	65.46 ± 0.72	0.747
		B1	None	304	312	65.62 ± 0.33	65.74 ± 0.32	0.802
	REA (cm^2)	A1	None	553	63	32.02 ± 0.34	33.34 ± 1.01	0.764
		B1	None	304	312	33.29 ± 0.46	32.82 ± 0.45	0.463
	HCW (kg)	A1	None	269	32	114.07 ± 1.77	116.92 ± 5.17	0.603
		B1	None	150	151	116.17 ± 2.37	112.59 ± 2.36	0.288
Intron 15 – exon 17	WBSF (kg)	A2	None	479	137	5.25 ± 0.05	5.30 ± 0.10	0.656
		B2	None	179	437	5.55 ± 0.08	5.13 ± 0.05	<0.001
		C2	None	207	409	5.04 ± 0.08	5.37 ± 0.06	0.001
		B2	C2	179	437	5.46 ± 0.08	5.13 ± 0.05	0.003
		C2	B2	207	409	5.10 ± 0.08	5.39 ± 0.06	0.008
	DLR (%)	A2	None	479	137	21.77 ± 0.25	21.34 ± 0.46	0.424
		B2	None	179	437	22.05 ± 0.41	21.51 ± 0.26	0.265
		C2	None	207	409	21.63 ± 0.38	21.65 ± 0.27	0.882
	CMR (%)	A2	None	479	137	65.74 ± 0.26	65.47 ± 0.49	0.622
		B2	None	179	437	65.27 ± 0.43	65.85 ± 0.27	0.253
		C2	None	207	409	65.98 ± 0.40	65.52 ± 0.28	0.359
	REA (cm^2)	A2	None	479	137	32.73 ± 0.36	34.18 ± 0.68	0.060
		B2	None	179	437	33.25 ± 0.60	32.97 ± 0.38	0.699
		C2	None	207	409	33.37 ± 0.55	32.89 ± 0.39	0.475
	HCW (kg)	A2	None	230	71	113.95 ± 1.92	115.74 ± 3.46	0.651
		B2	None	99	202	115.63 ± 2.92	113.36 ± 2.04	0.600
		C2	None	90	211	109.02 ± 3.05	116.66 ± 1.98	0.037

WBSF: Warner–Bratzler shear force; DLR: drip loss rate; CMR: cooked meat rate; REA: rib eye area; HCW: hot carcase weight.

^aEstimated marginal means and standard errors (SE).

^bp < 0.05 are in bold.

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 multi-haplotype 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 single-haplotype and multi-haplotype models (Table 6).

Table 5. Association of DGAT1 haplotypes with carcase and meat quality traits in Gannan yak.

Traits	Haplotype	Other haplotypes in model	Numbers		Mean ± SE^a		p Value^b
			Presence	Absence	Presence	Absence
WBSF (kg)	H1	None	326	193	5.18 ± 0.07	5.36 ± 0.08	0.088
	H2	None	114	405	5.34 ± 0.14	5.28 ± 0.10	0.346
	H3	None	141	378	4.99 ± 0.10	5.35 ± 0.06	0.002
	H4	None	141	378	5.27 ± 0.10	5.24 ± 0.06	0.849
	H5	None	50	469	5.87 ± 0.16	5.18 ± 0.05	<0.001
	H1	H3, H5	326	193	5.28 ± 0.13	5.45 ± 0.10	0.077
	H3	H1, H5	141	378	5.27 ± 0.15	5.57 ± 0.09	0.020
	H5	H1, H3	50	469	5.69 ± 0.18	5.15 ± 0.08	0.007
DLR (%)	H1	None	326	193	21.50 ± 0.30	21.90 ± 0.40	0.421
	H2	None	114	405	22.23 ± 0.52	21.49 ± 0.27	0.201
	H3	None	141	378	21.40 ± 0.47	21.75 ± 0.29	0.527
	H4	None	141	378	21.96 ± 0.47	21.54 ± 0.29	0.438
	H5	None	50	469	22.38 ± 0.78	21.57 ± 0.25	0.327
CMR (%)	H1	None	326	193	65.79 ± 0.32	65.37 ± 0.42	0.422
	H2	None	114	405	65.07 ± 0.54	65.79 ± 0.29	0.234
	H3	None	141	378	65.91 ± 0.49	65.53 ± 0.31	0.504
	H4	None	141	378	65.73 ± 0.49	65.60 ± 0.31	0.815
	H5	None	50	469	65.04 ± 0.82	65.80 ± 0.27	0.202
REA (cm^2)	H1	None	326	193	32.52 ± 0.43	33.57 ± 0.56	0.137
	H2	None	114	405	32.97 ± 0.73	32.89 ± 0.39	0.921
	H3	None	141	378	33.22 ± 0.66	32.79 ± 0.41	0.574
	H4	None	141	378	32.65 ± 0.66	33.00 ± 0.41	0.650
	H5	None	50	469	33.34 ± 1.11	32.86 ± 0.36	0.680
HCW (kg)	H1	None	155	97	114.77 ± 2.32	116.60 ± 2.95	0.628
	H2	None	63	189	116.02 ± 3.64	115.29 ± 2.10	0.862
	H3	None	59	193	109.68 ± 3.75	117.24 ± 2.07	0.079
	H4	None	66	186	119.58 ± 3.55	114.01 ± 2.11	0.179
	H5	None	31	221	120.58 ± 5.24	114.76 ± 1.94	0.301
	H3	H4	59	193	111.46 ± 3.75	117.82 ± 2.06	0.163
	H4	H3	66	186	116.44 ± 4.19	112.84 ± 2.26	0.410

WBSF: Warner–Bratzler shear force; DLR: drip loss rate; CMR: cooked meat rate; REA: rib eye area; HCW, hot carcase weight.

^aEstimated marginal means and standard errors (SE).

^bp < 0.05 are in bold.

Table 6. Association of DGAT1 haplotypes copy number with carcase and meat quality traits in Gannan yak.

Traits	Haplotype	Other haplotypes in model	Numbers			Mean ± SE^a			p-value^b
			Absence	One copy present	Two copies present	Absence	One copy present	Two copies present
WBSF (kg)	H1	None	193	214	112	5.36 ± 0.08	5.13 ± 0.08	5.29 ± 0.12	0.106
	H2	None	405	78	36	5.28 ± 0.10	5.31 ± 0.13	5.44 ± 0.16	0.312
	H3	None	378	103	38	5.35 ± 0.06^a	5.09 ± 0.11^a	4.71 ± 0.12^b	0.002
	H4	None	378	113	28	5.24 ± 0.06	5.21 ± 0.10	5.47 ± 0.22	0.564
	H5	None	469	41	9	5.18 ± 0.05^b	5.77 ± 0.18^a	6.32 ± 0.38^a	<0.001
	H1	H3, H5	193	214	112	5.45 ± 0.10	5.22 ± 0.21	5.35 ± 0.23	0.241
	H3	H1, H5	378	103	38	5.57 ± 0.09^a	5.37 ± 0.20^a	4.83 ± 0.28^b	0.006
	H5	H1, H3	469	41	9	5.15 ± 0.08^b	5.63 ± 0.23^ab	6.05 ± 0.41^a	0.026
DLR (%)	H1	None	193	214	112	21.90 ± 0.40	21.68 ± 0.37	21.15 ± 0.50	0.509
	H2	None	405	78	36	21.49 ± 0.27	22.16 ± 0.62	22.39 ± 0.54	0.432
	H3	None	378	103	38	21.75 ± 0.29	21.62 ± 0.54	20.82 ± 0.53	0.612
	H4	None	378	113	28	21.54 ± 0.29	21.82 ± 0.51	22.50 ± 0.67	0.623
	H5	None	469	41	9	21.57 ± 0.25	22.28 ± 0.86	22.83 ± 0.76	0.389
CMR (%)	H1	None	193	214	112	65.37 ± 0.42	65.77 ± 0.39	65.82 ± 0.54	0.723
	H2	None	405	78	36	65.79 ± 0.29	65.45 ± 0.65	64.22 ± 0.45	0.278
	H3	None	378	103	38	65.53 ± 0.31	65.52 ± 0.56	66.96 ± 0.93	0.338
	H4	None	378	113	28	65.60 ± 0.31	65.66 ± 0.54	65.99 ± 0.78	0.938
	H5	None	469	41	9	65.80 ± 0.27	65.05 ± 0.90	64.97 ± 0.92	0.311
REA (cm^2)	H1	None	193	214	112	33.57 ± 0.56	32.58 ± 0.53	32.38 ± 0.53	0.324
	H2	None	405	78	36	32.89 ± 0.39	32.85 ± 0.87	33.22 ± 0.78	0.969
	H3	None	378	103	38	32.79 ± 0.41	33.16 ± 0.76	33.40 ± 0.67	0.842
	H4	None	378	113	28	33.00 ± 0.41	32.73 ± 0.73	32.35 ± 0.46	0.878
	H5	None	469	41	9	32.86 ± 0.36	32.64 ± 1.22	36.52 ± 1.23	0.371
HCW (kg)	H1	None	97	99	56	116.60 ± 2.95	115.43 ± 2.93	113.62 ± 3.81	0.830
	H2	None	189	42	21	115.29 ± 2.10	115.14 ± 4.46	117.76 ± 3.30	0.930
	H3	None	193	41	18	117.24 ± 2.07	111.17 ± 4.50	106.28 ± 2.77	0.179
	H4	None	186	53	13	114.01 ± 2.11	120.62 ± 3.96	115.34 ± 4.99	0.341
	H5	None	221	24	7	114.76 ± 1.94	116.78 ± 5.95	133.41 ± 6.91	0.240

WBSF: Warner–Bratzler shear force; DLR: drip loss rate; CMR: cooked meat rate; REA: rib eye area; HCW: hot carcase weight.

^aEstimated marginal means and standard errors (SE).

^bp < 0.05 are in bold. Means within rows that do not share a superscript letter are different at p < 0.05.

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à 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 (Li et al. 2013). 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 fat-related traits in Iranian sheep. In this study, the presence of variant B₂ 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₁ (in intron 1) and C₂ (in intron 15) were associated with a decrease in WBSF and thus positively correlated with yak meat tenderness. In addition, variant C₂ 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 (An et al. 2018).

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.

Acknowledgements

The authors thank the co-workers of the Gansu Key Laboratory of Herbivorous Animal Biotechnology for the sample and data collection as well as analyses.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This research was funded by the Discipline Construction Fund Project of Gansu Agricultural University [grant number GAU-XKJS-2018-050].