Milk production is an economically important sector of agriculture. Selective breeding for high milk performance is beneficial to maximizing the interests of dairy animal farmers. Over the last several decades, numerous candidate genes and QTLs for milk production traits have been identified in dairy cattle (Fang et al. Reference Fang, Sahana, Su, Yu, Zhang, Lund and Rensen2017). As another major species of milk production, buffalo are lagging behind in research.
Diacylglycerol O-acyltransferase (DGAT1) is a microsomal enzyme which plays an important role in lipogenesis in many tissues, including mammary gland (Cases et al. Reference Cases, Smith, Zheng, Myers, Lear, Sande, Novak, Collins, Welch, Lusis, Erickson and Farese1998). The DGAT1 gene was identified as a strong functional candidate gene affecting milk performance. Defective lactation was observed in DGAT1 −/− female mice (Smith et al. Reference Smith, Cases, Jensen, Chen, Sande, Tow, Sanan, Raber, Eckel and Farese2000), and the polymorphism (K232A) of DGAT1 in exon 8 has been shown to be associated with higher fat percentage in dairy cattle (Naslund et al. Reference Naslund, Fikse, Pielberg and Lunden2008).
Cattle and buffalo species are in a close evolutionary relationship. Interestingly, a fixed K allele in the 8th exon of DGAT1 gene in both Indian water buffalo and Chinese water buffalo was discovered. Therefore, we hypothesized that there may be other SNPs in the coding region of the buffalo DGAT1 gene, which may affect the milk performance of buffalo. In this study, we chose to interrogate the SNPs in the DGAT1 gene coding region and to determine associations between SNPs and milk production traits in buffalo.
Materials and methods
Animals and phenotypic records
Our study used already stocked genomic DNA. A total of 971 buffaloes, including 419 Riverine buffaloes (2n = 50 chromosomes), 370 Swamp buffaloes (2n = 48 chromosomes) and 182 crossbred buffaloes (Riverine buffalo × Swamp buffalo, 2n = 49 chromosomes) were used to identify SNP of the buffalo DGAT1 gene. Riverine buffaloes were collected from Italy, Swamp buffaloes and crossbred buffaloes were collected from China. The phenotypic records of Riverine buffaloes for peak milk yield, total milk yield, protein, and fat percentages were collected from existing databases of the Italian Buffalo Breeders Association (ANASB), including 1294 lactation records from 419 Italian Mediterranean buffaloes with 1–6 parities born in 2000–2011 and located in four farms in Southern Italy. No milk yield records of Swamp buffaloes and crossbred buffaloes existed in the study population. Animal Care and Use Committee approval was not requested because no experimental work was performed.
Single nucleotide polymorphism (SNP) identification and genotyping
Fifty buffalo samples were chosen for detecting the polymorphism of the DGAT1 gene by direct polymerase chain reaction (PCR) sequencing. Twelve pairs of primers (online Supplementary Table S1) were designed for selective amplification of all exons using Primer Premier 5.0 (PREMIER Biosoft, Palo Alto, CA) based on the complete buffalo DGAT1 mRNA sequence (NCBI Reference Sequence: NW_005785136.1). PCR was performed in a Thermal Controller (Bio-Rad, Philadelphia) with a total volume of 30 µl reaction mixtures, containing 2·0 µl pooled genomic DNA template (100 ng), 2·0 µl primers mix (10 µm of each), 25 µl Premix LA Taq (Takara Biotechnology (Dalian) Co., Ltd., Dalian, China), and 11 µl double distilled water. Amplification was carried out with an initial denaturation at 94 °C for 3 min followed by 34 cycles of 94 °C for 30 s, respective annealing temperature for 35 s and 72 °C for 40 s, with a final extension for 8 min at 72 °C, cooling at 4 °C. 1·5% agarose gel was used to evaluate the amplified fragment. The PCR products were used for sequencing (Qing Ke Innovation Biology Technology Co., Ltd., Wuhan, China). Mutations were identified from the aligned sequencing results using the SeqMan software (DNASTAR Inc., Madison, WI).
All samples were genotyped using high resolution melting (HRM) method with the Light Cycler 480 System (Roche, no. 04909631001) according to the operation instructions. Gene Scanning Software was applied for analysis of the HRM data. The genotypes were evaluated based on melting curves of control samples.
Statistical analysis
Allelic frequencies, genotypic frequencies, and Hardy–Weinberg equilibrium was calculated for each SNP locus by the SHEsis online platform. The linkage disequilibrium (D′ and r 2) were estimated using Haploview software version 4.2. Association between DGAT1 polymorphisms for 419 Riverine buffaloes and milk production traits were investigated by fitting the following GLM procedure of SAS9.2 (SAS Inst. Inc., Cary, NC): y ijklm = μ + P i + Y j + S k + G l + e ijklm, where y ijklm is the observation of milk production traits; μ is the overall mean; P i is the fixed effect of the ith parity; Y j is the fixed effect of the jth lactation year; S k is the fixed effect of the kth lactation season-farm classes; G l is the fixed effect of the lth genotype; e ijklm is the random residual. Bonferroni adjusted test was applied to pairwise comparisons among different levels of fixed effects included in the model.
Results and discussion
The buffalo DGAT1 gene comprises 17 exons and the genomic structures are highly conserved with dairy cattle (Yuan et al. Reference Yuan, Zhou, Deng, Hu and Li2007). DNA sequencing analysis showed that two polymorphisms were identified in the coding sequences of DGAT1. The polymorphisms of g.8330T>C and g.9046T>C located in exon13 and exon17 of DGAT1, respectively (Fig. 1). Here, only g.9046T>C was a nonsynonymous switch resulting in the amino acid sequence changed from an arginine (Arg) to a histidine (His) at position 484. The identified SNP were genotyped using HRM method on 419 Riverine buffaloes, 182 crossbred buffaloes and 370 Swamp buffaloes. All the SNP in each studied group were consistent with Hardy–Weinberg equilibrium (P > 0·05). The polymorphisms of g.8330T>C were not found in Swamp buffaloes, and the distribution of g.9046T>C genotypes varied greatly in different buffalo breeds studied here (Fig. 1). The allele frequencies of the two SNPs at the DGAT1 gene in three buffalo breeds are given in online Supplementary Table S2. Mishra et al. (Reference Mishra, Tantia, Kumar and Vijh2007) reported 13 SNPs by comparing the DGAT1exons sequences of the Chinese buffalo (2n = 48) and Indian buffalo (2n = 50), which may due to the difference in the genetic evolution of the DGAT1 in the Riverine buffalo and the swamp buffalo. Moreover, all three buffalo groups shared the same K genotype with K232A mutation in the 8th exon of DGAT1, which is a well-recognized mutation related to milk fat and protein percentages (Spelman et al. Reference Spelman, Ford, McElhinney, Gregory and Snell2002; Naslund et al. Reference Naslund, Fikse, Pielberg and Lunden2008). All buffalo have a fixed K allele which may responsible for high milk fat percentage in buffalo than Bos indicus cattle (Özdil & İlhanh, Reference Özdil and İlhan2012).
Detection of SNP by direct sequencing and the distribution of genotypes in different buffalo groups.

The association between Riverine buffalo milk production traits and genetic variants are given in Table 1. The polymorphisms of g.8330T>C in exon 13 were found to be significantly associated with peak milk yield, milk yield and protein percentage (P < 0·05). Bonferroni test analysis revealed that the C variant was associated with an increase in milk yield and peak but less in protein percentage compared with the T variant. Polymorphisms of g.9046T>C in exon 17 have been previously identified in Chinese water buffalo population, but no association analysis with milk production traits was performed (Yuan et al. Reference Yuan, Zhou, Deng, Hu and Li2007). In our study these polymorphisms were significantly associated with fat percentage, and the buffaloes with TT genotype had a significantly higher fat percentage than those with CC genotype. Moreover, the results of linkage disequilibrium analysis showed that these two candidate loci were in a linkage balance.
Association analysis between genotypes of two SNP and four traits values in Riverine buffalo population.

The different letters of the superscript indicates P < 0·05
It is widely accepted that milk production characteristics are complex traits, more commonly influenced by variations in noncoding regions rather than coding regions (Glazier et al. Reference Glazier, Nadeau and Aitman2002). The findings described here reveal that variations in coding regions influence the variance of complex traits. These mutations may exert effects on associated traits by changing the expression of DGAT1 gene. The mRNA expression analysis of DGAT1 gene in different genotypes of buffalo population could be studied to provide more details about the molecular mechanism underlying the effect of DGAT1 mutations on buffalo milk production traits.
In conclusion, two polymorphisims located in the coding sequence of DGAT1 were identified in Riverine buffalo, Swamp buffalo and crossbred buffalo. Moreover, the SNPs of DGAT1 gene showed significant effects on milk productive traits in Riverine buffalo. The present study provided the evidence that polymorphisms in coding region of DGAT1 gene influence the milk performance of Riverine buffalo.
This work was supported by National Key R&D Program of China (2017YFD0501903) and the Earmarked Fund for Modern Agro-industry Technology Research System (CARS-37-04B).
Competing interests
The authors declare that they have no competing interests.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0022029918000468

