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

Keywords:

• Copy number types
• growth traits
• mRNA transcription

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 (Liu et al. 2014). 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 gene-specific SNPs for adaptation mechanisms to hypoxia and growth, and detection of the causal mutation of polledness in yak breeds have been reported (Liu et al. 2014; Wu et al. 2015; Lan et al. 2018). However, Henrichsen, Chaignat, et al. (2009) 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 (Henrichsen, Chaignat, et al. 2009). 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. A considerable number of studies have detected CNVs in genomic regions using different platforms; however, few or none have related this information to association and expression analysis. In the bovine genome, CNVs cover roughly 4.6% (Hou et al. 2011; Bagnato et al. 2015; Da Silva et al. 2016; Keel et al. 2016; Prinsen et al. 2016). In the human genome, CNVs contribute approximately 4.8–9.5% of variations, including the creation of evident phenotypic effects from the deletion of approximately 100 genes (Zarrei et al. 2015). Strikingly, in cattle, milk production is associated with CNVs (Xu, Cole, et al. 2014; Prinsen et al. 2017). Further studies performed in Chinese cattle revealed that CNVs of the leptin receptor (LEPR) (Shi et al. 2016), myosin heavy chain 3 (MYH3) (Xu, Shi, et al. 2014), potassium voltage-gated channel subfamily J member (KCNJ12) (Zhou et al. 2016), mitogen-activated protein kinase 10 (MAPK10) (Liu et al. 2016), MICAL-like 2 (MICAL-L2) (Xu et al. 2013), and phospholipase A2 group IID (PLA2G2D) (Zhang, Jia, et al. 2014) genes are associated with growth traits. In addition, Zhang et al. (2016) underlined that CNVs are associated with genetic diversity between individuals (Stankiewicz and Lupski 2010) and susceptibility to a complex human disease (Yang et al. 2007).

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. 2004, 2006, 2007, 2008, 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.

Materials and methods

Study location and body measurements

Five representative Chinese domestic yak breeds were studied: Polled yak (Qinghai province, Datong Yak Farm in Qinghai, China, elevation 3200 m), Tianzhu White yak (Gansu province, Tianzhu Tibetan autonomous county, China, elevation 3000 m), Plateau Yak (Qinghai province, northern and southern Qinghai, China, elevation 3700–4700 m), Datong Yak (Qinghai province, Datong Yak Farm in Qinghai, China, elevation 3200 m), Gannan yak (Gansu province, bordering Sichuan and Qinghai, China, elevation 3300–4400 m).

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^® Blood and Tissue kit (Qiagen^®) 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^® Blood and Tissue kit (Qiagen^®), according to the manufacturers’ instructions. The concentration and quality of RNA and DNA were examined using a NanoDrop™ 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 μg 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/primer-blast/index.cgi?LINK_LOC=blastHome) (Table1). A polymerase chain reaction (PCR) mixture containing 50 ng gDNA or cDNA, GoTaq^® Green master mix 2X, 10 μM primer, and nuclease-free water (ddH2O) was used for checking the amplification primers in a total volume of 25 μL (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^® 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 real-time 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™ Real-Time Detection System (Bio-Rad, Hercules, CA). The qPCR was performed in a 25 μL total reaction mixture containing 50 ng of gDNA or cDNA, 12.5 μL SYBR^® 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 no-template 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.

Table 1. Primers designed for qPCR based on bovine genome sequences.

Level	Accession no.	Gene name	Primers Sequence (5′->3′)	Product length, bp	Annealing temperature, °C
DNA	XM_014480020.1	GPC1	F: GCTTTGCTTTTCCTCACGGG	86	57.45
			R: ACATGGTTGGCGTGTTTGTG		55.40
	NM_001025314.2	BTF3	F:AACCAGGAGAAACTCGCCAA	166	55.40
			R: TTCGGTGAAATGCCCTCTCG		57.45
mRNA	XM_015462903.1	GPC1	F: CTGCAGAAGCTGGTCTCCGAA	90	59.76
			R: CTTCTCACTGCACAGCGTTCC		59.76
	NM_001034034.2	GAPDH	F: ATGAAAGGGCCATCACCATC	204	55.85
			R: GTGGTTCACGCCCATCACA		60.00

Statistical analysis

Threshold cycle (ΔΔC_t) values were used to investigate relative copy number (Yim et al. 2011). Copy number of DNA was calculated using the average threshold cycle (ΔCt) 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^–ΔΔCt (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 et al. 2013; 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 et al. 2013; 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: $${\mathrm{Y}}_{\text{ij}}=\mu +{\text{CNV}}_{\mathrm{i}}+\text{Aj}+{\epsilon }_{\text{ij}}$$ where Y_ij is the observation of the growth traits; μ 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 εij 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 (Xu et al. 2013).

The diversity between yak breeds was determined relative to the relative to the reference gene BTF3, using log₂ ratio (log₂^2-△△Ct) values; the copy number types were described as gain (>0.5), loss (<−0.5) and normal (<|±0.5|) (Xu et al. 2013; 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₂ 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₂ ratio signals give strong r² and p values; therefore, the correlation study was performed directly using data obtained from log₂ ratios as mentioned previously (Liu et al. 2016).

Results

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). Among them, about 80% of Tianzhu, 80% of Plateau and 97.7% of Gannan yaks exhibited copy number loss while 53% of Datong and 50% of Polled yaks contained normal copy number types of in the CNV of the GPC1 gene (Supplementary Table S2). Significant differences were found in the pairwise comparisons of Polled vs. Datong, Tianzhu, Gannan, and Plateau yak, and Datong vs. Tianzhu, Gannan, and Plateau yak (p < .05, Supplementary Table S1); low significant differences were found for Tianzhu vs. Gannan yak breeds, while there was no significant difference for Datong vs. Polled, Plateau vs. Tianzhu and Gannan, and Gannan vs. Plateau yak breeds (Supplementary Table S1).

Figure 1. Distributions of the CNV of the GPC1 gene in five yak breeds. (A) Scatter plot of CNV of the GPC1 gene in five yak breeds (n = 30). (B) The frequency of copy numbers (0 and 1 = loss, 2= normal and >2= gain copy number types) in five yak populations. CNV: copy number variation.

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.

Table 2. Analysis of the association between CNVs in the GPC1 gene and growth traits in Datong yak breeds.

Age	Growth trait	CNV type (Mean ± SE)			p-value
		Loss (n = 100)	Normal (n = 120)	Gain (n = 133)
6 months (n = 353)	Body height, cm	91.6 ± 0.40^ab	92.4 ± 0.30^a	90.6 ± 0.40^b	.0020
	Body length, cm	90.0 ± 0.30^a	89.4 ± 0.30^b	88.3 ± 0.40^b	.0030
	Chest girth, cm	120.8 ± 0.50^a	120.2 ± 0.60^a	118.3 ± 0.60^b	.0090
	Cannon width, cm	12.1 ± 0.03	12.1 ± 0.03	12.1 ± 0.03	.7700
	Body weight, kg	88.7 ± 1.10^b	95.7 ± 0.80^a	87.2 ± 1.10^b	.0001
		Loss (n = 17)	Normal (n = 26)	Gain (n = 140)
5 years (n = 183)	Body height, cm	111.8 ± 0.3^a	109.9 ± 0.5^b	108.9 ± 0.30^b	.0020
	Body length, cm	128.3 ± 2.1^a	122.7 ± 0.6^b	122.9 ± 0.30^b	.0001
	Chest girth, cm	161.4 ± 2.9^a	152.9 ± 1.3^b	154.6 ± 0.60^b	.0010
	Cannon width, cm	16.4 ± 0.4	16.2 ± 0.1	16.2 ± 0.07	.7300
	Body weight, kg	236.2 ± 13.6^a	204.3 ± 3.8^b	208.2 ± 2.10^b	.0010

Note: Different letters (^a,b) indicate significant differences (p < .05).CNV: copy number variation; SE: standard error.

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)).

Figure 2. Expression profiling of the GPC1 gene at two ages of Datong yak. (A) Expression pattern of the GPC1 gene in foetal yak. (B) Expression Pattern of the GPC1 gene in the adult yak. (C) Comparison of expression of GPC1 gen in two ages. The error bar indicated the standard deviation (SD). Different letters indicate significant difference (p < .05). CNV: copy number variation.

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) associated with higher mRNA expression levels than that of gain types, while no difference was observed between loss and normal copy number types.

Figure 3. Correlation analysis between the log₂ ratio of gene expression and CNVs of GPC1 in skeletal muscle of adult Datong yak. (A) Copy number variation of GPC1 gene (n = 41, M1-M41). (B) GPC1 gene expression level (n = 41, M1–M41). (C) Correlation between DNA copy number and mRNA GPC1 gene expression level. (4-6D) 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.

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₂ 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₂ 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 (Xu et al. 2016; 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 different scholars, and association between CNVs of the LEPR (Shi et al. 2016), MYH3 (Xu, Shi, et al. 2014), KCNJ12 (Zhou et al. 2016), MAPK10 (Liu et al. 2016), CYP4A11 (Yang et al. 2017), and MICAL-L2 (Xu et al. 2013) 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 (Liu et al. 2006) 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 (Velleman et al. 2013). 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 (Xu et al. 2013), 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 (Liu et al. 2016) 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.

Ethics of experimentation

The study was conducted between September 2016 and December 2017 at the Key Laboratory of Yak Breeding Engineering of Gansu Province, Lanzhou Institute of Husbandry and Pharmaceutical Sciences, and the data were obtained from yak breeding cooperatives of Gansu and Qinghai Province, China. All blood sample collection and body measurements were performed in strict accordance with the guide for the Care and Use of Laboratory Animals, Lanzhou Institute of Husbandry Animal and Pharmaceutical Sciences, China. Additionally, all the animals were slaughtered under anaesthesia, and all necessary efforts were made to minimise the risk of suffering. Thus, we were granted permission to perform this research study on yak, and the legal certificate number was SCXK (Gan) 2014-0002.

Acknowledgments

We thank Mr. Yong Feng Zhang for the guidance and technical assistance provided during the laboratory work and Dr. Hong Bo Wang, Ms. Li Xiaoxiao and Mr. Peng Tang for performing the blood sample collection.

Disclosure statement

The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

Additional information

Funding

This research was funded by the Innovation Project of Chinese Academy of Agricultural Sciences [grant number CAAS-ASTIP-2014-LIHPS-01] and the Program of National Beef Cattle and Yak Industrial Technology System [grant number CARS-37]. The APC was funded by the Lanzhou Institute of Animal Husbandry and Pharmaceutical Science.