Association study and expression analysis of GPC1 gene copy number variation in Chinese Datong yak (Bos grunniens) breed

Abstract Copy number variations (CNVs) can cause radical changes in phenotype variation, gene expression and evolutionary adaptation, through gene dosage effects, disruption of transcript arrangement and regulatory polymorphisms. Therefore, the current study was undertaken to investigate the distribution of glypican1 (GPC1) CNVs in five domestic yak breeds, and their association growth traits and gene expression. The data were analysed using real-time quantitative PCR (qPCR). Loss and normal copy number types had a significant (p < .05) associated with body height, weight and length and chest girth than gain of copy number types. Population studies indicated that loss of copy number was more frequently observed in the Tianzhu, Gannan and Plateau yaks than the Polled and Datong yaks. We elucidated a moderately negative and highly significant (p < .005) correlation among DNA copy numbers and mRNA transcription levels of GPC1. From these results, we hypothesise that GPC1 CNVs could alter mRNA transcription levels of skeletal muscles which impact on quantitative of growth traits. Our study makes data available for breeders and genomics studies focussing on the fundamental function of GPC1 CNVs in gene expression and growth traits and presents the foundation for future molecular marker application in designing yak breeding programmes. Highlights The GPC1 gene encompasses within CNVRs of yak individuals that overlap with growth traits. The copy number types of GPC1 had a significant associated with growth traits of yak breeds Negative and highly significant correlation among DNA copy numbers and mRNA transcription levels of GPC1 are obtained


Introduction
The yak (Bos grunniens) is one of the world's most significant domestic animals and survives in the extremely harsh environment of the Qinghai-Tibetan Plateau 'roof of the world' (Wiener et al. 2003). Long hair and large horns are particular phenotypic characteristics of the yak, which is an iconic symbol of animal husbandry at the Qinghai-Tibetan Plateau . Indeed, the yak has become strongly integrated into Tibetans' socio-cultural life (Liang et al. 2016) and occupies a significant economic role in the mountainous regions of Asia (Medugorac et al. 2017); they provide milk, meat, transportation ('boat of the Plateau'), hair, draught power and fuel for dwellers in high-altitude environments (Wiener et al. 2003). Notably, the publication of genome assembly version 1.1 (Yak 1.1.) of Bos grunniens, which is comparable with the cattle genome (UMD 3.1), is significantly important for envisaging gene and protein function and providing genomic sequence data (Hu et al. 2012). In previous studies, genetic variation in livestock has been described using microsatellites (Curi et al. 2005) and single-nucleotide polymorphisms (SNPs) (Utsunomiya et al. 2013). Until recently, SNPs were considered nucleotide variations detected during genome sequencing, and the most important cause of genetic variation by the livestock research community. For this reason, many studies have utilised SNP markers for genetic improvement of beef and dairy cattle (Hayes et al. 2009;VanRaden et al. 2009;Wiggans et al. 2009). The development of yak genespecific SNPs for adaptation mechanisms to hypoxia and growth, and detection of the causal mutation of polledness in yak breeds have been reported Wu et al. 2015;Lan et al. 2018). However,  reported that CNVs are significant sources of genetic variation and play a major role in phenotypic diversity and evolutionary adaptation (Zhang et al. 2009;Hou et al. 2011); variation in the number of gene copies in the genome (Redon et al. 2006;Zhang et al. 2009) affects the expression of genes as well as having a large-scale influence on the transcriptome ). Gene dosage, gene interruption, generation of a fusion gene, position effects, unmasking of recessive coding region mutations, and other functional SNPs are predominant molecular mechanisms by which CNVs convey diverse phenotypes and complex diseases (Lupski and Stankiewicz 2005;Redon et al. 2006;Stranger et al. 2007;Zhang et al. 2009). Stranger et al. (2007) reported that 83.6% and 17.7% of variation in gene expression was elucidated by SNPs and CNVs. According to Gamazon and Stranger (2015), these findings explained the underestimation of the impact of CNVs on the transcriptome because of the relatively greater completeness of SNP genotyping and the difficulty of genotyping CNVs due to their size. In addition, de novo structural changes of gene frequency could be greater in magnitude than the rate of change of specific nucleotides (Lupski 2007;Zhang et al. 2009). This implies that CNVs are significant sources of various in genomic able to explain the missing heritability when compared with SNP variants. Genomic copy number variations stretching from 1 kb to several Mb contain gains and losses compared with the designated reference genome sequence (Scherer et al. 2007). Furthermore, in past years, genomic detection technologies have provided advanced knowledge for construction of maps of CNVs for cattle (Liu et al. 2010), horses (Ghosh et al. 2014), goats (Fontanesi et al. 2010), sheep (Liu et al. 2013), pigs (Paudel et al. 2013), dogs (Alvarez and Akey 2012), chickens (Zhang, Du, et al. 2014), and yaks (Zhang et al. 2016), using different methods including whole-genome sequencing, array-based comparative genomic hybridisation (aCGH), single nucleotide polymorphism (SNP) array, and next-generation sequencing (NGS) and providing a very significant resource for determining how CNVs convey economically important phenotypic traits in livestock species.
In mammals, six classes of the glypican family have been identified; however, only glypican-1 (GPC1) (Campos et al. 1993), which is an extracellular heparin sulphate proteoglycan (Masuda et al. 2010;Harthan et al. 2013) is found in skeletal muscle. Velleman et al. (2006) reported alteration of proliferation and differentiation by the GPC1 gene through fibroblast growth factor 2 (FGF2). Similarly, Song et al. (2010) suggested that GPC1 plays a central role in cell proliferation, differentiation, and responsiveness to FGF2 in turkey myogenic satellite cells. Therefore, the above studies indicated that the GPC1 gene has a critical role among proteoglycans to differentially regulate muscle cell proliferation, differentiation, and cellular responsiveness to FGF2 (Brandan et al. 1996;Velleman et al. 2004Velleman et al. , 2006Velleman et al. , 2007Velleman et al. , 2008Velleman et al. , 2012. Interestingly, the CNVs of the GPC1 gene are related to meat production and quality, economically important traits that have been extensively considered for artificial selection in yak breeding (Zhang et al. 2016).
No previous studies have reported on the role of CNVs of the GPC1 gene in the domestic yak. Zhang et al. (2016) confirmed that the yak GPC1 gene was found within a CNV region (CNVR) that encompasses quantitative trait loci (QTLs) influence muscle development and may affect growth traits through modulation of gene expression. Thus, the current study was undertaken to investigate the distribution of GPC1 CNVs in five domestic yak breeds and their relationship with growth traits and gene expression.
The animals were in good physical condition and not genetically related individuals. All yaks were grazed on natural pasture without feeding supplementation, which provided similar feeding conditions and management. Data on the growth traits of body weight, body height, body length, chest girth, and cannon width were collected form 536 individual of Datong yaks at 6 months (n ¼ 353) and 5 years (n ¼ 183) of age for related studies. Body measurements were taken according to the method of Gilbert et al. (1993).

Data collection and isolation of ribonucleic acid (RNA) and deoxyribonucleic acid (DNA)
Fresh blood samples were collected from a total of 656 yaks (536 Datong yak, 30 Polled yaks, 30 Tianzhu yak, 30 Gannan yaks, and 30 Plateau yak) from the jugular vein into vacuum tubes (5 mL). From these 30 Datong, 30 Polled yaks, 30 Tianzhu yak, 30 Gannan yaks, and 30 Plateau yak were used for study of distribution of population. Tissue samples from skeletal muscle, heart, liver, lung, brain, spleen, kidney, and adipose fat (no adipose for foetal samples) were collected from three 90-day-old foetal (embryo) and three 3-year-old adult Datong yaks for total RNA isolation. In addition, 41 (n ¼ 41) skeletal muscles were collected for gene expression and CNV studies from the multiple individuals of adult Datong yaks. All blood and tissue samples were stored at À80 C before genomic DNA and total RNA purification.
Genomic DNA (gDNA) from blood samples was purified using a Clot Blood DNA kit (CWBIO, China, Beijing), and genomic DNA longer than 10 kb was obtained using a E.Z.N.A. MicroElute DNA Clean-Up kit (Omega Bio-tek, Norwalk, GA) following the manufacturer's instructions. Genomic DNA from skeletal muscles was extracted with a DNeasy V R Blood and Tissue kit (Qiagen V R ) following the manufacturer's recommended standard procedure. Total RNA from the tissues was also isolated using a TRIzol reagent (TriPure Isolation Reagent, Roche, carlsbad, CA, USA) and an RNeasy V R Blood and Tissue kit (Qiagen V R ), according to the manufacturers' instructions. The concentration and quality of RNA and DNA were examined using a NanoDrop TM BioPhotometer 2000 (Thermo Fisher Scientific, Inc., Waltham, MA) and through electrophoresis on ethidium bromide-stained 1% agarose gels. Complementary DNA (cDNA) was produced by reverse transcription from 1 lg RNA using a PrimeScript RT Reagent Kit with genomic DNA Eraser (TaKaRa Bio Inc., Shiga, Japan).

Primer design and PCR amplification
According to available bovine gene sequences, a pair of primers was designed for analysis of CNV and gene expression of the GPC1 gene using the National centre for Biotechnology information (NCBI) primer-BLAST webtool (https://www.ncbi.nlm.nih.gov/tools/primerblast/index.cgi?LINK_LOC¼blastHome) (Table1). A polymerase chain reaction (PCR) mixture containing 50 ng gDNA or cDNA, GoTaq V R Green master mix 2X, 10 lM primer, and nuclease-free water (ddH2O) was used for checking the amplification primers in a total volume of 25 lL (Promega, Madison, WI). Thermal cycling procedures were applied as follows: an initial denaturation at 95 C for 2 min, followed by 35 cycles of extension at 72 C for 5 min and holding tubes at 4 C forever. The PCR products were loaded directly into 1. 0% agarose gel because GoTaq V R Green Master Mix reaction buffer contains dyes with sufficient density to sink in the wells of agarose or non-denaturing TBE polyacrylamide gels. The amplification primers were also checked by melting curve analysis of qPCR.

Copy number variation and expression of the GPC1 gene
Primers designed for copy number variation and gene expression analyses are illustrated in Table 1. The relative expression and CNVs of the bovine GPC1 gene were considered in this study using quantitative realtime polymerase chain reaction (qPCR). Neither CNVs nor segmental duplication of bovine basic transcription factor 3 (BTF3) genes was found among variants in genomic databases (Bickhart et al. 2012); therefore, BTF3 was used as a diploid internal reference gene for genomic qPCR. However, for gene expression studies, the bovine glyceraldehyde-3-phosphate dehydrogenase gene (GAPDH) was chosen as a reference gene because GAPDH is commonly used by biological researchers as a control for qPCR; it is considered a housekeeping gene owing to its stable expression in most tissues and cells (Butte et al. 2001). The relative mRNA expression and CNV of the GPC1 were investigated using the Bio-Rad CFX 96 TM Real-Time Detection System (Bio-Rad, Hercules, CA). The qPCR was performed in a 25 lL total reaction mixture containing 50 ng of gDNA or cDNA, 12.5 lL SYBR V R Premix Ex Taq TM II(2X) (Tli RNase H Plus), and 10 pmol of primers (TaKaRa Bio, Inc., Shiga, Japan). The thermal cycling profile for the qPCR experiment included one cycle of 95 C for 1 min followed by 39 cycles of denaturation at 95 C for 10s, annealing at 60 C for 30s, and extension at 68 C for 10s. Melting curve analysis and notemplate control reactions were used for checking primers. All experiments were replicated three times, and the mean value and standard deviation were used for statistical analysis.

Statistical analysis
Threshold cycle (DDC t ) values were used to investigate relative copy number (Yim et al. 2011). Copy number of DNA was calculated using the average threshold cycle (DCt) value for the three replicates normalised against the reference gene, BTF3, by subtracting the BTF3 Ct value from the sample Ct value for each replicate; the final copy number value was determined using the formula 2 Â 2 -DDCt (Bae et al. 2010), where a copy number of 2 (diploid) was considered the normal DNA copy number (Yim et al. 2011). The C t values were converted to the nearest integer as described elsewhere Xu, Shi, et al. 2014;Liu et al. 2016;Shi et al. 2016;Zhou et al. 2016). The quantitative mRNA expression level of the target gene was determined by using the threshold cycle 2 -᭝᭝Ct method (Schmittgen and Livak 2008).
The association between CNV types of the GPC1 gene and growth traits in Datong yaks was determined using analysis of variance (ANOVA) with SPSS software (IBM SPSS 20,Hausmalt,CH,Switzerland). In the biostatistical model, the copy number types were grouped into loss (copy number 0 or 1), gain (copy number > 2), and normal (copy number ¼ 2) relative to the reference gene BTF3, as illustrated elsewhere Xu, Shi, et al. 2014;Liu et al. 2016;Shi et al. 2016;Zhou et al. 2016). The following analytical models were fitted for growth traits: where Y ij is the observation of the growth traits; l is the overall mean of each growth traits, CNV i is represents the fixed effects of the i th CNV type (i ¼ 1, 2, 3) of GPC1 gene, Aj is the fixed effect of the j th age (j ¼ 1, 2), and eij is the random residual error. Moreover, the Bonferroni correction, which presents a modified significance criterion (p/m, where m is the overall number of independent statistical tests conducted on the given data and p is the significance level of .05), was utilised to avoid type I errors derived from multiple comparisons .
The diversity between yak breeds was determined relative to the relative to the reference gene BTF3, using log 2 ratio (log 2 2-᭝᭝Ct ) values; the copy number types were described as gain (>0.5), loss (<À0.5) and normal (<j±0.5j) Xu, Shi, et al. 2014;Liu et al. 2016;Shi et al. 2016;Zhou et al. 2016). Pairwise comparison and scatter plots of yak breeds were carried out by using GraphPad Prism version 5.00. Pearson Product Moment Correlation in R 2.15.0 software was used to determine the relationship between the log 2 ratio of the DNA copy number and the messenger RNA (mRNA) GPC1 expression levels (Knezevic et al. 2007). According to Stranger et al. (2007), correlations generated using log 2 ratio signals give strong r 2 and p values; therefore, the correlation study was performed directly using data obtained from log 2 ratios as mentioned previously ).

Distribution of GPC1 CNV in populations
As shown in Figure 1, the GPC1 gene in the Tianzhu, Plateau, and Gannan yak individuals exhibited more copy number loss; however, the Datong and polled yak individuals exhibited gain copy numbers. As shown in Figure 1(B), across all breeds, there was an overall frequency of GPC1 copy number (Figure 1(B)) loss (copy number 0 or 1), normal level (copy number 2), and gain (more than two copies) of 55.3%, 26.7% and 18%, respectively (Supplementary Table S2 Table S1). The copy number of the GPC1 gene ranged from zero to seven and one to five copies in Tianzhu, Datong and Polled yaks, respectively, indicating that the consequence of CNVs of the GPC1 gene may be high genetic variability within a population of these breeds (Figure 1(B)).

Association between GPC1 CNV and growth
A one-way ANOVA was carried out to determine the relationship between GPC1 CNV and growth traits in Datong yaks. As elucidated in Table 2, there was a remarkable association between growth traits and GPC1 CNV. These results showed that the loss and normal copy number types of the GPC1 gene is significantly associated (p < .05) with higher body length, weight and height and chest girth than the gain copy number type in six-month-old yaks. On the other hand, the loss copy number types of the GPC1 gene is highly associated (p < .05) with body weight, length and height and chest girth than gain copy number type in five-year-old yaks. Therefore, these results suggest that GPC1 CNV could be substantially related to growth traits of Datong yaks.

Expression pattern analysis of the GPC1 gene in different tissues
During foetal periods, GPC1 had noticeably high expression in muscle and spleen, while moderate expression was observed in brain and lung, and weak expression in liver, kidney and heart tissues (Figure 2(A)). However, significantly (p < .05) high mRNA expression levels of GPC1 were found in adipose tissues of adult yaks (Figure 2(B)). Also, moderate mRNA expression levels were identified in the spleen, muscle and brain, while low mRNA expression levels were observed in kidney, lung, heart and liver tissues (Figure 2(B)). Moreover, a significant increase in GPC1 mRNA levels was determined in the muscle, spleen, and brain between foetal and adult yak samples (p < .05) (Figure 2(C)).

Correlation analysis of CNV and mRNA expression of the GPC1 gene
The genomic functions of the GPC1 CNVs and mRNA expression levels were examined in adult skeletal  muscle based on the scale of expression. Among 41 skeletal muscles, we identified nine normal copy number types, 25 gain types, and seven loss types. The relative copy numbers and mRNA expression levels varied from one to seven copies (Figure 3(A)) and 0.12-to 2.75-fold (Figure 3(B)), respectively. In the meantime, a moderately negative and highly significant (p ¼ .005) correlation (r ¼ À0.43) was detected between CNV and mRNA expression levels ( Figure  3(C)). On the other hand, in Figure 3(D), when comparing copy number types and gene expression, normal and loss copy number was significantly (p < .05) The impacts of DNA CNVs on GPC1 mRNA gene expression level (!1, n ¼ 7; 2, n ¼ 9; 3, n ¼ 25) were evaluated between different CNV types. The error bar showed the standard deviation (SD). Different letters indicate significant difference (p < .05). CNV: copy number variation.
associated with higher mRNA expression levels than that of gain types, while no difference was observed between loss and normal copy number types.

Discussion
In yak, CNV surveying has been restricted to the discovery of CNV regions (CNVRs) using aCGH (Zhang, Jia, et al. 2014), BovineHD Genotyping Bead Chip Array (Quanwei et al. 2014), and NGS (Zhang et al. 2016). However, no efforts have been made to study the association of CNV with gene expression in different breeds of Datong yak. Indeed, Xu et al. (2016) indicated that the biological impact of CNVs was a dynamic research area, due to their diversity, origin, and genetic properties. The qPCR detection method was highly powerful for CNV investigation (Liu et al. 2010;Cole et al. 2011;Bickhart et al. 2016;Zhang et al. 2016). The log 2 ratio indicated that greater copy number loss was observed in the GPC1 CNVs in the Tianzhu, Gannan and plateau yak populations, whereas copy number gain was observed in the Datong and polled yaks. Population differentiation in this study was discovered using the log 2 ratio to differentiate the diversity between yak breeds, and this method is more sensitive to detection of homozygous deletions than duplications (Olson 1996;Liu et al. 2016). On the other hand, gene loss was a major pattern of molecular evolution (Ko et al. 2012). Correspondingly, Zarrei et al. (2015) found that CNVs caused by deletions and duplications in the genome were the predominant causes of the variability among populations, with duplications and deletions occurring as results of positive and negative selective pressure in the individual populations, respectively. These results further explain the frequencies of different copy-number types in five breeds with variations that could originate from the diversity of breeds kept in different environments (Lehnert et al. 2007;Yang et al. 2017) and selection Yang et al. 2017). The present results show that the copy number gain events might be due to extra genetic material (artificial insemination, hybridisation and introgression) obtained during breed formation in Datong and Polled yak, since Datong yaks are produced by mating Huanhu yak cows with wild yaks (male) for meat production and Polled yak are produced from breeding Polled yak cows with horned yak bulls (Wiener et al. 2003;Liu et al. 2014); the other three yak breeds belong to the indigenous domestic yaks of China. These studies revealed that the CNVs of the GPC1 gene may affect the phenotypes of different breeds, hinting at discrepancy in genomic CNVs. This study was based on real-time polymerase chain reaction (qPCR) using five different domestic yaks; therefore, it is hard to judge the CNVs of genes in other studies due to the different populations, sizes, platforms, and algorithms applied for surveying.
Our results showed that the normal and loss of CNV types of the GPC1 gene are significantly (p < .05) associated with growth traits. This study may be revealed that GPC1 gene resides in the CNVR that influence on the quantitative traits of growth of yaks (Zhang et al. 2016). These results were confirmed by  (Xu, Shi, et al. 2014), KCNJ12 , MAPK10 , CYP4A11 (Yang et al. 2017), and MICAL-L2  genes and growth traits of Chinese cattle have been reported. Previous work have also revealed that most CNV genes are located within CNVRs that overlap with QTLs (Liu et al. 2010). Moreover, work revealed by Stothard et al. (2011) found that CNVs of the PLA2G2D gene was associated with milk production, health, and meat quality in beef and dairy cattle. Consistently, Cole et al. (2011) revealed that genes and chromosome regions of CNVs were associated with body shape traits. Comparable to these results, the CNV of PLA2G2D gene was associated with growth traits (Zhang, Jia, et al. 2014). Furthermore, Da Silva et al. (2016) reported that CNVs of the PLAG1 zinc finger (PLAG1) affect height. Therefore, this study was very interesting, because it clearly shows the association between CNVs and growth traits in the Datong yak breed; however, a larger yak population may be required for further elaboration.
The current study revealed that the CNVs of the GPC1 gene affect growth traits via altering mRNA transcription levels. Besides, GPC1 mRNA transcription levels were significantly (p < .05) high in muscle, spleen and brain among foetal and adult yaks (Figure 2(C)). However, the comparison showed that GPC1 mRNA transcription levels were significantly (p < .05) increased in muscle, spleen, brain and lung of foetal yak. Low mRNA expression levels were recorded in liver and heart tissues; this may indicate that GPC1 expression is highly regulated in liver and heart. In accordance with this result, a previous study indicated that the expression of GPC1 decreased with age, being increased in seven-week-old cells and decreased in 16-week-old cells compared with 1-day-old cells at 24 and 48 hours of differentiation in turkey (Harthan et al. 2013). Interestingly, the expression of GPC1 was age-dependant during proliferation and differentiation (Harthan et al. 2013), the origin of muscles (Zammit et al. 2002), genetics ) and tissues (Guryev et al. 2008). The high expression of the GPC1 mRNA in the muscles might explain the important of the GPC1 gene for muscle proliferation and differentiation in embryonic stages of yak ). on the other hand, the high expression of the GPC1 mRNA in the in the brain is essential for development of brain physiology, because loss of the GPC1 gene predominantly affects brain size and patterns at birth (Jen et al. 2009). In accordance with previous observations of Shiau et al. (2010) GPC1 mRNA transcription levels were significantly high in foetal Datong yaks, may be indicating the possible roles of GPC1 in the embryonic developments (Shiau et al. 2010). Comparable to these observations, high levels of GPC1 expression were found throughout the brain, kidney and skeletal system in mice (Litwack et al. 1998) and also a high expression of GPC1 in the skeletal muscle (Velleman and Song 2017).
Biologically, GPC1 mRNA expression levels in adult skeletal muscle were unquestionably influenced by CNVs. Meanwhile, a moderate negative correlation (r ¼ -0.42) was reported between DNA copy number and mRNA expressions of the GPC1 gene, and this association was highly significant (p ¼ .005). Consistent with these results, a negative correlation were reported between the levels of DNA copy number and MICAL-L2 gene in cattle , DNA copy number and the LEPR gene in cattle (Shi et al. 2016) and DNA copy numbers and the RHACD8 gene in chickens (Wright et al. 2009). However, a positive correlation was observed between the levels of DNA copy number and the MYH3 (Xu, Shi, et al. 2014) and MAPK10 ) genes in cattle. Furthermore, Xu et al. (2017) reported a positive association between the DNA copy number and the LEPR gene in the skeletal muscle and adipose tissues of Chinese cattle (Henrichsen, Vinckenbosch, et al. 2009). These results suggest that the mRNA expression and DNA copy numbers of GPC1 are negatively correlated due to gene dosage effects. Consistent with this study, Conrad et al. (2010) and Feuk et al. (2006) reported that CNV affects phenotypic variation directly through mechanisms of gene dosage and/or altering gene expression in the genomic regions. According to a study by Gamazon and Stranger (2015), expression of dosage-sensitive genes could be modulated by gene duplication and deletion. In line with this result, Merla et al. (2006) reported that gene transcription in tissues influenced by alteration of gene dosage and effects of neighbouring genomic regions (Stranger et al. 2007;Orozco and Cokus 2009). Correspondingly, Liu et al. (2016) suggested that CNVs of MAPK10 affect gene expression through gene dosage and gene-gene and gene-environment interactions. In contrast, McCarroll and Altshuler (2007) observed that mRNA expression levels affected by additional regulators or transcription factors that regulate gene expression in trans, and DNA methylation. Consistent with this finding, CNV is positively and/or negatively correlated with the mRNA expression landscape (McCarroll and Altshuler 2007;Guryev et al. 2008). Therefore, this observation revealed that GPC1 CNV may be a crucial factor accountable for the growth of yak breeds. The results of the present study showed that CNVs of GPC1 gene is associated with growth traits in the Datong yak breed.

Conclusions
The results elucidated that DNA copy number and mRNA expression were remarkably negatively correlated, showing the potential impact of GPC1 CNVs on growth traits in Datong yak breed. During the foetal stage, GPC1 genes were largely expressed in spleen, muscles, and brain tissues, underscoring their fundamental role in skeletal muscle development, muscle cell proliferation, differentiation and cellular responsiveness. This study reports on the crucial function of GPC1 CNVs in Chinese Datong yak breed, and it is hypothesis that a selection of the young animals on the base of their CNVs of the GPC1 gene important resources for future development of breeding programmes in yak.