Genome-wide association studies and genomic prediction of breeding values for calving performance and body conformation traits in Holstein cattle

Background Our aim was to identify genomic regions via genome-wide association studies (GWAS) to improve the predictability of genetic merit in Holsteins for 10 calving and 28 body conformation traits. Animals were genotyped using the Illumina Bovine 50 K BeadChip and imputed to the Illumina BovineHD BeadChip (HD). GWAS were performed on 601,717 real and imputed single nucleotide polymorphism (SNP) genotypes using a single-SNP mixed linear model on 4841 Holstein bulls with breeding value predictions and followed by gene identification and in silico functional analyses. The association results were further validated using five scenarios with different numbers of SNPs. Results Seven hundred and eighty-two SNPs were significantly associated with calving performance at a genome-wise false discovery rate (FDR) of 5%. Most of these significant SNPs were on chromosomes 18 (71.9%), 17 (7.4%), 5 (6.8%) and 7 (2.4%) and mapped to 675 genes, among which 142 included at least one significant SNP and 532 were nearby one (100 kbp). For body conformation traits, 607 SNPs were significant at a genome-wise FDR of 5% and most of them were located on chromosomes 5 (30%), 18 (27%), 20 (13%), 6 (6%), 7 (5%), 14 (5%) and 13 (3%). SNP enrichment functional analyses for calving traits at a FDR of 1% suggested potential biological processes including musculoskeletal movement, meiotic cell cycle, oocyte maturation and skeletal muscle contraction. Furthermore, pathway analyses suggested potential pathways associated with calving performance traits including tight junction, oxytocin signaling, and MAPK signaling (P < 0.10). The prediction ability of the 1206 significant SNPs was between 78 and 83% of the prediction ability of the BovineSNP50 SNPs for calving performance traits and between 35 and 79% for body conformation traits. Conclusions Various SNPs that are significantly associated with calving performance are located within or nearby genes with potential roles in tight junction, oxytocin signaling, and MAPK signaling. Combining the significant SNPs or SNPs within or nearby gene(s) from the HD panel with the BovineSNP50 panel yielded a marginal increase in the accuracy of prediction of genomic estimated breeding values for all traits compared to the use of the BovineSNP50 panel alone. Electronic supplementary material The online version of this article (10.1186/s12711-017-0356-8) contains supplementary material, which is available to authorized users.


Background
The profitability of dairy production depends on the ability of cows for high milk production, good health and fertility. Improving traits such as calving performance and body conformation reduces the culling rate, which in turn affects the profitability of the dairy cattle industry [1][2][3][4]. More than 50% of the first lactation heifers may require assistance at calving [5], which increases their culling risk by 18% and reduces their reproductive life [6]. In addition, dams with very difficult calving (score = 4) require eight more days to first service and 28 more days to conceive than those not needing assistance (score = 1) [2]. Furthermore, Dematawewa and Berger [7] reported that cows that experience extreme calving difficulty (score = 5) produce 703 kg of milk, 24 kg of fat and 21 kg of protein less per 305-d of lactation than unassisted cows (score = 1) over multiple parities. Body conformation traits are correlated to economically important traits in dairy cattle such as calving ease [8], longevity [9], lameness [10] and lifetime production efficiency [11]. Improving the accuracy of selection for calving performance and conformation traits would benefit the dairy industry as a whole and have a major impact on the profitability of individual farms and, thus, these traits are included in traditional and genomic breeding programs worldwide [12][13][14][15].
Since 2009, significant genetic improvement has been achieved for various dairy cattle traits through the use of high throughput single nucleotide polymorphism (SNP) panels for implementing genomic selection [16]. Increases in the rate of genetic progress were predicted to be as much as 100% [17], and data collected by the Canadian Dairy Network over the last few years support these predictions [18]. The Illumina Bovine SNP50 SNP chip (50 K; Illumina Inc., San Diego, USA) [19] has been used to identify significant regions that are associated with calving [20] and conformation traits [21]. Over the last four years, higher-density SNP chip panels have been developed, including the Illumina BovineHD BeadChip (HD; Illumina Inc., San Diego, USA) and, also, many animals have been sequenced [22]. These HD panels have not had a significant impact on genetic improvement programs, with very minor increases in prediction accuracies over the 50 K panel, either within or across breeds [23,24]. Nonetheless, HD panels can be used to improve the detection of regions that harbor causal mutations and, in which, genes of interest and new SNPs can be identified.
The objectives of our study were to identify genomic regions associated with 10 calving performance and 28 body conformation traits in Holstein cattle, to retrieve information from associated regions to increase the understanding of the underlying biology of these traits, and to test the predictive ability of these regions in young bulls.

Description of traits
The evaluated traits included 10 calving performance and 28 body conformation traits (Table 1). Calving performance traits were: maternal calving ease at first (heifer; CEh) and later calvings (cow; CEc), maternal calf survival at first (heifer; CSh) and later calvings (cow; CSc), direct calving ease at first (heifer; SCEh) and later calvings (cow; SCEc), direct calf survival at first (heifer; SCSh) and later calvings (cow; SCSc), calving ability index (CA), and daughter calving ability index (DCA). The phenotypes for calving ease were scored as 1 for unassisted, 2 for easy pull, 3 for hard pull and 4 for surgery. Calf survival was scored within 24 h after calving as 0 for dead and 1 for alive. The body conformation traits consisted of 23 linear descriptive traits scored on a 1 to 9 scale and four composite traits which were estimated from linear traits and overall score (www.holstein.ca). A detailed description of the traits and statistical methods for predicting the breeding values of the bulls was previously reported [9,13,[25][26][27].

Genotypes and imputation
Genotypes for HD (774,605 SNPs) and 50 K (40,666 SNPs) SNPs were provided by the Canadian Dairy Network (CND, Guelph, ON, Canada) for 2387 and 11,926 bulls, respectively, both registered and used in North America. The genotyping data was coded as 0, 1, and 2 for BB, AB and AA genotypes, respectively. For purposes of this study, reference population will refer to the animals genotyped with the HD panel used to create the library of haplotypes for the imputation from 50 K to HD genotypes, while target population refers to the bulls genotyped with the 50 K panel and imputed to HD genotypes.
A genotyping quality control (QC) on the HD genotypes removed 39,366 SNPs located on the sex chromosomes, 5427 individual SNPs with a call-rate lower than 90% and 31 SNPs with a heterozygosity that deviates by more than 15% from the expected value under Hardy-Weinberg equilibrium or that are not in Hardy-Weinberg equilibrium (P < 0.0000001). In addition, 3679 SNPs with a high Mendelian error rate (> 0.05) and 6902 SNPs that were assumed to be misplaced and identified based on linkage disequilibrium (LD) decay (see Additional file 1: Figure S1) were excluded. A total of 719,200 SNPs from the HD panel distributed over the 29 Bos taurus autosomes (BTA) passed all QC criteria and were used for imputation of the target population. These same QC criteria were applied to the 50 K genotypes and 39,148 SNPs remained for further analyses.
Genotype imputation from the 50 K to HD panel was performed by using the FImpute software and applying a population imputation approach [28]. SNPs with more that 10% of genotypes imputed by random filling based on the allele frequencies in the reference population had these imputed genotypes set to missing using the gp_thresh option in the FImpute software and, in a further QC step, these SNPs were excluded from further analyses. To estimate the accuracy of imputation, 387 animals were randomly selected from the reference population of animals with HD genotypes and their genotypes were masked down to the 50 K panel (39,148 SNPs). Then, the genotypes for these animals were imputed to HD by using the remaining 2000 animals in the reference population. The squared correlation (r 2 ) between imputed genotypes and true genotypes was used as a measure of imputation accuracy [29]. For the genome-wide association studies (GWAS) and genomic prediction of breeding values, imputation was carried out using all animals in the reference population (n = 2387).
After imputation, an additional QC was performed on the real and imputed genotypes. This step excluded 117,482 SNPs that had a minor allelic frequency (MAF) lower than 0.01 and one SNP for which heterozygosity deviated by more than 15% from the expected value under Hardy-Weinberg equilibrium. There were no animals and no individual SNPs with a call rate lower than 90%. A total of 601,717 SNPs distributed over the 29 bovine autosomes passed the QC procedure (see Additional file 2: Figure S2). To examine population stratification among the bulls in the training population, the HD genotypes were used to perform a multidimensional scaling analysis based on identical-by-state (IBS) pairwise distances, as implemented in PLINK software [30].

De-regressed estimated breeding values for the training population
Domestic official evaluations of 10 calving performance and 28 body conformation traits, and their associated reliabilities from the April 2009 genetic evaluation, were provided by the Canadian Dairy Network (CDN; www. cdn.ca) for 4841 progeny-tested Holstein bulls born between 1956 and 2007. These progeny-tested Holstein bulls had either a real HD (n = 287) or an imputed HD (n = 4554) genotype. The de-regressed estimated breeding values of bulls were computed by CDN following VanRaden et al. [31]: where De − EBV is the bull's de-regressed estimated breeding value for the trait of interest, PA is the parent average and Rel De−EBV is the reliability of the bull's EBV adjusted to remove the contribution of PA. For purposes of this study, training population refers to the population of bulls that had both pseudo-phenotypes (De − EBV ) and genotypes used for the estimation of the effects of SNPs, while validation population refers to the animals used for the validation of the effects of SNPs estimated by using the training population.
Pedigree data for 23,287 Holsteins individuals were obtained from CDN. The level of pedigree completeness of the training population was assessed using the pedigree completeness index (PCI) proposed by MacCluer et al. [32] and implemented in the CFC package [33]. The PCI was calculated as: where the K dam and K sire are the indexes for dams and the sires, respectively, given by: where a i is the proportion of known ancestors in generation i for either the dams or the sires, and g is the number of generations back (5 in our study).

Genome-wide association analyses
The GWAS were performed using an univariate single-SNP mixed linear model implemented in the snp1101 software [34]. The mixed model equations are described as: where Y is the vector of the De − EBV for each trait; 1 is a vector of ones, X is the vector of the animals' genotypes for a given SNP, coded as 0, 1 and 2; Z is the design matrix that assigns animals to De − EBV; μ is the overall mean, β is the linear regression coefficient (allele substitution effect) of the examined SNP; û is a vector of direct genomic values (DGV). Assumptions for the model were: u ∼ N 0, Gσ 2 g where G is the genomic relationship matrix, R is a diagonal matrix containing weights for the residual variance based on the reliabilities of the De − EBV. The diagonal elements of R are given by 1 1− 1 Rel i [35], where Rel i is the reliability of the ith De − EBV; σ 2 g is the additive genetic variance; and, σ 2 e is the residual variance. The values of σ 2 g and σ 2 e were estimated using the AI-REML algorithm implemented in the snp1101 software [34]. Allele frequencies used to calculate G were estimated from the observed genotype data. The G matrix was calculated as: where C is the genotypic coefficient matrix with dimensions equal to the number of genotyped animals by the number of SNPs and P i is the allele frequency of the ith SNP. The relationship between individuals j and k is G jk divided by the square root of the diagonal values of j (G jj ) and k (G kk ) individuals [35]. The significance of the effects of SNPs was determined by using a genome-wise false discovery rate (FDR) of 5% [36]. The more lenient FDR threshold of 5% was used to increase the power to detect SNPs with small effects since traits of interest may be controlled by many QTL with a small effect [37]. The genomic inflation factor (λ) and quantile-quantile (Q-Q) plots were applied to assess the inflation of the test statistics due to population stratification by comparing the genome-wide distribution of the test statistic (− log10 of P-values) with the expected Chi squared distribution. The significant SNPs at a genome-wise FDR of 5% were aligned to the publicly available QTL in the Animal QTLdb (Release 26, access date: April 27th, 2015) [38].

In-silico functional analyses
Candidate SNP enrichment analyses were performed on the GWAS results using the SNP2GO R package [39] for gene ontology (GO). To infer overrepresentation of candidate SNPs (i.e., the list of significant SNPs at a P value < 0.05), the Bos taurus annotations from Ensembl version 75 for the UMD 3.1 assembly, the associated GO terms, and the list of candidate and non-candidate (i.e., non-significant) SNPs were used as input for SNP2GO. In the SNP2GO analyses, significant SNPs were assigned to plus or minus 50 nucleotides up-and down-stream of the corresponding genes. To account for multiple testing in the SNP enrichment analysis, a FDR of 1% was applied.
Pathway analyses were performed on the list of genes obtained by mapping the significant SNPs from the GWAS to their corresponding genes or nearby genes at a distance of 10 kbp using NGS-SNP [40], the bovine genome assembly UMD3.1 and the Ensembl database 75_31. The 10 kbp distance was used to exploit the expected LD (r 2 ) between pairs of syntenic SNPs, which had an average r 2 of 0.58 [41]. The pathway analyses were performed using the web-based Database for Annotation, Visualization, and Integrated Discovery (DAVID) version 6.8 [42].

Validation study
We performed a forward validation study, in which the Canadian domestic estimated breeding values from the April 2009 genetic evaluation were used to estimate SNP associations and direct genomic breeding values (DGV).

Unproven young bulls in the April 2009 genetic evaluation with estimated breeding values in December 2014
were considered for the validation population ( Table 1). The DGV for the validation group were estimated by using the genomic best linear unbiased prediction method (GBLUP) [35,43,44] implemented in the gebv software [45]. In this method, the DGV were estimated using the following mixed model equations: where Y is the vector of De − EBV from the April 2009 genetic evaluation for the trait of interest for genotyped bulls, μ is the overall mean, 1 is a vector of ones, Z is the design matrix that assigns animals to De − EBV , û is a vector of DGV. Assumptions for this model are: g is the additive genetic variance; G is the genomic relationship matrix [35], A is the numerator relationship matrix and w is the weight put on G (0.8). As in Model (3), R is a diagonal matrix with elements defined as R ii = Different scenarios for the validation analyses were tested depending on the number of SNPs used in the prediction. Scenario 1 used all the SNPs in the 50 K panel that passed the QC (38,724 SNPs), Scenario 2 used all the SNPs in the HD panel that passed the QC (607,794 SNPs), Scenario 3 used only a set of significant (genomewise FDR of 5%) SNPs for all calving performance and conformation traits, Scenario 4 combined SNPs from the 50 K panel and the significant SNPs not in the 50 K panel, and Scenario 5 combined a set of significant SNPs for calving performance and conformation traits within or nearby the gene (± 100 kbp around the gene) plus the 50 K panel. The accuracy of genomic breeding values was calculated as the squared correlation (r 2 ) between DGV and official genetic evaluations from December 2014 for the validation group.

Population structure
The results from the multidimensional scaling analysis showed no divergent clusters within this population (see Additional file 3: Figure S3). In addition, the pedigree analysis indicated that the overall PCI considering five generations back was equal to 0.62. Although the PCI for the reference population was high, modeling the relationship (5) between animals using the pedigree information may not completely account for population structure because of missing pedigree records. The current Holstein cattle population was characterized by a high rate of inbreeding, genetic diversity loss and small effective population size [46,47] and would probably have cryptic relationships that are not described by the available pedigree. Thus, the population structure and cryptic relatedness not accounted for by an incomplete pedigree should be accounted for in GWAS to avoid spurious associations [48].

Accuracy of genome-wide imputation
After applying QC criteria, 719,200 SNPs distributed over the 29 bovine autosomes in the reference population (2387 individuals) were used for genome-wide imputation of the target population (11,926 animals) for which 39,148 SNPs were available. The squared correlation (r 2 ) between imputed genotypes and true genotypes of the target population was estimated at 99.29%, which indicated that the imputation was highly accurate. The allelic r 2 was used as a measure of imputation accuracy because it adequately reflects the imputation accuracy of SNPs with a low MAF [29]. In this study, the accuracy of Genomic inflation factor is around 1 indicating that there is no population stratification imputation was high for several reasons. First, the reference population with HD genotypes was relatively large (n = 2387) compared with other studies that reported high imputation accuracies by using the FImpute software and with less individuals in the reference group (e.g., 1733) than those used in our study [28]. In addition, the accuracy of imputation from 50 K to HD using FImpute was also high (99.96%) in a Holstein reference population of limited size (n = 1406) [24]. Second, the reference population in our study shares relatively longrange haplotypes [28] with high phasing accuracy which also contributes to higher imputation accuracies [49]. Third, this reference population is also closely related to the target population which helped increase the accuracy of imputation, even for SNPs with a low MAF (≤ 0.05) [28]. Furthermore, an internal score based on the accuracy of imputation for each SNP allele was used to set less accurate imputed alleles to missing and the SNPs with a call rate lower than 90% were excluded in the QC before GWAS. In addition, SNPs with a MAF lower than 1% were initially excluded, which resulted in more accurate genotypes for the GWAS. The use of imputation increased both the numbers of SNPs and animals included in the association analysis, which increased the power of QTL detection [50] and the accuracy of genomic prediction of breeding values [51,52]. Combining SNP panels with different densities and using imputed genotypes in QTL mapping is a common approach in human GWAS [53] and in dairy cattle genomic evaluations [24,54].

Genome-wide association for calving performance traits
In the GWAS, family and cryptic stratifications were accounted for by incorporating the full genomic covariance among individuals. The genomic inflation factor was used to assess bias in the test statistics. The average genomic inflation factor was equal to 0.98 and ranged Genomic inflation factor is around 1 indicating that there is no population stratification from 0.93 for SCEc to 1.00 for heel depth (see Additional file 4: Figure S4), which suggests that any potential bias due to population stratification was addressed [55][56][57]. Furthermore, accounting for the reliability of the de-regressed EBV allowed the use of more records in the estimation step by including de-regressed EBV with lower reliability, which subsequently improved the accuracy of GWAS [50,52]. For calving performance, 782 SNPs significant at a genome-wide FDR of 5% were identified. Most of these significant SNPs were on BTA18 (71.87%), 17 (7.42%), 5 (6.78%), 7 (2.43%), 19 (1.92%), 21 (1.66%), 29 (1.66%), 1 (1.02%) and 10 (1.02%). The significant SNPs were mapped to 675 genes, among which 142 included at least one significant SNP and 532 were nearby (100 kbp) a SNP. A strong peak on BTA18 that affected maternal and direct calving ease and calf survival in the first and later calvings was identified (Figs. 1, 2, 3, 4, 5). Previously, Sahana et al. [58] and Hoglund et al. [59] reported QTL on BTA18 that were associated with direct and maternal calving ease, calf size, stillbirth and birth index in Danish and Swedish Holstein cattle. In addition, Saatchi et al. [60] reported a QTL on BTA18 (at 54.37 Mb) that was associated with maternal calving ease in Simmental beef cattle. Also, the significant associations identified for most calving performance traits on BTA5, 6,7,9,13,17,18,19,20 and 23 are in agreement with previous studies [58][59][60]. The significant associations were located within the confidence interval of previously detected QTL (see Additional file 5: Table S1). Thirty significant (FDR of 5%) SNPs were common between maternal and direct calving performance traits (Fig. 6).

Fig. 6
Number of significant (at the genome-wise false discovery rate of 5%) SNPs that have a pleiotropy effect on calving performance using the Illumina BovineHD BeadChip in Holstein cattle. Calving performance traits are calving ability (CA); daughter calving ability (DCA); maternal calving ease at first calving (heifer; CEh); maternal calf survival at first calving (heifer; CSh); maternal calving ease at later calvings (cow; CEc); maternal calf survival at later calvings (cow; CSc); direct calving ease at first calving (heifer; SCEh); direct calf survival (heifer; SCSh); direct calving ease at later calvings (cow; SCEc); sire calf survival (cow; SCSc) Several genes that contain significant SNPs for calving performance traits were identified ( Table 2) and are involved in biological processes such as lipid metabolism, immunity, reproduction and anatomical structure development. For example, kallikrein-related peptidase 14 (KLK14) is involved in mating and reproduction processes in a multicellular organism [61,62]. KLK14 was reported to play an important role in the synergistic effects between estrogens and androgens [61] and at the onset of parturition [62]. Another important gene was killer cell immunoglobulin like receptor, two Ig domains and long cytoplasmic tail 5A (KIR2DL5A), which is attributed to Graft-versus-host disease and natural killer cell mediated cytotoxicity [63]. The membrane bound O-acyltransferase domain containing 7 (MBOAT7) gene is involved in the metabolism of lipids and lipoproteins, including the pathway of glycerophospholipids metabolism, which is linked to energy balance  metabolites including non-esterified fatty acids (NEFA), beta-hydroxybutyrate (BHBA) and glucose in animals [64]. It is well documented that the period around calving is associated with several biological processes related to fat, protein, glucose and mineral metabolism as well as complex hormonal changes [65]. Monitoring and nutritional modeling of the animal's energy balance to prevent negative energy balance during the transition period of dairy cows is an important management practice that can reduce the incidence of metabolic and infectious diseases [66]. Energy balance could be associated with the length of the dry period and subsequently with calving difficulty which is generally common under a short (< 60 days) compared to a long dry period [67].
In this study, we confirmed the relationship between the protective functions of the immune system and calving ease and survival by the detection of several genes such as Fc fragment of IgA receptor (FCAR). FCAR includes a significant SNP (rs134066287) which is associated with calving ease and survival at first and later calving (Table 2). FCAR is also involved in several immune defense processes such as phagocytosis and Staphylococcus aureus infection pathways [68]. The expression profile of fetal membranes, myometrium and cervix tissues showed that inflammatory response was associated with labor [69].

Genome-wide association for body conformation traits
Six hundred and seven SNPs were statistically significant at a genome-wise FDR of 5% for body conformation traits (Table 3). These significant SNPs were mapped to 553 genes, among which 89 genes included at least one significant SNP and 464 were nearby (100 kbp) one. Peaks indicating association were found on BTA5, 6,7,13,14,18,and 20 (Figs. 7,8,9,10,11,12,13). Four significant SNPs (FDR of 5%) overlapped between bone quality and pin width (Fig. 14). For example, the arresting domain-containing 3 (ARRDC3) gene on BTA7 was associated with body conformation (e.g., bone quality and chest width) and calving performance traits. The ARRDC3 gene was previously reported as a candidate gene with a pleiotropic effect on birth and weaning weights, direct and maternal calving ease and carcass traits in beef cattle [60]. Furthermore, this study identified other significant polymorphisms for body depth within the carnitine palmitoyl transferase 1C (CPT1C) gene, which is involved in PPAR signaling, adipocytokine signaling and fatty acid metabolism pathways linked to the animal's energy balance [70]. Interestingly, we identified a SNP within the diacylglycerol O-acyltransferase 1 (DGAT1) gene that was significantly associated with bone quality. Bone quality is an important trait that affects functional longevity of  Holstein cows [9]. Kaupe et al. [71] identified a SNP within DGAT1 that is associated with rump width and strength, which suggests a possible association between DGAT1 and functional longevity.

Pleiotropic polymorphisms for calving performance and body conformation traits
One hundred and eighty-three SNPs were associated with calving and conformation traits. Sixteen SNPs on BTA18 in the region between 56.9 and 59.9 Mbp were associated with calving performance and rump traits particularly pin width (Fig. 15). This region harbours the SH3 and multiple ankyrin repeat domains 1 (SHANK1) gene, the myosin heavy chain 14 (MYH14) gene, and the cytosolic thiouridylase subunit 1 (CTU1) gene, which are involved in tight junction, viral myocarditis, regulation of actin cytoskeleton, glutamatergic synapse, and sulfur relay system pathways [72]. The MYH14 gene is associated with calving and body conformation and with one of the most significant SNPs (rs137554975) identified for eight calving traits. MYH14 is also involved in the regulation of actin cytoskeleton and tight junction pathways which are critical for myometrial functions during parturition [73]. The pleiotropic associations between calving performance and rump traits support the known genetic correlation between characteristics of the pelvis area (e.g., rump traits) and calving ease and calf survival [8]. A cow with a wide pin, long sloping rump, and slight slope from pin bone to thurl is known to be able to calve easily [8,74]. An unfavorable genetic association between height of pin bones and slope to the cow's birth canal was observed, which suggests that higher pin bones may result in a tight birth canal causing calving difficulty [75].

Functional consequences of significant SNPs and in silico functional analyses
Most of the identified SNPs for calving and/or body conformation traits were intergenic (44.4%), intronic (31%), downstream (11.4%), or upstream (9.8%) of a gene. As shown in Table 4, we identified several SNPs with high functional interest (e.g., non-synonymous coding and splice site intronic variants). For instance, a stop-lost mutation was identified within the HORMA domain containing 2 (HORMAD2) gene and was significantly associated with CA, SCEc and SCEh. HORMAD2 is involved in the M-phase of the mitotic cell cycle and, in humans, a polymorphism within this gene is known to cause meiotic arrest leading to human azoospermia [76]. The SNP enrichment analyses provided potential GO terms for calving and conformation traits including biological mechanisms, molecular function and cellular component (see Additional file 6: Table S2). Biological processes including musculoskeletal movement (GO:0050881), meiotic cell cycle (GO:0051321), oocyte maturation (GO:0001556), protein localization to synapse (GO:0035418) and skeletal muscle contraction (GO:0003009) were enriched (FDR of 1%) for calving ease and survival traits. Also, biological processes such as steroid dehydrogenase activity (GO:0016229) and 3-beta-hydroxy-delta5-steroid dehydrogenase activity (GO:0016229) were over-represented (FDR of 1%) for calving ease in heifers. The steroids are very important for the control of the synthesis of oxytocin receptors and can facilitate the initiation of parturition [62].
Thirty-three genes were enriched (P < 0.10) in 13 pathways, in which tight junction, oxytocin signaling and MAPK signaling are associated with calving performance traits ( Table 5). The tight junction proteins are highly expressed in the final stages of cervical ripening and dilation in preparation for parturition [76,77]. Furthermore, mitogen-activated protein kinase (MAPK) signaling pathway was enriched (P = 0.065) based on six genes. MAPK signaling cascades are critical in the regulation of several mechanisms including uterine contractility and myometrial cell proliferation [73]. In agreement with our findings, previous GWAS in Angus and Hereford cattle Manhattan plot in which the genomic coordinates of SNPs are displayed along the horizontal axis, the negative logarithm of the association P-value for each SNP is displayed on the vertical axis, and the dark red line is the significance threshold line at the genome-wise false discovery rate of 5%. b Quantile-quantile (QQ) plot showing the late separation between observed and expected values. Genomic inflation factor is around 1 indicating that there is no population stratification detected the MAPK signaling pathway as having a pleiotropic effect on birth weight, calving ease (direct and maternal) [60]. In addition, oxytocin signaling pathway which comprises the main drivers of calving [77] was also enriched (P < 0.05) for five genes (see Additional file 7: Figure S5).

Validation of genomic predictions
In this study, there was no substantial difference between prediction accuracies of the DGV from the 50 K or HD panels for young bulls (Fig. 16). This finding is in agreement with the study of Erbe et al. [23], in which the use of 624,213 SNPs provided no extra gain in accuracy than the 50 K panel when using GBLUP in Jersey and Holstein cattle. The reason why the use of both panels yielded similar results may be that the imputed HD panel captures the same loci with a large effect than those captured by the 50 K panel and that other QTL with a small effect captured by the imputed HD panel explain just a small proportion of the genetic variance and, consequently, result in only a small increase in accuracy of prediction of DBV [24]. A slight improvement in the accuracy of DGV when combining the 50 K panel with the significant SNPs or the SNPs located nearby genes compared to using only the SNPs from the 50 K panel may indicate that the HD SNP panel contains SNPs that are potentially linked to causal mutations. However, a decrease of about 6.5% in prediction accuracy on average for calving traits was observed when only significant SNPs (1206) were used compared to the 50 K panel. Nonetheless, the small number of significant SNPs (1206) with a prediction accuracy up to 0.34 supported the informativeness of the candidate genes and the biological mechanisms associated with the traits.
In general, the validation of significantly associated SNPs in a different population is a necessary step before application in breeding programs in order to estimate the accuracy of the DGV properly and to test for potentially spurious associations in the original findings. The calving ability index and direct calving ease and survival traits were not evaluated in the validation analyses because of the limited number of young bulls. The GBLUP method was used because it is the official method for genomic evaluations in Canadian dairy cattle [35], it was recommended by earlier studies [78,79], and has been shown to be as accurate as other statistical methods [31]. The in which the genomic coordinates of SNPs are displayed along the horizontal axis, the negative logarithm of the association P-value for each SNP is displayed on the vertical axis, and the dark red line is the significance threshold line at the genome-wise false discovery rate of 5%. b Quantilequantile (QQ) plot showing the late separation between observed and expected values. Genomic inflation factor is around 1 indicating that there is no population stratification current study also identified several SNPs with pleiotropic effects, which are associated with known biological pathways involved in calving performance and conformation traits. The correlations between the breeding values of calving performance ranged from 0.11 to 0.82 [80], which suggests that common genes control more than one calving trait, as observed in this study.

Conclusions
We identified various SNPs that are significantly associated with calving performance on chromosomes 5, 18, and 19 and are located within or nearby genes with potential roles in tight junction, oxytocin signaling, and MAPK signaling. Sixteen SNPs within or nearby the SHANK1, MYH14, and CTU1 genes on BTA18 showed     a pleiotropic effect on calving performance and body conformation traits. In total, eight, six and four QTL were confirmed for direct calving ease, maternal calving ease, and direct stillbirth, respectively. Combining significant SNPs from the GWAS with SNPs in the 50 K panel for genomic evaluation slightly increased the prediction accuracy of genomic breeding values for calving and body conformation traits. Validation of the SNPs in independent populations, followed by the possible identification of the causal mutations within the validated candidate genes are the logical next research steps.