Implications
This study contributes to the understanding of the underlying genetic architecture and aetiology of dystocia and perinatal mortality both of which are quantitative traits that have a significant economic farm impact. Moreover, the detection of quantitative trait loci associated with only one calving trait can help resolve genetic antagonisms if targeted in genomic breeding programs – for example, selecting for decreased calving difficulty without impacting maternal calving difficulty.
Introduction
Dystocia and perinatal mortality are quantitative traits that significantly impact animal productivity and welfare (Bicalho et al., Reference Bicalho, Galvão, Cheong, Gilbert, Warnick and Guard2007; Bicalho et al., Reference Bicalho, Galvão, Warnick and Guard2008; Olsen et al., Reference Olsen, Hayes, Kent, Nome, Svendsen and Lien2010; Seidenspinner et al., Reference Seidenspinner, Tetens, Habier, Bennewitz and Thaller2011) and are common phenomena in both dairy and beef cattle (Hickey et al., Reference Hickey, Keane, Kenny, Cromie, Amer and Veerkamp2007; Mee et al., Reference Mee, Berry and Cromie2011). The contribution of dystocia and perinatal mortality to profit in cattle production systems is reflected by their inclusion in many national breeding programs. However, the genomic architecture influencing these traits has yet to be thoroughly elucidated.
Dystocia and perinatal mortality are influenced by both environmental factors such as age at first calving and season of calving (McGuirk et al., Reference McGuirk, Going and Gilmour1999; Berry and Cromie, Reference Berry and Cromie2009) and also by genetic factors (Berger et al., Reference Berger, Cubas, Koehler and Healey1992; Mee et al., Reference Mee, Berry and Cromie2011). Heritability estimates range from <0.01 to 0.17 for direct calving difficulty and <0.01 to 0.12 for the maternal components of calving difficulty (Steinbock et al., Reference Steinbock, Näsholm, Berglund, Johansson and Philipsson2003; Hickey et al., Reference Hickey, Keane, Kenny, Cromie, Amer and Veerkamp2007; Johanson et al., Reference Johanson, Berger, Tsuruta and Misztal2011). Heritability estimates for direct and maternal components of perinatal mortality are generally lower and range from <0.01 to 0.12 (Steinbock et al., Reference Steinbock, Näsholm, Berglund, Johansson and Philipsson2003; Johanson et al., Reference Johanson, Berger, Tsuruta and Misztal2011). Identification of genomic regions influencing lowly heritable traits typically requires a large, sometimes prohibitive, number of genotyped animals. Using however accurate measures of genetic merit of sires based on their large progeny group sizes, can partially overcome this requirement.
Previous quantitative trait loci (QTL) mapping and genome-wide association studies have identified several large genomic regions associated with calving difficulty and perinatal mortality. QTL regions on chromosomes 4, 6, 8, 12 and 18 have been associated with calving performance in several studies (Schrooten et al., Reference Schrooten, Bovenhuis, Coppieters and Van Arendonk2000; Kühn et al., Reference Kühn, Bennewitz, Reinsch, Xu, Thomsen, Looft, Brockmann, Schwerin, Weimann, Hiendleder, Erhardt, Medjugorac, Förster, Brenig, Reinhardt, Reents, Russ, Averdunk, Blümel and Kalm2003; Olsen et al., Reference Olsen, Hayes, Kent, Nome, Svendsen and Lien2010; Sahana et al., Reference Sahana, Guldbrandtsen and Lund2011). However, no sole causative QTL region seems to influence these traits, rather it is the cumulative effect of many of these regions and the allelic variation that lies within these that contributes to the observed phenotype.
The objective of this study was to identify regions of the bovine genome associated with dystocia and perinatal mortality in Holstein–Friesian dairy cows. Of particular interest is regions that are associated with only one trait thereby facilitating more targeted breeding that will help resolve genetic antagonisms, which are known to exist; for example, between direct and maternal calving difficulty (Steinbock et al., Reference Steinbock, Näsholm, Berglund, Johansson and Philipsson2003).
Material and methods
Phenotypic data
Predicted transmitting abilities (PTAs) and their associated reliabilities for 4683 Holstein–Friesian artificial insemination (AI) sires were obtained from the Irish national genetic evaluation in December 2012 for the direct and maternal component of calving difficulty and the direct component of perinatal mortality. In Ireland, calving difficulty is reported by producers on a subjectively scored linear scale of 1 to 4, where 1=no calving assistance; 2=slight assistance (assistance by one person, without needing to use a calf puller); 3=considerable assistance (assistance by one person using a calf puller or more than one person); 4=veterinary assistance (including caesarean). Perinatal mortality is recorded as a binary variable by producers, indicating whether the calf was dead within a 24 h period.
Parental contribution to the reliability of each PTA was removed using the method described by Harris and Johnson (Reference Harris and Johnson1998):
$$\tilde{R}={{R_{{\rm TRAD}} {\minus} R_{{\rm PA}} } \over {\left( {R_{{\rm TRAD}} {\times} R_{{\rm PA }} } \right){\plus}1{\minus}2R_{{\rm PA}} }}$$
where
$\tilde{R}$
is the reliability less the parental contribution, R
TRAD is the reliability from the traditional genetic evaluation from all relatives and R
PA is the parental average reliability. Only sires with an adjusted reliability of >40% for each trait were included in the analysis. Deregressed PTAs
$( \tilde y} )$
were obtained as
$$\tilde{y}=\dot{R}(\dot{R}^{{{\minus}1}} {\plus}A^{{{\minus}1}} )\hat{a}$$
where
$\hat{a}$
is a vector of PTAs from the domestic Irish national genetic evaluation, A is the numerator relationship of the matrix and
$\dot{R}$
is a diagonal matrix where the diagonal elements is
$$\dot{R}={{1{\minus}\tilde{R}} \over {\tilde{R}}}$$
where
$\tilde{R}$
is the animalʼs reliability less parental contribution. The relationship matrix was generated by tracing the pedigree of each animal back at least five generations where available. After edits, 1970 animals remained for direct calving dystocia, 1555 for maternal calving difficulty and 740 animals for direct perinatal mortality. The heritability of direct calving difficulty was 0.11 (s.e. 0.0137), maternal calving difficulty was 0.03 (s.e 0.0126) and perinatal mortality was 0.02 (s.e. 0.0156) (R. Evans, ICBF, personal communication, April 2013).
Genotypic data
Single nucleotide polymorphisms (SNPs) genotypes consisting of 54 001 SNPs from the BovineSNP50 BeadChip (Illumina In., San Diego, CA, USA) were available on the Holstein–Friesian AI bulls. Animals with genotype call rates <90% were removed. SNPs with >0.5% mendelian inconsistencies between parents and offspring, with call rates <95%, monomorphic SNPs and SNPs that deviated from Hardy–Weinberg equilibrium (P<1×10−9) were all discarded. SNPs also had to have a GCscore >0.15 and a GTscore >0.55 where GCscore and GTscore are quality measures of the genotype calls from the genotyping assay. All SNPs with unknown position and those located on the X chromosome were discarded. After edits a total of 43 204 SNPs remained across 1970 animals. Missing genotypes were imputed using Beagle Version 3.1.0 (Browning and Browning, 2007). Default parameters were used which included 10 iterations to impute missing data.
Whole genome association
Association analyses were undertaken by fitting all SNP simultaneously in a Bayesian framework (Habier et al., Reference Habier, Fernando, Kizilkaya and Garrick2011). The Bayesian approach was a two-step procedure: (1) Bayes-Cπ (Habier et al., Reference Habier, Fernando, Kizilkaya and Garrick2011) followed by (2) Bayes B (Meuwissen et al., Reference Meuwissen, Hayes and Goddard2001). Bayes B is sensitive to the previous genetic variance and an inaccurate estimate could affect the results. However, Bayes Cπ is more tolerant to the prior genetic variance and can provide an accurate estimation of this for prior use in the Bayes B model. The general statistical model can be written as
$$y=X\beta {\plus}u{\plus}\mathop{\sum}\limits_{K=1}^K {z_{k} a_{k} {\plus}e}$$
where y is a vector of phenotypes, X is the incidence matrix of the fixed effects, β is the vector of fixed effects, K is the number of SNPs, z k is the column vector representing the genotypes at SNP k, a k is the additive effect of that SNP and e is a vector of residual effects (Habier et al., Reference Habier, Fernando, Kizilkaya and Garrick2011). Both algorithms were undertaken using GenSel, a web-based program (http://bigs.ansci.iastate.edu/bigsgui/login.html) developed by Fernando and Garrick (Reference Fernando and Garrick2011). The chain length for both Bayesian analyses was 50 000 iterations, with the first 10 000 excluded as burn in. Deregressed PTA phenotypes were weighted by the respective reliability of the PTA.
Bayes Cπ was initially used to estimate π, the probability of a SNP having no association with the phenotype under investigation. As proposed by Habier et al. (Reference Habier, Fernando, Kizilkaya and Garrick2011), the priors of all SNP effects were assumed to have a common variance and the effect of an SNP fitted with probability (1−π) followed a mixture of multivariate Student’s t-distributions. Following the estimation of the posterior mean of π from the Bayes Cπ algorithm, the Bayes B algorithm was then invoked using the posterior mean of π. This Bayes B method fits a different variance for every covariate (SNP) in the model and is dependent on the prior (π).
Summary statistics of SNPs described hereafter are from the Bayes B analyses. The proportion of post-burn-in iterations that included the SNP in the model was calculated as the model frequency and this was used as evidence for an association. Bayes factors were used to quantify the strength of the posterior QTL probabilities for all loci (Kass and Raftery, Reference Kass and Raftery1995). Bayes factors can indicate the strength of an association based on their range of values; a Bayes factor >3.1 is indicative of ‘substantial evidence’ that the SNP is associated with a QTL, whereas Bayes factors >10.1 indicate ‘strong evidence’ and values >30.1 indicate ‘very strong evidence’ for association with a QTL (Jeffreys, Reference Jeffreys1961). The Bayes factor, BF, for comparing opposing hypotheses, was calculated as
$$BF={{{{\Pr (H_{1} \,\mid\,y)} \over {\Pr (H_{2} \,\mid\,y)}}} \mathord{\left/ {\vphantom {{{{\Pr (H_{1} \,\mid\,y)} \over {\Pr (H_{2} \,\mid\,y)}}} {{{\Pr (H_{1} )} \over {\Pr (H_{2} )}}}}} \right. \kern-\nulldelimiterspace} {{{\Pr (H_{1} )} \over {\Pr (H_{2} )}}}}$$
where y is the data, H 1 is the hypothesis that the marker is linked to a QTL and H 2 is the hypothesis that the marker is not linked to the QTL. The prior probability of the first hypothesis and alternative hypothesis was Pr(H 1) and Pr(H 2), respectively, and Pr(H 1|y) and Pr(H 2|y) were the posterior probability of H 1 and H 2. This can also be expressed in the simpler version of
$$BF={{{{\Pr (H_{1} \,\mid\,y)} \over {1{\minus}\Pr (H_{1} \,\mid\,y)}}} \mathord{\left/ {\vphantom {{{{\Pr (H_{1} \,\mid\,y)} \over {1{\minus}\Pr (H_{1} \,\mid\,y)}}} {{{\Pr (H_{1} )} \over {1{\minus}\Pr (H_{1} )}}}}} \right. \kern-\nulldelimiterspace} {{{\Pr (H_{1} )} \over {1{\minus}\Pr (H_{1} )}}}}$$
The Pr(H 1) was the posterior mean of (1−π) from the Bayes Cπ analysis, which was used as the prior for the Bayes B analysis.
The posterior probability of association (PPA) can be thought of as the Bayesian analogue of a P-value and this was calculated for all SNPs (Stephens and Balding, Reference Stephens and Balding2009). The PPA combines the evidence in the observed data (i.e. the BF) with the prior probability that the SNP is truly associated with the SNP. To calculate the PPA, the posterior odds (PO) was calculated first, where π was the value calculated from the Bayes Cπ analyses and therefore 1−π is the proportion of SNPs that are believed to be associated with the trait:
$${\rm PO}=\left( {{\rm BF}{\times}\left( {{\rm 1}{\minus}\pi } \right)} \right)/\left( {{\rm 1}{\minus}\left( {{\rm 1}{\minus}\pi } \right)} \right)$$
and the PPA was subsequently calculated as
$${\rm PPA}={\rm PO}/\left( {{\rm 1}{\plus}{\rm PO}} \right)$$
The closer the PPA value is to one the more convincing the evidence for an association.
Beacuse of the existence of genomic regions with strong linkage disequilibrium (LD) in the bovine genome (McKay et al., Reference McKay, Schnabel, Murdoch, Matukumalli, Aerts, Coppieters, Crews, Dias Neto, Gill, Gao, Mannen, Stothard, Wang, Van Tassell, Williams, Taylor and Moore2007), the effect of any particular QTL and the SNP model frequency may be distributed across numerous SNPs that are in LD with the QTL. This may hinder the identification of true QTL (Fan et al., Reference Fan, Onteru, Du, Garrick, Stalder and Rothschild2011; Onteru et al., Reference Onteru, Fan, Nikkilä, Garrick, Stalder and Rothschild2011). To accurately identify the presence of QTL, a sliding genomic window of five consecutive SNPs based on genomic location was used to predict the genomic breeding value of the region. This was achieved by summing the contribution of the SNPs to the genetic variance per window. This approach also helps control against false-positives arising from single putative SNPs (Fan et al., Reference Fan, Onteru, Du, Garrick, Stalder and Rothschild2011; Onteru et al., Reference Onteru, Fan, Nikkilä, Garrick, Stalder and Rothschild2011).
LD blocks were constructed for a 500 kb region up and downstream surrounding each of the 10 most associated SNPs per trait. LD was estimated using the solid spine of LD option in Haploview, which is an internally developed method that searches for a ‘spine’ of strong LD running from one marker to another (Barrett et al., Reference Barrett, Fry, Maller and Daly2005). Gene search using ensemble (http://ensembl.org) on the Btau 4.0 genome build was completed for those regions that were located in LD blocks and also neighbouring LD blocks if they were side by side with the LD blocks containing the SNP of interest and strong LD was seen between SNPs. QTL regions were compared with previously reported QTL in CattleQTLdb for calving difficulty or perinatal mortality as well as likely underlying correlated traits such as birth weight and stature.
Genes within 500 kb of SNPs that showed ‘substantial evidence’ of an association (i.e. BF>3.1) were also identified using Ensembl Build BTA4.0 for all three calving performance traits. Biomart was then used to map the genes in each of these three data sets to their human orthologs using the mapping available from version hg19 of the human genome and GO terms were also identified. For each data set, the R package GOseq, without the correction for gene length bias, was used to identify the KEGG pathways, which were significantly over-represented by the set of genes (P<0.05) compared against a background of all genes in the human genome. Functional annotation of the genes was undertaken using DAVID (http://david.abcc.ncifcrf.gov).
Results
Direct calving difficulty
The mean proportion of SNPs estimated not to be associated with direct calving difficulty (i.e. π) was 0.94. The proportion of genetic variance in direct calving difficulty accounted for by all SNPs using the Bayes B algorithm was 0.93. The Bayesian posterior probabilities (model frequency) per SNP varied (Table 1) with one SNP included in the model 97.41% of the time (Figure 1). Table 2 details the Bayes factors, PPAs and minor allele frequency (MAF) for all 45 SNPs that showed at least 'substantial evidence' (Jeffreys, Reference Jeffreys1961) of an association with direct calving difficulty. The Bayes factor for the strongest SNP, ss86324977, was 580.21, suggesting ‘very strong evidence’ of a QTL (Jeffreys, Reference Jeffreys1961). The SNP ss86324977 accounted for 2.10% of the genetic variance in direct calving difficulty and the surrounding SNPs contributed no more to the genetic variance with the five SNP window containing ss86324977 accounting for 2.11% of the genetic variance. The distance between the flanking SNPs in the window varied either side of ss86324977, with only a distance of ∼30 kb between BTA-97501-no-rs and ss86324977 but ARS-BFGL-NGS-6380 was over 200 kb downstream from ss86324977. However, strong LD existed between these flanking SNPs, with a D′ value of one (r 2-value of 0.194) between BTA-97501-no-rs and ss86324977 and a D′ value of 0.959 (r 2-value of 0.035) between ss86324977 and ARS-BFGL-NGS-6380. The D′ was one between these SNPs due to only three of the four possible haplotypes existing in the population. LD analysis undertaken on a 500 kb region either side of ss86324977 revealed several small LD blocks (Figure 2). The SNP ss86324977 was located in an LD block with the SNP BTA-97501-no-rs. Moderate to strong levels of LD also existed between some of the SNPs in neighbouring LD blocks downstream of ss86324977, with a D′ value of 0.711 (r 2-value of 0.013) between ARS-BFGL-NGS-23169 and ARS-BFGL-NGS-100080. Over 83% of the study population were homozygous CC for the SNP ss86324977, whereas only 1% were homozygous AA. The allele substitution effect for the A allele was +0.576 units (Table 2). Animals that were homozygous CC had a lower mean PTA for direct calving difficulty (−3.20 units; s.e. 0.06), compared with heterozygote (−1.89 units; s.e. 0.08) and homozygous AA (−0.59 units; s.e. 0.08).
The Bayesian posterior probability of model frequency from a Bayes B analysis for calving traits. (a) Manhattan plot for direct calving difficulty; chromosome 18 had a single nucleotide polymorphism (SNP) with the highest model frequency appearing in the Bayesian model 97.4% of the time. (b) Manhattan plot for maternal calving difficulty; chromosome 28 had an SNP with a model frequency value 73.4%. (c) Manhattan plot for direct perinatal mortality; chromosome 14 had an SNP with a model frequency value 19.3%.
LD analysis of a 500 kb region up and downsteam surrounding ARS-BFGL-NGS-109285 (ss86324977) on chromosome 18. Seven LD blocks are outlined with block 4 containing ARS-BFGL-NGS109285, which was the most significantly associated SNP with direct calving difficulty. LD=linkage disequilibrium.
Proportion of genetic variance explained by SNPs and the number of SNPs found in three different Bayesian posterior probability ranges for all traits
SNP=single nucleotide polymorphism.
Bayes Factors and PPA for all SNPs that showed 'substantial evidence' (Bayes factor >3.1) of an association with direct calving difficulty ordered by PPA
PPA=posterior probability of associations; SNP=single nucleotide polymorphism.
LD structure upstream and downstream of the top 10 SNPs most often included in the model can be seen in Supplementary Figure S1 and genes identified within these LD blocks are listed in Supplementary Table S1. In total, over 20 genes were identified overlapping with all LD blocks. These genes ranged in function from transcription factors, transporters and membrane proteins. The ss86324977 SNP was located in an intron in the gene sialic acid binding Ig-like lectin 5 (SIGLEC5).
Gene annotation on all 45 SNPs that showed 'substantial evidence' for an association with direct calving difficulty (Table 2) revealed that 538 genes or gene products were identified within a 500 kb distance up and downstream. A total of 454 genes mapped to 428 human orthologs. The most significantly over represented KEGG pathway was the riboflavin metabolism signalling pathway (P=1.54×10−3). The protein binding GO term, GO:0005515, was the most associated term (12.6%) for genes associated with direct calving difficulty.
Maternal calving difficulty
The proportion of SNPs not associated with maternal calving diffculty (i.e. π) was 0.94. The proportion of the genetic variance in maternal calving difficulty accounted for by all SNPs was 0.96 (Table 1). The Bayes factors, PPAs and MAF for all 44 SNPs that showed 'substantial evidence' (Jeffreys, Reference Jeffreys1961) of an association with maternal calving difficulty are in Table 3. A moderate level of heterozygosity and moderate MAF (ranging from 0.13 to 0.50) existed for all 44 SNPs associated with maternal calving difficulty. Only one SNP on chromosome 28, rs41652463, exhibited ‘very strong evidence’ for association with maternal calving difficulty, whereas three of the remaining SNPs showed ‘strong evidence’ (Jeffreys, Reference Jeffreys1961; Figure 1). The SNP rs41652463 was included in the Bayesian model for 73% of the iterations, accounting for 0.69% of the genetic variance in maternal calving difficulty with the additional four flanking SNPs in the window containing rs41652463 contributing no additional genetic variance. LD analysis carried out on a 500 kb region flanking rs41652463 revealed that, although the SNP was not in a LD block, three surrounding LD blocks were assigned (Figure 3). Weak LD existed between SNPs flanking rs41652463, suggesting the possibility of a recombination hotspot but also that the causal mutation is likely to be very near the rs41652463 SNP.
LD analysis of a 500 kb region up and downstream surrounding BTA-63756-no-rs (rs41652463) on chromosome 28. BTA-63756-no-rs was found to have low LD with adjacent SNPs but is surrounded by two LD blocks in the vicinity. LD=linkage disequilibrium; SNP=single nucleotide polymorphism.
Bayes Factors and PPA for all SNPs that showed 'substantial evidence' (Bayes factor >3.1) of an association with maternal calving difficulty ordered by PPA
PPA=posterior probability of associations; SNP=single nucleotide polymorphism.
Two SNPs, ss61478550 and rs41628739 on chromosome 12 were in close proximity, just over 80 kb apart and were both associated with maternal calving difficulty; the LD D′ value between both SNPs was of 1 (r 2-value was 0.954). Combined, both of these SNPs were included in the model 84.38% of the time suggesting a QTL influencing maternal calving difficulty to be located in the vicinity. Both SNPs were located within the same LD block (Supplementary Figure S2), however, no genes were identified within this region but one microRNA did exist, bta-mir-1256, a stemloop microRNA involved in post-transcriptional regulation of gene expression. Several other genes were identified within LD blocks of the remaining seven most often included SNPs in the model and are detailed in Supplementary Table S2.
Gene annotation on the 44 SNPs that showed 'substantial evidence' for an association with maternal calving difficulty (Table 3) revealed that 313 genes were identified within a 500 kb distance up and downstream. A total of 290 genes mapped to 283 human orthologs. The glycosaminoglycan degradation pathway was the most significantly over represented KEGG pathway (P=1.57×10−3) and the protein binding GO term, GO:0005515, was the most associated term again (9.8%).
Perinatal mortality
The mean proportion of SNPs not associated with perinatal mortality (i.e. π) was 0.96. The proportion of genetic variance accounted for by all SNPs was 0.96. Bayesian posterior probabilities were low for perinatal mortality with the highest value of 0.19 for SNP ss86296129 on chromosome 14 (Figure 1). Bayes Factors and PPAs were consequently low for this trait with no SNP showing a Bayes Factor value >10.1, indicating there was no evidence for a ‘strong association’ (Jeffreys, Reference Jeffreys1961) for this trait (Table 4). Although chromosome 14 had the SNP with the highest posterior probability, accounting for 0.0287% of the genetic variance in perinatal mortality, a five-SNP window on chromosome 18 contributed the greatest to the genetic variance (i.e. 0.0338% of the genetic variation). Excluding the SNP on chromosome 18 (MAF of 0.07), a moderate to high level of heterozygosity and consequently high MAF (ranging from 0.17 to 0.49) existed for all 17 SNPs that showed ‘substantial evidence’ for an association with perinatal mortality. LD structure of the region surrounding the ss86296129 SNP on chromosome 14 revealed that it was located in a LD block (Figure 4). Strong LD existed between the SNP ss86296129 and both ARS-BFGL-NGS-8890 and BTA-110405-no-rs from the neighbouring LD block, with a D′ value of 1; therefore, gene annotation was undertaken on both LD blocks. No genes were in the immediate vicinity of ss86296129 with the nearest gene (SLC26A7) located >300 kb from this SNP. LD structure for the remaining nine SNPs most often included in the model can be seen in Supplementary Figure S3 and several genes identified within these LD blocks and their functions are listed in Supplementary Table S3.
LD analysis of a 500 kb region surrounding ARS-BFGL-NGS-16146 (ss86296129) on chromosome 14. Three LD blocks are outlined with block 2 containing ARS-BFGL-NGS-16146, which was the most significantly associated SNP with perinatal mortality. LD=linkage disequilibrium; SNP=single nucleotide polymorphism.
Bayes factors and PPA for the ten most significant SNPs associated with direct perinatal mortality ordered by PPA
PPA=posterior probability of associations; SNP=single nucleotide polymorphism.
A total of 200 genes were identified within a 500 kb distance up and downstream for all SNPs that showed ‘substantial evidence’ of an association with perinatal mortality. A total of 164 genes mapped to 146 human orthologs. The most significantly over-represented KEGG pathway was the systemic lupus erythematosus pathway (P=1.64×10−5) and again the protein binding GO term, GO:0005515, was the most associated gene ontology term (9.6%).
Discussion
In the present study, we identified several SNPs that were ‘strongly associated’ (Jeffreys, Reference Jeffreys1961) with calving performance traits in Irish dairy cows. Associations existed on 16 different chromosomes, some of which have previously been documented to be associated with calving traits (i.e. ss86324977), whereas others are novel previously unreported associations but have been reported to be associated with growth, skeletal development, birth weight and stature (Cole et al., Reference Cole, VanRaden, O'Connell, Van Tassell, Sonstegard, Schnabel, Taylor and Wiggans2009; McClure et al., 2010; Meyers et al., Reference Meyers, McDaneld, Swist, Marron, Steffen, O'Toole, O'Connell, Beever, Sonstegard and Smith2010). QTL regions that were solely associated with one trait were also identified facilitating direct genomic selection for these regions thereby possibly aiding in resolving known genetic antagonisms between the calving performance traits (Luo et al., 1999; Steinbock et al., Reference Steinbock, Näsholm, Berglund, Johansson and Philipsson2003). However, caution has to be taken when interpreting these QTL regions, as the power of detecting QTLs associated with perinatal mortality in comparison with the other traits was reduced because of the lower quantity of records.
Chromosome 18
The SNP, ss86324977 on chromosome 18, which was most strongly associated with direct calving difficulty in the present study, has previously been documented to be associated with calving difficulty in United States (Cole et al., Reference Cole, VanRaden, O'Connell, Van Tassell, Sonstegard, Schnabel, Taylor and Wiggans2009) and both Danish and Swedish (Sahana et al., Reference Sahana, Guldbrandtsen and Lund2011; Höglund et al., Reference Höglund, Guldbrandtsen, Lund and Sahana2012) dairy cattle populations. This SNP is located in an intron of the sialic acid binding Ig-like lectin 5 gene (Siglec-5). Cole et al. (Reference Cole, VanRaden, O'Connell, Van Tassell, Sonstegard, Schnabel, Taylor and Wiggans2009) reported this gene to be associated with a variety of calving traits such as sire and daughter calving ease, rump width, stature and body depth. In humans, a similar gene Siglec-6 is expressed in the placenta and may play a role in the initiation of parturition (Brinkman-Van der Linden et al., Reference Brinkman-Van der Linden, Hurtado-Ziola, Hayakawa, Wiggleton, Benirschke, Varki and Varki2007). Cole et al. (Reference Cole, VanRaden, O'Connell, Van Tassell, Sonstegard, Schnabel, Taylor and Wiggans2009) hypothesized that animals homozygous ‘AA’ for the SNP ss86324977 may experience a delayed parturition because of a lectin deficiency caused by high levels of sialic acid-binding Ig-like lectin 5. The high frequency in the present population of the favourable allele associated with reduced difficulty calving implies past breeding programs have selected against the AA genotype of this SNP or that it has been subjected to natural selection pressure. Qanbari et al. (Reference Qanbari, Gianola, Hayes, Schenkel, Miller, Moore, Thaller and Simianer2011) confirmed the presence of a strong selection signature using both iHS and F ST analyses in the vicinity of the sialic acid-binding Ig-like lectin 5 gene. Qanbari et al. (Reference Qanbari, Gianola, Hayes, Schenkel, Miller, Moore, Thaller and Simianer2011) suggested that this selection signature is the outcome of strong and recent selective pressure because of the long run of homozygosity observed in this region. This SNP ss86324977 was also shown to be associated with perinatal mortality in the present study thereby substantiating the known genetic association between both traits (Luo et al.,1999); the allele substitution effect for the A allele of this SNP was positive for both direct calving difficulty and perinatal mortality. This suggests this genomic region harbours either a gene with pleiotrophic effects on both traits or harbours two or more linked (i.e. co-inherited) genes affecting both traits.
Chromosome 17
Four of the top 10 SNPs for direct and maternal calving difficulty were located on chromosome 17. The SNP ss117972188, the second most associated SNP with direct calving difficulty, was located in the gene HHIP whose orthologue has been shown to be associated with human stature (Liu et al., Reference Liu, Medland, Wright, Henders, Heath, Madden, Duncan, Montgomery, Martin and McRae2010). However, this HHIP gene has not previously been identified as a significant SNP controlling stature in cattle (Pryce et al., Reference Pryce, Hayes, Bolormaa and Goddard2011). Stature has been previously correlated with the direct effects of calving difficulty, with larger calves more likely to experience a difficult birth (Bennett and Gregory, Reference Bennett and Gregory2001). The remaining three SNPs (ss86289496, ss86307824 and ss105262945) have not previously been documented to be associated with calving traits in cattle. However, SNP ss86289496 was located in close proximity to zinc finger protein 827 and these proteins, which are well-known transcription factors, have recently been reported to play an important role in the specific regulation of mammalian bone and skeletal development with zinc deficiencies known to result in skeletal growth retardation (Ganss and Jheon, Reference Ganss and Jheon2004); this therefore is a likely positional candidate gene.
The allele effects for these four SNPs (ss117972188, ss86289496, ss86307824 and ss105262945) corroborate the negative genetic correlation that is known to exist between direct and maternal calving difficulty (Steinbock et al., Reference Steinbock, Näsholm, Berglund, Johansson and Philipsson2003) with the allele effects of ss117972188 and ss86289496 having a positive effect on direct calving difficulty and a negative effect on maternal calving difficulty and vice versa for SNP ss86307824.
Chromosomes 2, 16 and 12
Birth weight is known to be non-linearly associated with dystocia and perinatal mortality, with lighter and heavier than average calves tending to have an increased risk of mortality and dystocia (Berger et al., Reference Berger, Cubas, Koehler and Healey1992). Putative QTL regions affecting birth weight are well documented on chromosome 2; Grosz and MacNeil (Reference Grosz and MacNeil2001) reported that ∼2.8 kg difference of birth weight could be accounted for by a QTL residing at the 114 cM region. In the present study, the SNP ss86303585 on chromosome 2 associated with perinatal mortality was in close proximity to this QTL region described by Grosz and MacNeil (Reference Grosz and MacNeil2001), suggesting this SNP may be in a region influencing birth weight and therefore perinatal mortality. Gene analysis also revealed that the insulin-like growth factor-binding protein 2 precursor gene (IGFBP-2) was located just 140 kb downstream of the ss86316303 SNP. These IGF-binding proteins inhibit insulin like growth factors with IGFBP-2 believed to play an important role as a growth-promoting hormone during gestation (Fowden, Reference Fowden2003) therefore making it a logical candidate gene. McClure et al. (2010) previously associated several QTL regions with birth weight, two of which overlapped with regions identified in the present study; the QTL region around 72 to 74 Mb on chromosome 16 associated with all three traits and the 57 Mb QTL region on chromosome 12 associated only with maternal calving difficulty.
Chromosomes 23 and 28
Evidence of a strong QTL for maternal calving difficulty containing the SNP rs41652463 existed on chromosome 28, in close proximity to a QTL documented previously to be associated with direct but not maternal calving difficulty by Olsen et al. (Reference Olsen, Hayes, Kent, Nome, Svendsen and Lien2010). Similarly, a QTL associated with direct calving difficulty on chromosome 23 in the present study was located 1 Mb downstream of a QTL documented in Danish and Swedish cattle populations (Sahana et al., Reference Sahana, Guldbrandtsen and Lund2011; Höglund et al., Reference Höglund, Guldbrandtsen, Lund and Sahana2012) that influenced maternal calving traits such as the calving index, maternal calving ease and maternal still birth in first-parity animals. An earlier study by Seidenspinner et al. (Reference Seidenspinner, Bennewitz, Reinhardt and Thaller2009) also identified a QTL influencing direct calving difficulty within this region. This suggests QTL influencing calving traits are located within these genomic regions on chromosomes 23 and 28. The QTL on chromosome 28 was not located in an LD block but several genes associated with the nucleus and transportation were located within a 500 kb distance with a small nucleolar RNA belonging to the SNORA36 family and the Graves disease carrier protein being in the closest proximity. Similarly on chromosome 23, no genes were identified in the LD block containing SNP ss86328687, however, several genes associated with mRNA enzymes and triggering receptors on myeloid cells were identified within 500 kb of the QTL region surrounding both ss86328687 and ss896339077.
Chromosomes 10 and 14
Three possible candidate genes (SLC26A7, TMEME55A and TMEME64) were identified within the LD block containing ss86296129 on chromosome 14, the SNP most strongly associated with perinatal mortality. Two of these genes, SLC26A7 and TMEM64 have been documented to be associated with cartilage and bone development in humans (Lohi et al., Reference Lohi, Kujala, Makela, Lehtonen, Kestila, Saarialho-Kere, Markovich and Kere2002; Kim et al., Reference Kim, Kim, Jeong, Cho, Han, Takegahara, Negishi-Koga, Takayanagi, Lee, Sul, Prasad, Lee and Choi2013) suggesting that these genes may be responsible for the many morphological abnormalities associated with perinatal mortality. In cattle, a deletion in a similar gene (SLC4A2) has been shown to be associated with osteoporosis in Red Angus cattle (Meyers et al., Reference Meyers, McDaneld, Swist, Marron, Steffen, O'Toole, O'Connell, Beever, Sonstegard and Smith2010). Osteoporosis is a bone disorder observed in humans and many other animals and can be inherited as a recessive defect. Loss of SLC4A2 function induces premature cell death and cytoplasmic alkalinization of osteoclasts (Meyers et al., Reference Meyers, McDaneld, Swist, Marron, Steffen, O'Toole, O'Connell, Beever, Sonstegard and Smith2010). The death of osteoclasts, which are large multinucleated cells that reabsorb bone, can result in the formation of extremely dense fragile bones and affected calves are typically born stillborn, slightly premature and morphologically small with a flat skull (Meyers et al., Reference Meyers, McDaneld, Swist, Marron, Steffen, O'Toole, O'Connell, Beever, Sonstegard and Smith2010). Selecting on the causative mutation tagged by the ss86296129 SNP could aid in reducing both direct calving difficulty and perinatal mortality; a moderate MAF exists for this SNP signifying scope for selection.
The two QTL on chromosome 10 for perinatal mortality were not in the vicinity of QTL previously reported to be associated with perinatal mortality on this chromosome (Kühn et al., Reference Kühn, Bennewitz, Reinsch, Xu, Thomsen, Looft, Brockmann, Schwerin, Weimann, Hiendleder, Erhardt, Medjugorac, Förster, Brenig, Reinhardt, Reents, Russ, Averdunk, Blümel and Kalm2003; Schnabel et al., Reference Schnabel, Sonstegard, Taylor and Ashwell2005; Sahana et al., Reference Sahana, Guldbrandtsen and Lund2011; Seidenspinner et al., Reference Seidenspinner, Tetens, Habier, Bennewitz and Thaller2011; Höglund et al., Reference Höglund, Guldbrandtsen, Lund and Sahana2012). Two candidate genes (RAB11A and MEGF11) were identified in the LD block of rs43616983. The gene RAB11A that encodes the Rab11A small GTPase, has been shown to play an important role in mediating transferrin recycling back to the plasma membrane and may play a role in innate immunity (Ullrich et al., Reference Ullrich, Reinsch, Urbé, Zerial and Parton1996). A weak innate immune system has previously been documented to contribute to perinatal mortality (Kirkbride, Reference Kirkbride1993).
Pathway analysis
Several KEGG pathways were associated with calving performance traits of which the systemic lupus erythematosus pathway was the most overrepresented for perinatal mortality. Systemic lupus erythematosus is an autoimmune disease that affects connective tissues, suggesting that the genes associated with perinatal mortality are primarily of an immune function. Systemic lupus erythematosus has been associated with perinatal mortality in humans where an increased risk of spontaneous abortion, intrauterine foetal death, preclampsia, intrauterine growth retardation and preterm birth has been documented (Clowse et al., Reference Clowse, Magder, Witter and Petri2005). The riboflavin metabolism pathway was the most overrepresented pathway for direct calving difficulty. Riboflavin, an important micronutrient, is required for the metabolism of fats, ketone bodies, carbohydrates and proteins, and consequently has been shown to impact growth (Chan et al., Reference Chan, Deng, Mikael, Yan, Pickell, Wu, Caudill and Rozen2010). This impact on growth may be why it was the most overrepresented KEGG pathway for this trait.
Glycosaminoglycans have been documented to be associated with the reproductive cycle in many species. Their composition and concentration vary at different stages throughout the cycle and they have been documented to be important regulators of cervical function as they contribute to the visco-elastic properties of the connective tissue present there (Osmers et al., Reference Osmers, Rath, Pflanz, Kuhn, Stuhlastz and Szeverényi1993). The glycosaminoglycan degradation pathway was the most overrepresented pathway for maternal calving difficulty possibly because of the continuous role they play throughout pregnancy particularly with cervical ripening and dilatation being associated with vast changes in glycosaminoglycan pattern (Osmers et al., Reference Osmers, Rath, Pflanz, Kuhn, Stuhlastz and Szeverényi1993).
Conclusion
This study confirmed several different QTL regions previously documented to be associated with calving traits but also several novel QTL regions. Identification of QTLs that are only associated with one calving trait have also been identified allowing for direct selection through genomic breeding programs while overcoming the genetic antagonism that exists between calving traits. Two genes of particular interest identified in this study, Siglec-5 and SLC26A7, have themselves or similarly been shown to impact calving traits, suggesting that there is the existence of influential calving trait genes.
Acknowledgements
This work was supported by Science Foundation Ireland principal investigator award grant number 09/IN.1/B2642 and Research Stimulus Fund 11/5/112.
