Skip to main content
  • Research Article
  • Open access
  • Published:

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

Abstract

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.

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.

Methods

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

Table 1 Heritability and number of records used in the association and validation analyses

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 (r2) 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]:

$${\text{De}} - {\text{EBV}} = {\text{PA}} + \left( {{\text{EBV}} - {\text{PA}}} \right)/{\text{Rel}}_{{{\text{De}} - {\text{EBV}}}} ,$$
(1)

where \({\text{De}} - {\text{EBV}}\) is the bull’s de-regressed estimated breeding value for the trait of interest, \({\text{PA}}\) is the parent average and \({\text{Rel}}_{{{\text{De}} - {\text{EBV}}}}\) is the reliability of the bull’s EBV adjusted to remove the contribution of \({\text{PA}}\). For purposes of this study, training population refers to the population of bulls that had both pseudo-phenotypes (\({\text{De}} - {\text{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 (\({\text{PCI}}\)) proposed by MacCluer et al. [32] and implemented in the CFC package [33]. The \({\text{PCI}}\) was calculated as:

$${\text{PCI}} = \frac{{K_{dam} K_{sire} }}{{K_{dam} + K_{sire} }},$$
(2)

where the K dam and K sire are the indexes for dams and the sires, respectively, given by:

$$K = \frac{1}{g}\mathop \sum \limits_{i = 1}^{g} a_{i}$$

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:

$$\left[ {\begin{array}{*{20}c} {{\bf{1}} '{\mathbf{R}}^{ - 1} {\bf{1}}} & {{\bf{1}} '{\mathbf{R}}^{ - 1} {\mathbf{X}}} & {{\bf{1}} '{\mathbf{R}}^{ - 1} {\mathbf{Z}}} \\ {{\mathbf{X}} '{\mathbf{R}}^{ - 1} {\bf{1}}} & {{\mathbf{X}} '{\mathbf{R}}^{ - 1} {\mathbf{X}}} & {{\mathbf{X}} '{\mathbf{R}}^{ - 1} {\mathbf{Z}}} \\ {{\mathbf{Z}} '{\mathbf{R}}^{ - 1} {\bf{1}}} & {{\mathbf{Z}} '{\mathbf{R}}^{ - 1} {\mathbf{X}}} & {{\mathbf{Z^{\prime}R}}^{ - 1} {\mathbf{Z}} + {\mathbf{G}}^{ - 1} *\frac{{\upsigma_{\text{e}}^{2} }}{{\upsigma_{\text{g}}^{2} }}} \\ \end{array} } \right] \times \left[ {\begin{array}{*{20}c} {{\hat{\upmu }}} \\ {{\hat{\upbeta }}} \\ {{\hat{\mathbf{u}}}} \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {{\bf{1}} '{\mathbf{R}}^{ - 1} {\mathbf{Y}}} \\ {{\mathbf{X}} '{\mathbf{R}}^{ - 1} {\mathbf{Y}}} \\ {{\mathbf{Z}} '{\mathbf{R}}^{ - 1} {\mathbf{Y}}} \\ \end{array} } \right],$$
(3)

where \({\mathbf{Y}}\) is the vector of the \({\text{De}} - {\text{EBV}}\) for each trait; 1 is a vector of ones, \({\mathbf{X}}\) is the vector of the animals’ genotypes for a given SNP, coded as 0, 1 and 2; \({\mathbf{Z}}\) is the design matrix that assigns animals to \({\text{De}} - {\text{EBV}}\); \({\hat{\upmu }}\) is the overall mean, \({\hat{\upbeta }}\) is the linear regression coefficient (allele substitution effect) of the examined SNP; \({\hat{\mathbf{u}}}\) is a vector of direct genomic values (DGV). Assumptions for the model were: \({\mathbf{u}}\sim N\left( {0, {\mathbf{G}}\upsigma_{\text{g}}^{2} } \right)\) where \({\mathbf{G}}\) is the genomic relationship matrix, \({\mathbf{R}}\) is a diagonal matrix containing weights for the residual variance based on the reliabilities of the \({\text{De}} - {\text{EBV}}\). The diagonal elements of \({\mathbf{R}}\) are given by \(\left( {\frac{1}{{1 - \frac{1}{{{\text{Rel}}_{\text{i}} }}}}} \right)\) [35], where \({\text{Rel}}_{\text{i}}\) is the reliability of the \({\text{i}}\)th \({\text{De}} - {\text{EBV}}\); \(\upsigma_{\text{g}}^{2}\) is the additive genetic variance; and, \(\upsigma_{\text{e}}^{2 }\) is the residual variance. The values of \(\upsigma_{\text{g}}^{2}\) and \(\upsigma_{\text{e}}^{2}\) were estimated using the AI-REML algorithm implemented in the snp1101 software [34]. Allele frequencies used to calculate \({\mathbf{G}}\) were estimated from the observed genotype data. The \({\mathbf{G}}\) matrix was calculated as:

$${\mathbf{G}} = {\mathbf{CC^{\prime}}}/2\sum {\text{P}}_{\text{i}} (1 - {\text{P}}_{\text{i}} ) ,$$
(4)

where \({\mathbf{C}}\) is the genotypic coefficient matrix with dimensions equal to the number of genotyped animals by the number of SNPs and \({\text{P}}_{\text{i}}\) is the allele frequency of the \({\text{i}}\)th SNP. The relationship between individuals \({\text{j}}\) and \({\text{k}}\) is \({\text{G}}_{\text{jk}}\) divided by the square root of the diagonal values of \({\text{j}}\) (\({\text{G}}_{\text{jj}}\)) and \({\text{k}}\) (\({\text{G}}_{\text{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 (r2) between pairs of syntenic SNPs, which had an average r2 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:

$$\left[ {\begin{array}{*{20}c} {{\bf{1}}'{\mathbf{R}}^{ - 1} {\bf{1}}} \\ {{\mathbf{Z^{\prime}R}}^{ - 1} {\bf{1}}} \\ \end{array} \begin{array}{*{20}c} {{\bf{1}} '{\mathbf{R}}^{ - 1} {\mathbf{Z}}} \\ { {\mathbf{Z}}^{ '} {\mathbf{R}}^{ - 1} {\mathbf{Z}} + {\mathbf{G}}_{\text{w}}^{ - 1} \frac{{\upsigma_{\text{e}}^{2} }}{{\upsigma_{\text{g}}^{2} }}} \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {{\hat{\upmu }}} \\ {{\hat{\mathbf{u}}}} \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {{\bf{1}} '{\mathbf{R}}^{ - 1} {\mathbf{Y}}} \\ {{\mathbf{Z}} '{\mathbf{R}}^{ - 1} {\mathbf{Y}}} \\ \end{array} } \right],$$
(5)

where \({\mathbf{Y}}\) is the vector of \({\text{De}} - {\text{EBV}}\) from the April 2009 genetic evaluation for the trait of interest for genotyped bulls, \({\hat{\upmu }}\) is the overall mean, 1 is a vector of ones, \({\mathbf{Z}}\) is the design matrix that assigns animals to \({\text{De}} - {\text{EBV}}\), \({\hat{\mathbf{u}}}\) is a vector of DGV. Assumptions for this model are: \({\mathbf{u}}{\sim}N\left( {0, {\mathbf{G}}_{\text{w}} {{\upsigma }}_{\text{g}}^{2} } \right)\), where \({\mathbf{G}}_{\text{w}} = {\text{w}}^{ *} {\mathbf{G}} + \left( {1 - {\text{w}}} \right)^{{*}} {\mathbf{A}}\) and \(\upsigma_{\text{g}}^{2}\) is the additive genetic variance; \({\mathbf{G}}\) is the genomic relationship matrix [35], \({\mathbf{A}}\) is the numerator relationship matrix and \({\text{w}}\) is the weight put on \({\mathbf{G}}\) (0.8). As in Model (3), \({\mathbf{R}}\) is a diagonal matrix with elements defined as Rii = \(\left( {\frac{1}{{1 - \frac{1}{{{\text{Rel}}_{\text{i}} }}}}} \right)\) where \({\text{Rel}}_{\text{i}}\) is the reliability of the \({\text{i}}\)th \({\text{De}} - {\text{EBV}}\); \(\upsigma_{\text{e}}^{2}\) is the error variance.

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 (genome-wise 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 (r2) between DGV and official genetic evaluations from December 2014 for the validation group.

Results and discussion

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 \({\text{PCI}}\) considering five generations back was equal to 0.62. Although the \({\text{PCI}}\) for the reference population was high, modeling the relationship 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 (r2) 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 r2 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 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 long-range 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 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. 1
figure 1

Manhattan (a) and Q–Q (b) plots from GWAS for calving ability index and daughter calving ability index using Illumina BovineHD BeadChip in Holstein cattle. a 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 (Q–Q) plot showing the late separation between observed and expected values. Genomic inflation factor is around 1 indicating that there is no population stratification

Fig. 2
figure 2

Manhattan (a) and Q–Q (b) plots from GWAS for maternal calving ease at first calving and maternal calving ease at later calvings using Illumina BovineHD BeadChip in Holstein cattle. a 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 (Q–Q) plot showing the late separation between observed and expected values. Genomic inflation factor is around 1 indicating that there is no population stratification

Fig. 3
figure 3

Manhattan (a) and Q–Q (b) plots from GWAS for maternal calf survival at first calving and maternal calf survival at later calvings using Illumina BovineHD BeadChip in Holstein cattle. a 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 (Q–Q) plot showing the late separation between observed and expected values. Genomic inflation factor is around 1 indicating that there is no population stratification

Fig. 4
figure 4

Manhattan (a) and Q–Q (b) plots from GWAS for direct calving ease at first calving and direct calving ease at later calvings using Illumina BovineHD BeadChip in Holstein cattle. a 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 (Q–Q) plot showing the late separation between observed and expected values. Genomic inflation factor is around 1 indicating that there is no population stratification

Fig. 5
figure 5

Manhattan (a) and Q–Q (b) plots from GWAS for direct calf survival at first calving and direct calf survival at later calvings using Illumina BovineHD BeadChip in Holstein cattle. a 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 (Q–Q) plot showing the late separation between observed and expected values. Genomic inflation factor is around 1 indicating that there is no population stratification

Fig. 6
figure 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].

Table 2 The 10 most significant SNPs associated with each calving performance trait at the genome-wise false discovery rate of 5%

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.

Table 3 The most highly significant SNPs associated with each body conformation trait at the genome-wise false discovery rate of 5%
Fig. 7
figure 7

Manhattan (a) and Q–Q (b) plots from GWAS for dairy strength and pin setting using the Illumina BovineHD BeadChip in Holstein cattle. a 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

Fig. 8
figure 8

Manhattan (a) and Q–Q (b) plots from GWAS for pin width and udder depth using the Illumina BovineHD BeadChip in Holstein cattle. a 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

Fig. 9
figure 9

Manhattan (a) and Q–Q (b) plots from GWAS for fore udder attachment and front teat placement using the Illumina BovineHD BeadChip in Holstein cattle. a 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

Fig. 10
figure 10

Manhattan (a) and Q–Q (b) plots from GWAS for rear attachment height and teat length using the Illumina BovineHD BeadChip in Holstein cattle. a 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

Fig. 11
figure 11

Manhattan (a) and Q–Q (b) plots from GWAS for bone quality and stature using the Illumina BovineHD BeadChip in Holstein cattle. a 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

Fig. 12
figure 12

Manhattan (a) and Q–Q (b) plots from GWAS for chest width and body depth using the Illumina BovineHD BeadChip in Holstein cattle. a 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

Fig. 13
figure 13

Manhattan (a) and Q–Q (b) plots from GWAS for angularity using the Illumina BovineHD BeadChip in Holstein cattle. a 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

Fig. 14
figure 14

Number of significant (at the genome wise false discovery rate of 5%) SNPs that have a pleiotropic effect on body conformation traits using the Illumina BovineHD BeadChip in Holstein cattle

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

Fig. 15
figure 15

Number of significant (at the genome-wise false discovery rate of 5%) SNPs that have a pleiotropic effect on calving performance and body conformation traits using the Illumina BovineHD BeadChip in Holstein cattle. Rump is the rump traits including pin setting and pin width

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

Table 4 Significant SNPs and their genes associated with calving and conformation traits in Holstein cattle and causing potential functional consequences

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

Table 5 KEGG biological pathways enriched with the identified candidate genes for calving performance and body conformation traits

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.

Fig. 16
figure 16

Squared correlation (r2) between direct genomic breeding values (DGV) and the corresponding de-regressed EBV for validation bulls using different SNP sets. a Calving performance traits are daughter calving ability (DCA); maternal calving ease at first calving (heifer; CEh); maternal calf survival at first calving (heifer; CSh) and maternal calf survival at later calvings (cow; CSc). b Body conformation traits are conformation score (CONF); rump (RU); mammary system (MS); feet & legs (FL); dairy strength (DS); rump angle (RAN); pin setting (PS); pin width (PW); loin strength (LS); udder depth (UD); udder texture (UT); median suspensory ligament (MSL); fore udder attachment (FA); front teat placement (FTP); teat length (TL); foot angle (FAN); heel depth (HD); bone quality (BQ); leg side view (LSV); set of rear legs (SRL); leg rear view (LRV); stature (ST); height at front end (FE); chest width (CW); body depth (BD) and angularity (ANG). Set 1 (38,405 SNPs) = set of SNPs from the 50K SNP panel; set 2 (601,717 SNPs) = set of SNPs from the HD SNP panel; set 3 (1206 SNPs) = set of significant (genome-wise false discovery rate of 5%) SNPs for calving and body conformation traits; set 4 (39,564 SNPs) = set of combined SNPs from sets 1 and 3; set 5 (39,408 SNPs) = set of combined SNPs from set 1 and a set of significant SNPs for all traits within or nearby genes (± 100 kbp around the gene); reliability = average reliabilities of the EBV from the December 2014 genetic evaluation for validation population

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

References

  1. Beaudeau F, Seegers H, Ducrocq V, Fourichon C, Bareille N. Effect of health disorders on culling in dairy cows: a review and a critical discussion. Ann Zootech. 2000;49:293–311.

    Article  Google Scholar 

  2. Eaglen SAE, Coffey MP, Woolliams JA, Mrode R, Wall E. Phenotypic effects of calving ease on the subsequent fertility and milk production of dam and calf in UK Holstein-Friesian heifers. J Dairy Sci. 2011;94:5413–23.

    Article  CAS  PubMed  Google Scholar 

  3. McGuirk BJ, Forsyth R, Dobson H. Economic cost of difficult calvings in the United Kingdom dairy herd. Vet Rec. 2007;161:685–7.

    Article  CAS  PubMed  Google Scholar 

  4. Eaglen SAE, Coffey MP, Woolliams JA, Wall E. Direct and maternal genetic relationships between calving ease, gestation length, milk production, fertility, type, and lifespan of Holstein-Friesian primiparous cows. J Dairy Sci. 2013;96:4015–25.

    Article  CAS  PubMed  Google Scholar 

  5. Lombard JE, Garry FB, Tomlinson SM, Garber LP. Impacts of dystocia on health and survival of dairy calves. J Dairy Sci. 2007;90:1751–60.

    Article  CAS  PubMed  Google Scholar 

  6. López de Maturana E, Ugarte E, González-Recio O. Impact of calving ease on functional longevity and herd amortization costs in Basque Holsteins using survival analysis. J Dairy Sci. 2007;90:4451–7.

    Article  PubMed  Google Scholar 

  7. Dematawewa CM, Berger PJ. Effect of dystocia on yield, fertility, and cow losses and an economic evaluation of dystocia scores for Holsteins. J Dairy Sci. 1997;80:754–61.

    Article  CAS  PubMed  Google Scholar 

  8. Ali TE, Burnside EB, Schaeffer LR. Relationship between external body measurements and calving difficulties in Canadian Holstein-Friesian cattle. J Dairy Sci. 1984;67:3034–44.

    Article  CAS  PubMed  Google Scholar 

  9. Sewalem A, Kistemaker GJ, Miglior F, Van Doormaal BJ. Analysis of the relationship between type traits and functional survival in Canadian Holsteins using a Weibull proportional hazards model. J Dairy Sci. 2004;87:3938–46.

    Article  CAS  PubMed  Google Scholar 

  10. Chapinal N, Koeck A, Sewalem A, Kelton DF, Mason S, Cramer G, et al. Genetic parameters for hoof lesions and their relationship with feet and leg traits in Canadian Holstein cows. J Dairy Sci. 2013;96:2596–604.

    Article  CAS  PubMed  Google Scholar 

  11. Sawa A, Bogucki M, Krężel-Czopek S, Neja W. Relationship between conformation traits and lifetime production efficiency of cows. ISRN Vet Sci. 2013;2013:124690.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. VanRaden PM, Sanders AH, Tooker ME, Miller RH, Norman HD, Kuhn MT, et al. Development of a national genetic evaluation for cow fertility. J Dairy Sci. 2004;87:2285–92.

    Article  CAS  PubMed  Google Scholar 

  13. Jamrozik J, Fatehi J, Kistemaker GJ, Schaeffer LR. Estimates of genetic parameters for Canadian Holstein female reproduction traits. J Dairy Sci. 2005;88:2199–208.

    Article  CAS  PubMed  Google Scholar 

  14. Weigel KA, Rekaya R. Genetic parameters for reproductive traits of Holstein cattle in California and Minnesota. J Dairy Sci. 2000;83:1072–80.

    Article  CAS  PubMed  Google Scholar 

  15. Van Doormaal BJ, Kistemaker G, Fatehi J, Miglior F, Jamrozik J, Schaeffer LR. Genetic evaluation of female fertility in Canadian dairy breeds. Interbull Bull. 2004;32:86–9.

    Google Scholar 

  16. Bouquet A, Juga J. Integrating genomic selection into dairy cattle breeding programmes: a review. Animal. 2013;7:705–13.

    Article  CAS  PubMed  Google Scholar 

  17. Schaeffer LR. Strategy for applying genome-wide selection in dairy cattle. J Anim Breed Genet. 2006;123:218–23.

    Article  CAS  PubMed  Google Scholar 

  18. Van Doormaal BJ. Impact of genomics on genetic selection and gain. Canada: Canadian Dairy Network; 2014. https://www.cdn.ca/document.php?id=340.

  19. Matukumalli LK, Lawley CT, Schnabel RD, Taylor JF, Allan MF, Heaton MP, et al. Development and characterization of a high-density SNP genotyping assay for cattle. PLoS One. 2009;4:e5350.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Purfield DC, Bradley DG, Kearney JF, Berry DP. Genome-wide association study for calving traits in Holstein-Friesian dairy cattle. Animal. 2014;8:224–35.

    Article  CAS  PubMed  Google Scholar 

  21. Cole JB, Wiggans GR, Ma L, Sonstegard TS, Lawlor TJ Jr, Crooker BA, et al. Genome-wide association analysis of thirty one production, health, reproduction and body conformation traits in contemporary U.S. Holstein cows. BMC Genomics. 2011;12:408.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Daetwyler HD, Capitan A, Pausch H, Stothard P, van Binsbergen R, Brondum RF, et al. Whole-genome sequencing of 234 bulls facilitates mapping of monogenic and complex traits in cattle. Nat Genet. 2014;46:858–65.

    Article  CAS  PubMed  Google Scholar 

  23. Erbe M, Hayes BJ, Matukumalli LK, Goswami S, Bowman PJ, Reich CM, et al. Improving accuracy of genomic predictions within and between dairy cattle breeds with imputed high-density single nucleotide polymorphism panels. J Dairy Sci. 2012;95:4114–29.

    Article  CAS  PubMed  Google Scholar 

  24. VanRaden PM, Null DJ, Sargolzaei M, Wiggans GR, Tooker ME, Cole JB, et al. Genomic imputation and evaluation using high-density Holstein genotypes. J Dairy Sci. 2013;96:668–78.

    Article  CAS  PubMed  Google Scholar 

  25. Sewalem A, Miglior F, Kistemaker GJ, Sullivan P, Van Doormaal BJ. Relationship between reproduction traits and functional longevity in Canadian dairy cattle. J Dairy Sci. 2008;91:1660–8.

    Article  CAS  PubMed  Google Scholar 

  26. Van Doormaal BJ, Kistemaker GJ, Miglior F. Implementation of reproductive performance genetic evaluations in Canada. Interbull Bull. 2007;37:129–33.

    Google Scholar 

  27. Miglior F, Muir BL, Van Doormaal BJ. Selection indices in Holstein cattle of various countries. J Dairy Sci. 2005;88:1255–63.

    Article  CAS  PubMed  Google Scholar 

  28. Sargolzaei M, Chesnais JP, Schenkel FS. A new approach for efficient genotype imputation using information from relatives. BMC Genomics. 2014;15:478.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Browning BL, Browning SR. A unified approach to genotype imputation and haplotype-phase inference for large data sets of trios and unrelated individuals. Am J Hum Genet. 2009;84:210–23.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81:559–75.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. VanRaden PM, Van Tassell CP, Wiggans GR, Sonstegard TS, Schnabel RD, Taylor JF, et al. Invited review: reliability of genomic predictions for North American Holstein bulls. J Dairy Sci. 2009;92:16–24.

    Article  CAS  PubMed  Google Scholar 

  32. MacCluer JW, Boyce AJ, Dyke B, Weitkamp LR, Pfenning DW, Parsons CJ. Inbreeding and pedigree structure in Standardbred horses. J Heredity. 1983;74:394–9.

    Article  Google Scholar 

  33. Sargolzaei M, Iwaisaki H, Colleau J. CFC: a tool for monitoring genetic diversity. In: Proceedings of the 8th world congress on genetics applied to livestock production, 13–18 August 2006; Belo Horizonte; 2006.

  34. Sargolzaei M. snp1101 User’s Guide; 2014.

  35. VanRaden PM. Efficient methods to compute genomic predictions. J Dairy Sci. 2008;91:4414–23.

    Article  CAS  PubMed  Google Scholar 

  36. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B. 1995;57:289–300.

    Google Scholar 

  37. Purfield DC, Bradley DG, Evans RD, Kearney FJ, Berry DP. Genome-wide association study for calving performance using high-density genotypes in dairy and beef cattle. Genet Sel Evol. 2015;47:47.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Hu ZL, Park CA, Wu XL, Reecy JM. Animal QTLdb: an improved database tool for livestock animal QTL/association data dissemination in the post-genome era. Nucleic Acids Res. 2013;41:D871–9.

    Article  CAS  PubMed  Google Scholar 

  39. Szkiba D, Kapun M, von Haeseler A, Gallach M. SNP2GO: functional analysis of genome-wide association studies. Genetics. 2014;197:285–9.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Grant JR, Arantes AS, Liao X, Stothard P. In-depth annotation of SNPs arising from resequencing projects using NGS-SNP. Bioinformatics. 2011;27:2300–1.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Sargolzaei M, Schenkel FS, Jansen GB, Schaeffer LR. Extent of linkage disequilibrium in Holstein cattle in North America. J Dairy Sci. 2008;91:2106–17.

    Article  CAS  PubMed  Google Scholar 

  42. da Huang W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4:44–57.

    Article  CAS  Google Scholar 

  43. Nejati-Javaremi A, Smith C, Gibson JP. Effect of total allelic relationship on accuracy of evaluation and response to selection. J Anim Sci. 1997;75:1738–45.

    Article  CAS  PubMed  Google Scholar 

  44. Hayes BJ, Visscher PM, Goddard ME. Increased accuracy of artificial selection by using the realized relationship matrix. Genet Res (Camb). 2009;91:47–60.

    Article  CAS  Google Scholar 

  45. Sargolzaei M, VanRaden PM, Kistemaker G, Schenkel F, Chesnais J. gebv User’s Guide, Release 1.0; 2013.

  46. Melka MG, Stachowicz K, Miglior F, Schenkel FS. Analyses of genetic diversity in five Canadian dairy breeds using pedigree data. J Anim Breed Genet. 2013;130:476–86.

    Article  CAS  PubMed  Google Scholar 

  47. Ma L, Wiggans GR, Wang S, Sonstegard TS, Yang J, Crooker BA, et al. Effect of sample stratification on dairy GWAS results. BMC Genomics. 2012;13:536.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Kadri NK, Guldbrandtsen B, Sørensen P, Sahana G. Comparison of genome-wide association methods in analyses of admixed populations with complex familial relationships. PLoS One. 2014;9:e88926.

    Article  PubMed  PubMed Central  Google Scholar 

  49. Marchini J, Howie B. Genotype imputation for genome-wide association studies. Nat Rev Genet. 2010;11:499–511.

    Article  CAS  PubMed  Google Scholar 

  50. Goddard ME. Uses of genomics in livestock agriculture. Anim Prod Sci. 2012;52:73–7.

    Article  CAS  Google Scholar 

  51. VanRaden PM, O’Connell JR, Wiggans GR, Weigel KA. Genomic evaluations with many more genotypes. Genet Sel Evol. 2011;43:10.

    Article  PubMed  PubMed Central  Google Scholar 

  52. de Roos APW, Hayes BJ, Goddard ME. Reliability of genomic predictions across multiple populations. Genetics. 2009;183:1545–53.

    Article  PubMed  PubMed Central  Google Scholar 

  53. Li Y, Willer C, Sanna S, Abecasis G. Genotype imputation. Annu Rev Genom Hum Genet. 2009;10:387–406.

    Article  CAS  Google Scholar 

  54. Segelke D, Chen J, Liu Z, Reinhardt F, Thaller G, Reents R. Reliability of genomic prediction for German Holsteins using imputed genotypes from low-density chips. J Dairy Sci. 2012;95:5403–11.

    Article  CAS  PubMed  Google Scholar 

  55. Wu C, DeWan A, Hoh J, Wang Z. A comparison of association methods correcting for population stratification in case–control studies. Ann Hum Genet. 2011;75:418–27.

    Article  PubMed  PubMed Central  Google Scholar 

  56. Price AL, Zaitlen NA, Reich D, Patterson N. New approaches to population stratification in genome-wide association studies. Nat Rev Genet. 2010;11:459–63.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Bermingham ML, Bishop SC, Woolliams JA, Pong-Wong R, Allen AR, McBride SH, et al. Genome-wide association study identifies novel loci associated with resistance to bovine tuberculosis. Heredity (Edinb). 2014;112:543–51.

    Article  CAS  Google Scholar 

  58. Sahana G, Guldbrandtsen B, Lund MS. Genome-wide association study for calving traits in Danish and Swedish Holstein cattle. J Dairy Sci. 2011;94:479–86.

    Article  CAS  PubMed  Google Scholar 

  59. Höglund JK, Guldbrandtsen B, Lund MS, Sahana G. Analyzes of genome-wide association follow-up study for calving traits in dairy cattle. BMC Genet. 2012;13:71.

    Article  PubMed  PubMed Central  Google Scholar 

  60. Saatchi M, Schnabel RD, Taylor JF, Garrick DJ. Large-effect pleiotropic or closely linked QTL segregate within and across ten US cattle breeds. BMC Genomics. 2014;15:442.

    Article  PubMed  PubMed Central  Google Scholar 

  61. Paliouras M, Diamandis EP. Androgens act synergistically to enhance estrogen-induced upregulation of human tissue kallikreins 10, 11, and 14 in breast cancer cells via a membrane bound androgen receptor. Mol Oncol. 2008;1:413–24.

    Article  PubMed  PubMed Central  Google Scholar 

  62. Ventolini G. Activation of simultaneous pathways in the initiation of parturition in humans. Med Hypotheses. 2013;1:11.

    Google Scholar 

  63. Storset AK, Slettedal IÖ, Williams JL, Law A, Dissen E. Natural killer cell receptors in cattle: a bovine killer cell immunoglobulin-like receptor multigene family contains members with divergent signaling motifs. Eur J Immunol. 2003;33:980–90.

    Article  CAS  PubMed  Google Scholar 

  64. Ha NT, Gross JJ, van Dorland A, Tetens J, Thaller G, Schlather M, et al. Gene-based mapping and pathway analysis of metabolic traits in dairy cows. PLoS One. 2015;10:e0122325.

    Article  PubMed  PubMed Central  Google Scholar 

  65. Ingvartsen KL. Feeding- and management-related diseases in the transition cow: physiological adaptations around calving and strategies to reduce feeding-related diseases. Anim Feed Sci Techn. 2006;126:175–213.

    Article  Google Scholar 

  66. Goff JP, Horst RL. Physiological changes at parturition and their relationship to metabolic disorders. J Dairy Sci. 1997;80:1260–8.

    Article  CAS  PubMed  Google Scholar 

  67. El-Tarabany MS. Effects of non-lactating period length on the subsequent calving ease and reproductive performance of Holstein, Brown Swiss and the crosses. Anim Reprod Sci. 2015;158:60–7.

    Article  PubMed  Google Scholar 

  68. Wehrli M, Cortinas-Elizondo F, Hlushchuk R, Daudel F, Villiger PM, Miescher S, et al. Human IgA Fc receptor FcαRI (CD89) triggers different forms of Neutrophil death depending on the inflammatory microenvironment. J Immunol. 2014;193:5649–59.

    Article  CAS  PubMed  Google Scholar 

  69. Bollopragada S, Youssef R, Jordan F, Greer I, Norman J, Nelson S. Term labor is associated with a core inflammatory response in human fetal membranes, myometrium, and cervix. Am J Obstet Gynecol. 2009;200:104.e1–11.

    Article  Google Scholar 

  70. Dowell P, Hu Z, Lane MD. Monitoring energy balance: Metabolites of fatty acid synthesis as hypothalamic sensors. Annu Rev Biochem. 2005;74:515–34.

    Article  CAS  PubMed  Google Scholar 

  71. Kaupe B, Brandt H, Prinzenberg EM, Erhardt G. Joint analysis of the influence of CYP11B1 and DGAT1 genetic variation on milk production, somatic cell score, conformation, reproduction, and productive lifespan in German Holstein cattle. J Anim Sci. 2007;85:11–21.

    Article  CAS  PubMed  Google Scholar 

  72. Kanehisa M, Goto S, Furumichi M, Tanabe M, Hirakawa M. KEGG for representation and analysis of molecular networks involving diseases and drugs. Nucleic Acids Res. 2010;38:D355–60.

    Article  CAS  PubMed  Google Scholar 

  73. Breuiller-Fouche M, Germain G. Gene and protein expression in the myometrium in pregnancy and labor. Reproduction. 2006;131:837–50.

    Article  CAS  PubMed  Google Scholar 

  74. Sawa A, Bogucki M, Krężel-Czopek Neja W. Association between rump score and course of parturition in cows. Anim Breed. 2013;56:816–22.

    Google Scholar 

  75. Atkins G. Using conformational anatomy to identify functionality in dairy cows. In: Proceedings of the 24th South Western Ontario Dairy Symposium, 22 February 2007; Woodstock; 2007. p. 37–48.

  76. Miyamoto T, Tsujimura A, Miyagawa Y, Koh E, Namiki M, Horikawa M, et al. Single-nucleotide polymorphisms in HORMAD1 may be a risk factor for azoospermia caused by meiotic arrest in Japanese patients. Asian J Androl. 2012;14:580–3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Gimpl G, Fahrenholz F. The oxytocin receptor system: structure, function, and regulation. Phys Rev. 2001;81:629–83.

    CAS  Google Scholar 

  78. Goddard M. Genomic selection: prediction of accuracy and maximisation of long term response. Genetica. 2009;136:245–57.

    Article  PubMed  Google Scholar 

  79. Goddard ME, Hayes BJ, Meuwissen THE. Using the genomic relationship matrix to predict the accuracy of genomic selection. J Anim Breed Genet. 2011;128:409–21.

    Article  CAS  PubMed  Google Scholar 

  80. Miglior F. Genetic evaluation of reproductive performance in Canadian dairy cattle. Ital J Anim Sci. 2010;6:29–37.

    Google Scholar 

Download references

Authors’ contributions

MKA designed the studies, performed the statistical analyses, and drafted the manuscript. MS helped to collect data, performed analyses and edited the manuscript. FSS, SPM, SSM, GP and PS, developed the project, designed the studies, helped perform analyses and edited the manuscript. DAG and SN helped perform analyses and edit the manuscript. LFB edited and helped to finalize the manuscript. All authors read and approved the final manuscript.

Acknowledgements

This research was funded by a Collaborative Research and Development grant from NSERC and L’Alliance Boviteq (LAB). We thank the Canadian Dairy Network for providing some of the data used in this study.

Competing interests

The authors declare that they have no competing interests.

Availability of data and materials

All the data supporting the results of this article are included within the article and in its supplementary files. The raw data cannot be made available, as it is property of the Canadian dairy cattle producers and this information is commercially sensitive.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Flavio S. Schenkel.

Additional files

12711_2017_356_MOESM1_ESM.docx

Additional file 1: Figure S1. Distribution of linkage disequilibrium (r2) calculated between pairs of SNPs on each chromosome before (A) and after (B) exclusion of the misplaced SNPs. The physical distances between pairs of SNPs are displayed along the horizontal axis, while the r2 value for each pair is displayed on the vertical axis. The mean r2 over successive intervals of 0.1 Mb is plotted in white

12711_2017_356_MOESM2_ESM.docx

Additional file 2: Figure S2. Distribution of 601,717 SNPs in the high-density panel across the bovine genome. Genomic coordinates of SNPs are displayed along the horizontal axis (Mb) and chromosome numbers are displayed on the vertical axis

12711_2017_356_MOESM3_ESM.docx

Additional file 3: Figure S3. The first three principal components of genetic co-ancestry based on Illumina BovineHD BeadChip (611,146 SNPs) genotypes for the 4848 bulls in the training population.

12711_2017_356_MOESM4_ESM.docx

Additional file 4: Figure S4. Genomic inflation factor (lambda) used for the test statistic (P-value) from the GWAS methods using the generalized linear mixed model accounting for reliability values associated with the de-regressed EBV. The traits are direct calf survival (cow) (SCSc); direct calf survival (heifer) (SCSh); direct calving ease (cow) (SCEc); direct calving ease (heifer) (SCEh); maternal calf survival (cow) (CSc); set of rear legs (SRL): calf survival (heifer) (CSh); daughter calving ability (DCA); heel depth (HD); calving ease (cow) (CEc); pin setting (PS); foot angle (FAN); calving ease (heifer) (CEh); leg rear view (LRV); median suspensory ligament (MSL); udder texture (UT); feet & legs (FL); rear attachment width (RAW); chest width (CW); rump (RU); rear attachment height (RAH); leg side view (LSV); mammary system (MS); loin strength (LS); angularity (ANG); height at front end (FE); conformation (CONF); fore udder attachment (FA); teat length (TL); rear teat placement (RTP); calving ability (CA); bone quality (BQ); front teat placement (FTP); body depth (BD); pin width (PW); dairy strength (DS); rump angle (RAN); udder depth (UD); and stature (ST).

12711_2017_356_MOESM5_ESM.xlsx

Additional file 5: Table S1. List of SNPs associated with calving performance and conformation traits at a FDR of 5% and located within the confidence interval of a previously reported quantitative trait loci (QTL).

12711_2017_356_MOESM6_ESM.xlsx

Additional file 6: Table S2. Gene ontology terms enriched at a FDR of 1% using SNP enrichment analysis for various calving performance and body conformation traits.

12711_2017_356_MOESM7_ESM.docx

Additional file 7: Figure S5. The oxytocin signaling pathway enriched for genes (in red boxes) located near the identified significant SNPs for calving performance and body conformation traits.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Abo-Ismail, M.K., Brito, L.F., Miller, S.P. et al. Genome-wide association studies and genomic prediction of breeding values for calving performance and body conformation traits in Holstein cattle. Genet Sel Evol 49, 82 (2017). https://doi.org/10.1186/s12711-017-0356-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12711-017-0356-8

Keywords