Implication
In this study, we sequenced the coding region of the MC4R gene for 354 yaks, and five single nucleotide polymorphisms (SNPs) were identified, among which SNP4 (1069G>C) and SNP5 (1206 G>C) were missense mutations and resulted in amino acid substitutions (286V>L and 331R>S). SNP4 was assumed to correlate with influential changes in the seventh transmembrane domain of the MC4R protein and with the functional deterioration or even incapacitation of the MC4R gene. Among the three genotypes (GG, GC and CC) of SNP4, CC genotype was associated with increased BW of animals that were 18 months old and with the average daily gain. The data from this study suggested that 1069G>C SNP of the MC4R gene can act as a useful marker for the selection of growth traits in the Maiwa yak breed.
Introduction
BW is an important trait for meat production in livestock animals, and genetic factors play important roles in regulating energy balance. Several genes have been identified to be involved in the neural signaling pathway of energy homeostasis (Barsh et al., Reference Barsh, Farooqi and Rahilly2000; Zhang et al., Reference Zhang, Wang, Chen, Lan, Lei and Fang2009). Among these genes, the MC4R (melanocortin 4 receptor) gene encodes a G-protein-coupled receptor expressed in the appetite-regulating areas of the brain, which takes part in leptin signaling pathways and leads to a decrease in feed intake because of agonist stimulation (Seeley et al., Reference Seeley, Drazen and Clegg2004). Mutations of this gene were found to be the most common cause for hereditary human obesity (Pérusse et al., Reference Pérusse, Rankinen, Zuberi, Chagnon, Weisnagel, Argyropoulos, Walts, Snyder and Bouchard2005). Moreover, significant associations were found between MC4R genotypes and backfat, growth rates and feed intake in some breeds or lines of livestock animals. A functional mutation (Asp298Asn) was associated with daily food intake, growth and fat deposition traits in some commercial pig lines (Hernandez-Sanchez et al., Reference Hernandez-Sanchez, Visscher, Plastow and Haley2003; Houston et al., Reference Houston, Cameron and Rance2004; Van den Maagdenberg et al., Reference Van den Maagdenberg, Stinckens, Claeys, Seynaeve, Clinquart, Georges, Buys and De Smet2007). Variants (−293C>G and −129A>G) in the 5′-untranslated region of the bovine MC4R gene were associated with two growth traits in Nanyang cattle breed (Zhang et al., Reference Zhang, Wang, Chen, Lan, Lei and Fang2009). Therefore, the MC4R gene could be an effective candidate gene for growth-related traits selection in livestock breeding.
Bos grunniens (yak) is the most important domesticated species with outstanding adaptability to the alpine climates of the Qinghai-Tibetan Plateau, and the number of yaks in China account for over 90% of those distributed all over the world (Cai et al., Reference Cai, Mipam, Zhang and Yue2011). As in other animals, the yak MC4R gene is a single-exon gene containing the entire 332 amino acid coding sequence (999 bp) and exhibits abundant variations with other mammalian species (Cai et al., Reference Cai, Mipam, Zhang and Yue2011; Cai and Mipam, Reference Cai and Mipam2013). However, the variation patterns of the MC4R gene in the yak population and their association with growth traits have not been reported thus far. In this study, we aimed to investigate the single nucleotide polymorphisms (SNPs) in the MC4R gene and analyze their association with growth traits in the Maiwa yak population.
Material and methods
Animals and DNA sample
In total, 354 yaks were sampled and each individual was 12 months old (BW: 72±4 kg). The sampled yaks were genetically unrelated and randomly selected from a breeding population of Maiwa yak in Hongyuan county, Sichuan province of China. Most of the yaks included in this study population were selected from different pastures at least 50 km away from each other. These yaks were of similar body size and grazed on alpine meadow grassland from May to October, during which abundant feed grass is available. The BW and average daily gain (from the 12th to the 18th month) for these 18-month-old yaks were determined, and ear tissue samples were collected and maintained in an ice box before being transported to laboratory. Genomic DNA of each individual was extracted from the ear muscle tissue using the normal method of phenol–chloroform extraction (Sambrook et al., Reference Sambrook, Fritsch and Maniatis1989). The sampling and protocol for this research were approved by animal use and ethics committee of the Southwest University of Science and Technology and Southwest University for Nationalities.
MC4R gene sequencing and mutation analysis
PCR was performed using the primers MC4R-F (5′-TGGGACATTTATTCACAGCAG-3′) and MC4R-R (5′-CCTACACAGAAGAAAAAGCT-3′), designed by Cai et al. (Reference Cai, Mipam, Zhang and Yue2011) to amplify a fragment of 1238 bp in length covering the ORF region of the yak MC4R gene. PCR was performed in a 50 μl reaction mixture containing 200 ng of genomic DNA, 1.5 mM MgCl2, 10 mM Tris–HCl (pH 8.3), 50 mM KCl, 0.6 units of Taq polymerase (TaKaRa), 10 μM of each primer and 0.2 mM of each dNTP. Thermal cycling was performed in a PTC-200 thermocycler (MJ Research, Inc., Watertown, MA, USA) under the following conditions: 4 min denaturation at 94°C, followed by 30 cycles of 30 s at 94°C, 30 s at 55°C, 30 s at 72°C and a final 7-min extension at 72°C, before cooling to 12°C for 10 min. PCR products were purified using a Qiagen QIAquick PCR purification kit and were sequenced in an ABI 3730 automated sequencer at BeiJing Zixi Bio Tech Co. Ltd., Changping, Beijing, China. Both strands of the PCR product were completely sequenced.
The MC4R gene sequences of 354 yaks were aligned and edited in Clustal X with parameters set to default (Thompson et al., Reference Thompson, Gibson, Plewniak, Jeanmougin and Higgins1997). Polymorphic sites of nucleotide and amino acid sequences of the coding region of the MC4R gene were explored and analyzed using MEGA 4 (Kumar et al., Reference Kumar, Dudley, Nei and Tamura2008), and Tajima’s D test (Tajima Reference Tajima1989) for neutrality was performed in DnaSP Version 5.0 (Librado and Rozas, Reference Librado and Rozas2009).
SNP genotyping and analysis
The genotypes of each yak from different SNPs of the MC4R gene were determined by examining the sequencing chromatographs. Gene frequency, genotype frequency and deviation from the Hardy–Weinberg equilibrium (HWE) for each SNP site were analyzed using the online program SNPStats (http://bioinfo.iconcologia.net/snpstats/start.htm). This program was also used to measure the pairwise linkage disequilibrium (LD) pattern between SNPs based on multi-allelic D′ (coefficient of LD) and r 2 (correlation coefficient) and estimate dominant haplotypes for all SNPs identified in the gene MC4R.
Statistical association analysis
The associations between MC4R SNP genotypes/major haplotypes and growth traits were analyzed using mixed models in the software SPSS 17.0. The models are as follows:
Model 1:
Model 2:
Model 3:
In these models, Y ijklm or Y ijklmn represents the phenotypic value of the growth trait; μ is the least square mean; Sex i indicates the fixed effect for the i th sex; Age j is the fixed effect for the j th age; Family k indicates the fixed effect for the k th family; Pasture l denotes the fixed effect for the l th pature; Genotype m denotes the fixed effect for the m th genotype; Additive m indicates the fixed effect for the m th additive effect; Dominance n indicates the fixed effect for the n th dominance effect; Haplotype m denotes the fixed effect for the m th haplotype; and e ijklm or e ijklmn is the random error.
Results
Variations and SNP genotypes in the MC4R gene
Yak MC4R is a single-exon gene containing the entire 999 bp sequence coding 332 amino acids (Cai et al., Reference Cai, Mipam, Zhang and Yue2011). In total, 354 MC4R gene sequences of yak were examined, and mutation analysis revealed five SNP sites. Among the five SNPs, SNP1, SNP2 and SNP3 were synonymous C>T, G>T and C>A substitutions at the positions 273, 321 and 864 (AF265221.1; Haegeman et al., Reference Haegeman, Coopman, Jacobs, Mattheeuws, Van Zeveren and Peelman2001), respectively (Table 1). However, SNP4 and SNP5 were both G>C missense mutations at positions 1069 and 1206 (AF265221.1), respectively, and these two sites resulted in amino acid substitutions (V286L and R331S), determined according to yak MC4R (HM051376) (Cai et al., Reference Cai, Mipam, Zhang and Yue2011; Table 1).
SNPs, alleles and corresponding genotypes detected in the MC4R gene of Bos grunniens

SNP=single nucleotide polymorphism.
* The location of each SNP was nominated according to the bovine MC4R gene (GenBank accession number: AF265221.1).
** Amino acid exchanges were determined according to yak MC4R (GenBank accession number: HM051376).
*** Exact test for Hardy–Weinberg equilibrium.
Genotyping the 354 yaks in the five SNP sites revealed that homozygotes were dominant in four sites, which could be displayed by the homozygote frequency as high as 96% (G/G), 85% (C/C), 86% (G/G) and 94% (G/G) in SNP2, SNP3, SNP4 and SNP5, respectively. In SNP1, however, the frequency of the heterozygote (44%) was higher than that of each homozygote (Table 1). The genotypes determined from the sequencing chromatographs of SNP1 and SNP4 are shown in Figure 1. All these SNPs indicated no significant deviation from HWE (P>0.05); therefore, the yak population sampled in this study were in an equilibrium state.
Sequencing chromatographs showing the genotypes determined from SNP1 and SNP4. a=homozygote C/C; b=heterozygote C/T; c=homozygote T/T; d=homozygote G/G; e=heterozygote G/C; and f=homozygote C/C. SNP=single nucleotide polymorphism.
Linkage disequilibrium of the SNPs
Pairwise LD analysis based on D′ (coefficient of LD) and r 2 (correlation coefficient) showed that different degrees of LD were present in the five SNP locations (Figure 2). Higher degrees of LD existed between SNP2 and SNP5, as indicated by the higher r 2 values (r 2=0.81027). Similarly, the higher r 2 values (r 2=0.53816) also illustrated that higher degree of LD existed between SNP4 and SNP5. However, lower extent of LD was discovered between other pairwise SNPs (r 2<0.20346), including SNP1-SNP2, SNP1-SNP3, SNP1-SNP4, SNP1-SNP5, SNP2-SNP3, SNP2-SNP4, SNP3-SNP4 and SNP3-SNP5, which suggested that there might be higher mutation rates or multiple recombination hotspots in the yak MC4R gene.
Linkage disequilibrium (LD) plot for the five single nucleotide polymorphisms (SNPs) of the yak MC4R gene. Color scheme is in accordance with D, D′ and r 2 scheme from SNPStats program. The numerical data in each block donates pairwise value of D, D′ and r 2. The blank block represents the equal value of D, D′ and r 2 between the same pairs of SNPs. The five SNPs are SNP1.273, SNP2.321, SNP3.864, SNP4.1069 and SNP5.1206, in which the number after each SNP represents the nucleotide site in the bovine MC4R gene (GenBank accession number: AF265221.1).
In total, seven haplotypes were identified for all SNPs in the yak MC4R gene (Figure 3), among which TGCGG and CGCGG were dominant with a higher frequency of 46.3% and 44.3%, respectively. However, other types of haplotype were rare as exemplified by low frequencies (<6.64%).
Haplotypes and estimated frequencies for the five single nucleotide polymorphisms (SNPs) of the yak MC4R gene.
Association of SNP genotypes/major haplotypes with growth traits
Statistical association analyses between MC4R genotypes and growth traits of Maiwa yak revealed that SNP1, SNP2, SNP3 and SNP5 showed no effects on BW of animals at 18 months of age and on average daily gain as well. However, missense mutations G>C at position 856 of the SNP4 gene may suggest that the change in MC4R protein and the CC genotype determined from SNP4 were associated with increased BW of animals aged 18 months (CC-113.62±8.63 v. GC-106.87±5.67 and GG-101.95±3.24; P<0.01) and with average daily gain (CC-0.22 v. GC-0.19 and GG-0.17; P<0.05) (Table 2). Two major haplotypes for all SNPs, CGACG and CTCCC were associated with significantly increased BW of animals aged 18 months (P<0.01) and with the average daily gain (P<0.05) (Table 3).
MC4R genotypes corresponding to the least square means, and estimated genetic effects of growth traits in the Maiwa yak population
* Least square means significantly different among the genotypes (P<0.05).
** Least square means significantly different among the genotypes (P<0.01).
Association between the major haplotypes and growth traits in the Maiwa yak population

* Least square means significantly different among the haplotypes (P<0.05).
** Least square means significantly different among the haplotypes (P<0.01).
Analyses of the genetic effects of different SNP genotypes indicated that additive effects were significant or extremely significant at SNP1 (P<0.05) and SNP4 (P<0.01) (Table 2), but there were no significant dominance effects at these loci. The negative dominance effects in SNP4 may indicate negative heterosis. All absolute values of the additive effect were higher than those of the dominance effect at the two loci, which suggested that additive effects played more important roles than dominance effects.
Discussion
Current situation and difficulties in yak breeding
Yaks are one of the most important domesticated species in China; however, studies on comprehensive and accurate assessment of yak genetic resources have been greatly restricted. The lack of reliable DNA markers from yak genome was one of the important reasons, which also resulted in stagnation of marker-assisted selection in yak breeding. On the other hand, yaks are semi-grazing and semi-feeding animals in different periods of the year and it is difficult to record and obtain reliable data of quantitative traits from these animals. Considering these conditions, we selected 12-month-old yaks with similar body size (BW: 72±4 kg) and these animals were grazed on grasslands from May to October, during which abundant feed grass was available.
The missense mutation 1069 G>C in the MC4R gene and its association with growth traits
The missense mutation 1069 G>C in the MC4R gene and its association with growth traits have been studied in different cattle breeds. Five polymorphisms were detected at position 19 (C>A), 20 (A>T), 83(T>C), 128 (G>A) and 1069 (G>C) in the MC4R gene of eight Chinese cattle breeds, among which 1069 (G>C) was significantly associated with backfat thickness (Huang et al., Reference Huang, Gao, Li, Ren, Chen and Xu2010). The missense mutation 1069 G>C in the MC4R gene of Chinese Qinchuan cattle was significantly associated with live weight, carcass weight, backfat thickness and marbling score (Liu et al., Reference Liu, Tian, Zan, Wang and Cui2010). In the MC4R gene of Korean cattle (Hanwoo), 1069 G>C significantly affected backfat thickness (Seong et al., Reference Seong, Suh, Park, Lee and Kong2012). However, no significant association was found between the 1069 G>C locus of MC4R and residual feed intake, feed conversion ratio and average daily gain in a Simmental bull population (Du et al., Reference Du, Chen, Yuan, Zhang, Chen, Wang, Gao, Zhang, Gao, Li and Xu2013).
The yak MC4R is a single-exon gene containing the entire 332 amino acid coding sequence (999 bp; Cai et al., Reference Cai, Mipam, Zhang and Yue2011). Analysis of 354 sequences of the MC4R gene showed five SNP sites, representing 0.5% of the total gene sequence analyzed (999 bp). Therefore, the sequences of MC4R of yaks showed low levels of polymorphisms in this study, which was due to the important function of MC4R in regulating food intake behavior. Low levels of variation in the MC4R gene could also be explained by the view that MC4R might have been subjected to high levels of continuous purifying selection during yak evolution and variations that would influence the activity of metabolism also had been subjected to purifying selection mechanism during yak evolution (Hughes et al., Reference Hughes, Hinney, Brumm, Wermter, Biebermann, Hebebrand and Stoneking2009; Cai et al., Reference Cai, Mipam, Zhang and Yue2011).
Pairwise LD analysis of the five SNPs showed that higher degrees of LD only existed between SNP2 and SNP5 (r 2=0.81027), similarly higher degrees of LD was also examined between SNP4 and SNP5 (r 2=0.56715). However, lower extent of LD was discovered between other pairwise SNPs, including SNP1-SNP2, SNP1-SNP3, SNP1-SNP4, SNP1-SNP5, SNP2-SNP3, SNP2-SNP4, SNP3-SNP4 and SNP3-SNP5. Domestication and breeding improvement were involved in the selection for specific alleles at candidate genes, which resulted in reduced genetic diversity and increased LD relative to unselected genes (Yamsaki et al., Reference Yamasaki, Tenaillon, Bi, Schroeder, Sanchez and Doebley2005; Zhou et al., Reference Zhou, Lei, Rao, Nie, Zeng and Xia2008). In this study, low levels of LD together with a negative Tajima’s D (−1.08498, P>0.10) further indicated low levels of polymorphisms in the yak MC4R gene and yaks may have been subjected to purifying (negative) selection during their evolution. This was also correlated with the fact that yaks have been grazed livestock for thousands of years and less domestication, artificial selection or breeding improvements were imposed on and they have undergone a long history of natural selection during their pasture grazing and extensive management.
Missense mutations caused amino acid substitution in the protein and secondary or spatial structural changes in the protein, which frequently gave rise to functional changes of the protein. Among the five SNPs identified in this study, two (SNP4 and SNP5) were missense mutations and resulted in amino acid substitutions (V286L and R331S), determined according to yak MC4R (HM051376) (Cai et al., Reference Cai, Mipam, Zhang and Yue2011). Association analyses between SNP genotypes and the growth traits indicated that SNP4 may have caused the functional change of MC4R, because CC genotype determined from SNP4, CGACG and CTCCC haplotypes for all SNPs were associated with increased BW of animals aged 18 months and with average daily gain. The secondary structure (http://npsa-pbil.ibcp.fr/cgi-bin/nspa-automat. pl.) and transmembranal region prediction (http://www.cbs.dtu.dk/services/TMHMM/) of the yak MC4R protein revealed that the replacement of valine with leucine at position 286 of the MC4R protein was within the seventh transmembrane domain (Figure 4). In this sense, we assumed that SNP4 resulted in influential changes in the seventh transmembrane domain of the MC4R protein as well as in the functional deterioration or even incapacitation of the MC4R protein, which might contribute to the increased feed intake, BW and average daily gain of the yaks with CC genotypes. Therefore, all these data suggested that 1069 G>C SNP of the MC4R gene could be used in marker-assisted selection of growth traits in Maiwa yak breed.
Prediction of the secondary structure and transmembrane region of the yak melanocortin 4 receptor (MC4R) protein. a=the secondary structure prediction of the yak MC4R gene; the blue, red and purple vertical lines stand for α-helix, extended chains and random coils, respectively. b=the transmembrane region prediction for the yak MC4R; the red, blue and purple areas stand for transmembrane regions, inside regions of membrane and outside regions of membrane, respectively.
Acknowledgments
This work was supported by the grants from China Postdoctoral Science Foundation funded project (No.20090461336) and China Postdoctoral Special Science Foundation funded project (No.201003703).

