Implications
Traditional genetic selection for fertility is made difficult by the low heritability of traditional fertility phenotypes. This limitation can be minimised by increasing the heritability of the trait(s) or by using genomic information that explains a large proportion of the genetic variation in fertility. In the present study, both approaches were investigated. Heritability estimates greater than some traditional fertility traits were observed for the post-partum commencement of luteal activity quantified using milk progesterone concentration. Specific regions of the genome associated with fertility were also observed, albeit they only explained a small proportion of the total genetic variance in the traits.
Introduction
The well-documented deterioration in reproductive performance in international populations of Holstein–Friesian dairy cattle (Royal et al., Reference Royal, Darwash, Flint, Webb, Woolliams and Lamming2000) has led to the broadening of international breeding goals to now include functional traits such as fertility (Miglior et al., Reference Miglior, Muir and Van Doormaal2005). However, the low heritability of many traditional fertility traits (Pryce and Veerkamp, Reference Pryce and Veerkamp2001; Berry et al., Reference Berry, Buckley, Dillon, Evans, Rath and Veerkamp2003; Wall et al., Reference Wall, Brotherstone, Woolliams, Banos and Coffey2003), coupled with the long time period necessary to acquire some of the phenotypes (e.g. calving interval, CI), imply that selection for these traits using traditional methods is more difficult than selection, for example, for milk yield traits. Fertility phenotypes defined using information on hormonal profiles are, in general, more heritable (Veerkamp et al., Reference Veerkamp, Oldenbroek, Van Der Gaast and Van Der Werf2000; Royal et al., Reference Royal, Flint and Woolliams2002; Petersson et al., Reference Petersson, Berglund, Strandberg, Gustafsson, Flint, Wooliams and Royal2007), partly because of a possible reduction in data recording errors or minimising the impact of preferential treatment of some animals such as post-partum voluntary waiting periods to service. The generation of such hormonal profiles on sufficient animals for routine use in genetic evaluations is, nonetheless, currently not economically feasible. Use of genomic information in genetic evaluation of animals, however, has the potential to overcome some of the limitations of traditional methods of genetic evaluations.
The recent availability of fast-throughput genotyping platforms with dense genome-wide markers is reducing the cost of genotyping considerably, resulting in accurate phenotypes on large numbers of animals now being the limiting factor in association analyses. Most studies to date that attempted to identify quantitative trait loci (QTLs) in dairy cattle have used specific study designs such as (grand-) daughter designs (Weller et al., Reference Weller, Kashi and Soller1990; Holmberg and Andersson-Eklund, Reference Holmberg and Andersson-Eklund2006; Höglund et al., Reference Höglund, Buitenhuis, Guldbrandtsen, Su, Thomsen and Lund2009), and have therefore been generally limited to routinely available phenotypes. We are not aware of any genome-wide association study for traditional measures of fertility in dairy cattle that has been undertaken using actual cow phenotypic data. Furthermore, we are not aware of any study that has attempted to locate genomic regions associated with fertility phenotypes derived from progesterone profiles, most likely because of the expense associated with collecting such data on a sufficiently large sample population. In this study, we collate cow genotypic and phenotypic information from research herds in Ireland, the United Kingdom, the Netherlands and Sweden and relate genetic markers across the genome to the fertility phenotypes using a Bayesian stochastic search variable selection statistical approach.
Material and methods
The data were from Teagasc, Moorepark, Ireland; Scottish Agricultural College, United Kingdom; Wageningen UR Livestock Research, the Netherlands and the Jälla herd of the Swedish University of Agricultural Science. In the present study, these datasets will be referred to as ‘Ireland’, ‘UK’, ‘the Netherlands’ and ‘Sweden’, respectively. Phenotypic data were available on 2031 Irish, 1018 UK, 725 Dutch and 225 Swedish primiparous Holstein cows. A detailed description of the experimental treatments imposed on the animals in the different countries is provided elsewhere (Veerkamp et al., Reference Veerkamp, Oldenbroek, Van Der Gaast and Van Der Werf2000; Horan et al., Reference Horan, Mee, O'Connor, Rath and Dillon2005; Petersson et al., Reference Petersson, Strandberg, Gustafsson and Berglund2006; Pollot and Coffey, Reference Pollot and Coffey2008).
Sampling procedures for milk and quantification of the progesterone concentration in milk have previously been described in detail for the data originating from Ireland (Horan et al., Reference Horan, Mee, O'Connor, Rath and Dillon2005), the UK (Pollot and Coffey, Reference Pollot and Coffey2008), the Netherlands (Veerkamp et al., Reference Veerkamp, Oldenbroek, Van Der Gaast and Van Der Werf2000) and Sweden (Petersson et al., Reference Petersson, Strandberg, Gustafsson and Berglund2006). Milk sampling for progesterone was undertaken on Mondays, Wednesdays and Fridays in both Ireland and the UK, and data between the years 2001 and 2004 in Ireland and between 2003 and 2005 in the UK were included in this study. Between the years 1991 and 1998, inclusive milk sampling for progesterone in the Netherlands was undertaken on Mondays and Fridays and between 2003 and 2004, inclusive milk sampling was undertaken on Mondays and Thursdays. In Sweden, milk sampling for progesterone was undertaken twice weekly between 1989 and 2009. Quantification of milk progesterone concentration across countries was undertaken on raw milk.
Definition of fertility phenotypes
Two distinct types of fertility traits were defined: (1) traditional fertility traits and (2) a fertility trait derived from progesterone profiles. The traditional fertility traits in the present study were days from calving to first observed heat (CFH) or first service (CFS), calving interval (CI), number of services (NS) and pregnancy rate to first service (PRFS). CFH and CFS were restricted to be between 2 and 150 days and between 10 and 200 days, respectively. CI records had to be between 300 and 800 days. Cows were coded as pregnant to first service if they had no second service, were not diagnosed as ‘not pregnant’ by transrectal ultrasound and if the number of days from first service to next calving was between 265 and 300 days. Because of the experimental treatments imposed, no traditional fertility information other than CFH was available from the Netherlands.
Post-partum interval to the commencement of luteal activity (CLA) was defined as the number of days from calving to the first occurrence of two consecutive test-day records with a milk progesterone concentration of ⩾3 ng/ml (Royal et al., Reference Royal, Darwash, Flint, Webb, Woolliams and Lamming2000). Number of records per trait are summarised in Table 1.
Heritability estimated from the entire dataset as well as the number of records, mean and phenotypic standard deviation for each of the fertility traits within each country
CFH = calving to first heat; CFS = calving to first service; CI = calving interval; NS = number of services; PRFS = pregnancy rate to first service; CLA = commencement of luteal activity.
†No data available on calving to first service, number of services, pregnancy rate to first service or calving interval for the Netherlands.
‡No data available on calving to first heat for Sweden.
Genotypic data
All animals had phenotypic information on a range of performance traits and were genotyped with the Illumina BovineSNP50 BeadChip (Illumina Inc., San Diego, CA, USA) containing 54 001 single nucleotide polymorphisms (SNPs). A total of 65 cows with genotype call rates <95% were not included. SNPs were included in the association analyses if they fulfilled all of the following quality criteria: (1) GCscore >0.20 and GTscore >0.55; (2) call rate >0.95%; (3) minor allele frequency >0.01 in each country; and (4) no extreme deviation from Hardy Weinberg Equilibrium (i.e. χ 2 test statistic <600). The GCscore and GTscore are quality measures on the genotype calls from the genotyping assay. Following these edits, 37 590 SNPs from 1956 cows remained. Checks for Mendelian inconsistencies between pedigree and SNP data were first performed for all genotyped parent–offspring pairs and among sibs (Calus et al., Reference Calus, Mulder and Bastiaansen2011b). Following the removal of animals that did not pass parentage verification, a total of 1887 cows from Ireland, the UK, the Netherlands and Sweden remained. Chromosome number and positions of the SNPs on the BovineSNP50 were obtained from the UMD3.0 bovine genome assembly from the University of Maryland.
Variance component estimation
Genetic and residual (co)variances for the different fertility traits were estimated using animal linear mixed models in ASREML (Gilmour et al., Reference Gilmour, Cullis, Welham and Thompson2009). Fixed effects included in the model for all traits were the contemporary groups of country-experimental treatment-year and country-year-season of calving. For the estimation of variance components for PRFS, days since calving at the time of insemination was also included as a class effect (<30 days, 31 to 60 days, 61 to 90 days, 91 to 120 days, >120 days) in the model to account for differing pregnancy rate by days post calving (Berry et al., Reference Berry, Evans and Mc Parland2011). Genetic and residual covariances were estimated using 15 bivariate analyses. Pedigree information on each genotyped animal was traced back to at least four generations where available; a total of 7619 animals were included in the pedigree file.
Genome-wide association analysis
The genome-wide association analysis was conducted using a model that estimates effects for all SNPs simultaneously. The general model used to estimate SNP effects was
![\[--><$$>\eqalign {{{y}_{ij}}\, = &\,{{\mu }_j}\, + \,fixed\;effect{{s}_{ij}}\, + \,anima{{l}_{ij}}\, \cr &+ \,\mathop{\sum}\limits_{k\, = \,1}^{nloc} {\mathop{\sum}\limits_{l\, = \,1}^2 {SN{{P}_{ijkl}}\, + \,{{e}_{ij}}} } \eqno<$$><!--\]](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20160920214515743-0706:S1751731112000067:S1751731112000067_eqnU1.gif?pub-status=live)
where yij is the phenotypic record of animal i, μj is the overall mean for trait j, fixed effectsij are the fixed effects pertinent to animal i for trait j, animalij is the random polygenic effect of animal i for trait j, SNPijkl is a random effect for allele l on trait j at locus k of animal i and eij is a random residual for animal i. The fixed effects included in the models were the same as those used in the estimation of the variance components. This model was used both for univariate and bivariate analyses.
The polygenic effect was fitted with a numerator pedigree relationship matrix. SNP effects, denoted as SNPijkl , were estimated as qijkl × vjk , as described by Meuwissen and Goddard (Reference Meuwissen and Goddard2004), where qijkl is the size of the effect of allele l at locus k and vjk is a scaling factor in the direction vector for locus k that scales the effect at locus k for trait j. The variance of the direction vector v.k is denoted as V... The prior distribution for V in this case was
![\[--><$$>p{\rm{(}}{\bf {{V}}}..{\rm{)}}\,{\rm{ = }}\,{{\chi }^{{\rm{ - 2}}}} \,{\rm{(}}{{{\bf {{S}}}}_{\rm{0}}}{\rm{(}}..{\rm{), 10)}}\eqno<$$><!--\]](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20160920214515743-0706:S1751731112000067:S1751731112000067_eqnU2.gif?pub-status=live){kind=small}
where S0 ( no ) is chosen such that it reflects the total genetic (co)variance between traits n and o, divided by the total number of SNPs. V.. was sampled from the following conditional m variate-inverted Wishart distribution with (nloc + 10) degrees of freedom:
![\[--><$$> {\bf {{V}}}..{\rm{|}}{\bf {{v}}}..{\rm{,}}\,I{\rm{.}}\,{\rm{\sim }}\,I{{W}_m}\,{\rm{( }}{{{\bf {{S}}}}_{\rm{0}}}{\rm{(}}..{\rm{)}}\,{\rm{ + }}\,{\bf {{S}}}{{{\bf {{Z}}}}_{{\rm{(}}..{\rm{)}}}}\,{\rm{,}}\,{}nloc\,{\rm{ - }}\,m\,{\rm{ - }}\,{\rm{1}}\,{\rm{ + }}\,{\rm{10)}} \eqno<$$><!--\]](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20160920214515743-0706:S1751731112000067:S1751731112000067_eqnU3.gif?pub-status=live){kind=small}
where 
, nloc = number of evaluated marker loci and 10 is the number of degrees of freedom for the prior distribution.
In the model, a QTL indicator (Ik ) was sampled for each locus k, which either had the value of 0 or 1. S0 ( .. ) is chosen such that it reflected the total genetic (co)variances of traits n and o, divided by the total number of expected QTLs instead of the number of SNPs. V.. was sampled from an inverted Wishart distribution as follows:
![\[--><$$>{\bf {{S}}}{{{\bf {{Z}}}}_{(..)}}\, = \,\mathop{\sum}\limits_{k\, = \,1}^{nloc} {{{{v^{\prime}}}_{.k}}{{v}_{.k}}({{I}_k}\, + \,(1\,{\rm{ - }}\,{{I}_k})\,\times \,100)} \eqno<$$><!--\]](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20160920214515743-0706:S1751731112000067:S1751731112000067_eqnU5.gif?pub-status=live){kind=small}
where the QTL indicator Ik was sampled from
![\[--><$$> \left[ {\frac{{\phi \,\left( {{{v}_{.k}};\,{\bf {{0}}},\,{\bf {{V}}}} \right)\,\times \,{{p}_k}}}{{\phi \,\left( {{{v}_{.k}};\,{\bf {{0}}},\,{\bf {{V}}}} \right)\,\times \,{{p}_k}\, + \,\phi \,\left( {{{v}_{.k}};\,{\bf {{0}}},\,\frac{{\bf {{V}}}}{{ 100}}} \right)\,\times \,\left( {1\,{\rm{ - }}\,{{p}_k}} \right)}}} \right]\eqno<$$><!--\]](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20160920214515743-0706:S1751731112000067:S1751731112000067_eqnU7.gif?pub-status=live){kind=small}
where pk is the prior QTL probability, that is, the prior probability that Ik is equal to 1. Prior QTL probabilities used in the analyses reflected the prior assumption that 100 QTL underlie each of the traits. The model is referred to as Bayesian Stochastic Search Variable Selection (BSSVS; e.g. Verbyla et al., Reference Verbyla, Hayes, Bowman and Goddard2009).
The residual variance was assumed to be normally distributed N(0, R), where R is the m × m residual covariance matrix and m is the number of traits (1 or 2). The polygenic variance was assumed to be normally distributed N(0, A⊗GA ), where A is the numerator relationship matrix and GA is the m × m polygenic covariance matrix. Matrices R and GA were both sampled in the Gibbs sampler from an inverted Wishart distribution, with a uniform prior distribution. The BSSVS model was performed using Gibbs sampling. All univariate models were run for 50 000 cycles, discarding the initial 10 000 as burn-in. All bivariate models were run for 100 000 cycles, discarding the initial 20 000 as burn-in. To limit the required computational effort, only one Markov chain Monte Carlo chain was run for all models, except for those including CLA. For those, four chains were run for all univariate models and three chains were run for all bivariate models. For CLA, results were averaged across chains.
The BSSVS model yields posterior QTL probabilities for each SNP locus that are calculated as the average of the QTL indicator Ik across all cycles after the burn-in. Bayes factors (Kass and Raftery, Reference Kass and Raftery1995) were used to assess the posterior QTL probabilities for all SNP loci. The Bayes Factor, B12 , was calculated as follows. Let y be the data, H1 is the hypothesis that the marker is linked to a QTL and H2 is the hypothesis that the marker is not linked to a QTL. Then Pr(H1 ) is the prior probability of the first hypothesis H1 and Pr(H2 ) is the prior probability of the alternative hypothesis, H2 . Similarly, Pr(H1| y ) is the posterior probability of H1 and Pr(H2| y ) is the posterior probability of H2 . Using Bayes theorem, a Bayes factor comparing hypotheses H1 and H2 , for a single SNP, is
![\[--><$$>{{B}_{12}}\,{{ = \,\frac{{Pr\left( {{{H}_1}\left| y \right.} \right)}}{{Pr\left( {{{H}_2}\left| y \right.} \right)}}} \mathord{\left/ {\vphantom {{ = \,\frac{{Pr\left( {{{H}_1}\left| y \right.} \right)}}{{Pr\left( {{{H}_2}\left| y \right.} \right)}}} {\frac{{Pr\left( {{{H}_1}} \right)}}{{Pr\left( {{{H}_2}} \right)}}}}} \right. \kern-\nulldelimiterspace} {\frac{{Pr\left( {{{H}_1}} \right)}}{{Pr\left( {{{H}_2}} \right)}}}}\eqno<$$><!--\]](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20160920214515743-0706:S1751731112000067:S1751731112000067_eqnU8.gif?pub-status=live){kind=small}
As this is the ratio of the posterior odds to the prior odds, it can be expressed as
![\[--><$$>{{B}_{12}}\, = \,{{\frac{{Pr\left( {{{H}_1}\left| y \right.} \right)}}{{1\,{\rm{ - }}\,Pr\left( {{{H}_1}\left| y \right.} \right)}}} \mathord{\left/ {\vphantom {{\frac{{Pr\left( {{{H}_1}\left| y \right.} \right)}}{{1\,{\rm{ - }}\,Pr\left( {{{H}_1}\left| y \right.} \right)}}} {\frac{{Pr\left( {{{H}_1}} \right)}}{{1\,{\rm{ - }}\,Pr\left( {{{H}_1}} \right)}}}}} \right. \kern-\nulldelimiterspace} {\frac{{Pr\left( {{{H}_1}} \right)}}{{1\,{\rm{ - }}\,Pr\left( {{{H}_1}} \right)}}}}\eqno<$$><!--\]](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20160920214515743-0706:S1751731112000067:S1751731112000067_eqnU9.gif?pub-status=live){kind=small}
Bayes factors >3.1 indicate ‘substantial evidence’ in favour of the first hypothesis, that is, that the SNP is associated with a QTL. Furthermore, Bayes factors >10.1 were assumed to indicate ‘strong evidence’ for association with a QTL, and Bayes factors >30.1 were assumed to indicate ‘very strong evidence’ (Jeffreys, Reference Jeffreys1961).
Results
Summary statistics for each of the fertility traits within each country are in Table 1, together with heritability estimates across all data. Heritability estimates for the traditional fertility traits varied from 0.03 (PRFS) to 0.16 (interval from CFH). The heritability of CLA was 0.13.
Genetic correlations among traits were generally stronger than their respective phenotypic correlation and were, in general, of the same sign (Table 2). The interval traits (i.e. CFH, CFS, CI and CLA) were all strongly genetically correlated (0.37 to 0.99) with each other; the respective phenotypic correlations varied from 0.12 to 0.63. Genetically longer CIs were associated with poor PRFS and more services (Table 2).
Genetic (above diagonal; standard errors in parenthesis) and phenotypicFootnote † (below diagonal) correlations among the fertility traits
CFH = calving to first heat; CFS = calving to first service; CI = calving interval; NS = number of services; PRFS = pregnancy rate to first service; CLA = commencement of luteal activity.
† Standard errors of the phenotypic correlations were ⩽0.052.
Genome-wide associations
The posterior QTL probabilities at each SNP were <0.021, with the exception of CLA (Figures 1 and 2). The maximum posterior probability for CFH was 0.020 for SNP ARS-BFGL-NGS-104801 located on BTA 8; the next highest posterior probability for CFH was 0.014 for SNP Hapmap50203-BTA-107053 located on BTA 4. The maximum posterior probability for CFS was 0.010 for SNP BTB-01412031 located on BTA 20, whereas the next highest posterior probability for CFS was 0.009 for SNP ARS-BFGL-NGS-13749 located on BTA 6. The maximum posterior probability for NS was 0.006 for Hapmap39421-BTA-65907 located on BTA 4, whereas the maximum posterior probability for CI was 0.005 for ARS-BFGL-NGS-14867 located on BTA 19. The maximum posterior probability for PRFS was 0.005 for BTB-01390173 located on BTA 13. The proportion of the genetic variance explained by the top two SNPs in each of the traits varied from 0.002% to 0.012%.
Posterior quantitative trait locus probabilities, from top to bottom, for post-partum commencement of luteal activity, calving to first heat interval and calving to first service interval for each of the autosomes from left (BTA 1) to right (BTA 29) and the X-chromosome at the far right.
Posterior quantitative trait locus probabilities, from top to bottom, for calving interval, number of services and pregnancy rate to first service for each of the autosomes from left (BTA 1) to right (BTA 29) and the X-chromosome at the far right.
Posterior QTL probabilities of >0.04 were observed for CLA (Table 3) on BTA 2 (BTA-49769-no-rs; posterior probability of 0.060) and BTA 21 (BTA-12468-no-rs; posterior probability of 0.045). The SNP on BTA 2 explained 0.51% of the genetic variance in CLA, whereas the SNP on BTA 21 explained 0.35% of the genetic variance in CLA. The Bayes factors of BTA-49769-no-rs and BTA-12468-no-rs were 24 and 18, respectively, suggesting strong evidence of a QTL.
Ten SNPs with greatest posterior probabilities for associations with CLA from the univariate analysis, as well as their respective posterior probabilities when included in each of the five bivariate analyses. Also included is the chromosome (BTA) and position on the chromosomes for each SNP
SNPs = single nucleotide polymorphisms; CFH = calving to first heat; CFS = calving to first service; CI = calving interval; NS = number of services; PRFS = pregnancy rate to first service; CLA = commencement of luteal activity.
Figure 3 illustrates the QTL probabilities at each SNP for the five bivariate analyses that included CLA. The posterior probabilities of the 10 most significant SNPs from the univariate analysis of CLA when included in each of the bivariate analyses are detailed in Table 3. The large posterior QTL probabilities for CLA on BTA 2 and BTA 21 remained on the same SNPs, but they were larger compared with the univariate analysis of CLA. The posterior probability of 0.060 for CLA at SNP BTA-49769-no-rs was increased to 0.094, 0.121, 0.162, 0.662 and 0.162 when CLA was included in a bivariate analysis with CFH, CFS, NS, CI and PRFS, respectively (Figure 3 and Table 3). The posterior probability of 0.045 for CLA at SNP BTA-12468-no-rs on BTA 21 was increased to 0.052, 0.152, 0.072, 0.123 and 0.135 when CLA was included in a bivariate analysis with CFH, CFS, NS, CI and PRFS, respectively (Figure 3 and Table 3). The Bayes factor for BTA-49769-no-rs and BTA-12468-no-rs when included in the bivariate analysis with CI was 880 and 52, respectively (Figure 4), indicating a very strong support of a QTL.
Posterior quantitative trait locus probabilities from the bivariate analyses with commencement of luteal activity, from top to bottom, for calving to first heat interval, calving to first service interval, calving interval, number of services and pregnancy rate to first service for each of the autosomes from left (BTA 1) to right (BTA 29) and the X-chromosome at the far right.
Bayes factor from the bivariate analyses with commencement of luteal activity, from top to bottom, for calving to first heat interval, calving to first service interval, calving interval, number of services and pregnancy rate to first service for each of the autosomes from left (BTA 1) to right (BTA 29) and the X-chromosome at the far right. Horizontal line represents a Bayes factor of 10.
Table 4 summarises the correlations between SNP effects and genomic breeding values of animals from the pairwise bivariate analyses. The correlations between SNP effects were similar to the respective correlations between animal genomic breeding values. Correlations were also similar to the respective correlations estimated using traditional linear models using the numerator relationship matrix based on recorded pedigree (Table 2).
Correlations between estimated SNP effects (above diagonal) and animal genomic breeding values (below diagonal) from the pairwise bivariate analyses
SNPs = single nucleotide polymorphisms; CFH = calving to first heat; CFS = calving to first service; CI = calving interval; NS = number of services; PRFS = pregnancy rate to first service; CLA = commencement of luteal activity.
Discussion
The success of genome-wide association studies is a function of, amongst others, the heritability of the trait under investigation and the number of phenotypic records on that trait (Daetwyler et al., Reference Daetwyler, Villanueva and Woolliams2008), as well as the genetic architecture of the trait (i.e. influenced by a large or small number of genes). Large population sizes for genome-wide association analysis are required for lower heritability traits such as traditional measures of fertility (Pryce and Veerkamp, Reference Pryce and Veerkamp2001; Berry et al., Reference Berry, Buckley, Dillon, Evans, Rath and Veerkamp2003; Wall et al., Reference Wall, Brotherstone, Woolliams, Banos and Coffey2003); this can be somewhat overcome by using fertility-estimated breeding values on genotyped sires as the phenotype under investigation; nonetheless, phenotypes on a large number of cows are still required to generate accurate sire-estimated breeding values. The low heritability of fertility traits can be due to a multitude of factors including: (1) the trait is truly under poor additive genetic control relative to non-additive and non-genetic factors; (2) the statistical model used to estimate the variance components is too simplistic relative to reality; (3) recording errors adding to random noise; and (4) preferential treatment that cannot be fully accounted for in the statistical model. The use of hormonal profiles, related to fertility performance, can aid in circumventing some of the limitations of traditional measures, and previous research has shown that heritability estimates for fertility traits derived from progesterone profiles tend to be greater than traditional measures of fertility (Veerkamp et al., Reference Veerkamp, Oldenbroek, Van Der Gaast and Van Der Werf2000; Royal et al., Reference Royal, Flint and Woolliams2002; Petersson et al., Reference Petersson, Berglund, Strandberg, Gustafsson, Flint, Wooliams and Royal2007). Nonetheless, the quantity of data on progesterone profiles at an individual animal level is limited mainly because of the expense of collecting these data. The objective of this study, therefore, was to collate fertility information, including that derived from routine progesterone assays, on Holstein–Friesian cows from multiple countries and attempt to locate genomic regions associated with these traits.
Phenotypic performance
Mean phenotypic fertility performance in the present study was similar to previous reports from Ireland (Horan et al., Reference Horan, Mee, O'Connor, Rath and Dillon2005), the UK (Pollot and Coffey, Reference Pollot and Coffey2008), the Netherlands (Veerkamp et al., Reference Veerkamp, Oldenbroek, Van Der Gaast and Van Der Werf2000) and Sweden (Petersson et al., Reference Petersson, Strandberg, Gustafsson and Berglund2006), all of which included some of the data used in the present study. However, the mean CLA was longer in the present study compared with the aforementioned studies, but this is more than likely attributable to only first parity animals being included in the present study unlike the other studies, with the exception of Veerkamp et al. (Reference Veerkamp, Oldenbroek, Van Der Gaast and Van Der Werf2000), which included pluriparae as well. Petersson et al. (Reference Petersson, Strandberg, Gustafsson and Berglund2006) reported longer CLA in primiparae compared to pluriparae.
Genetic parameters
The low heritability of the traditional fertility traits in the present study corroborates the many previous studies on genetic parameters for traditional fertility measures (Pryce and Veerkamp, Reference Pryce and Veerkamp2001; Berry et al., Reference Berry, Buckley, Dillon, Evans, Rath and Veerkamp2003; Wall et al., Reference Wall, Brotherstone, Woolliams, Banos and Coffey2003). The heritability estimates of CFH and CFS, however, are slightly larger than that generally observed from field data, which may be attributable to an expected greater accuracy of recording in experimental herds. The heritability estimate for CLA in the present study (0.13) is slightly lower than estimates of 0.16 reported elsewhere (Veerkamp et al., Reference Veerkamp, Oldenbroek, Van Der Gaast and Van Der Werf2000; Royal et al., Reference Royal, Flint and Woolliams2002; Petersson et al., Reference Petersson, Berglund, Strandberg, Gustafsson, Flint, Wooliams and Royal2007), although some degree of similarity is expected as the data used in the present study was also included in some previous estimates (Veerkamp et al., Reference Veerkamp, Oldenbroek, Van Der Gaast and Van Der Werf2000).
The phenotypic correlations among the different fertility traits are in general agreement with previous studies. Using data independent of the present study, Royal et al. (Reference Royal, Flint and Woolliams2002) reported identical correlations to those in the present study of 0.12 and −0.03 between CLA and both CFS interval and PRFS, respectively. Genetic correlations among the traditional fertility variables were generally similar to those estimated on primiparae in larger datasets (Wall et al., Reference Wall, Brotherstone, Woolliams, Banos and Coffey2003). Nonetheless, the genetic correlations between the traditional fertility traits and CLA were different from those reported elsewhere, albeit the standard errors in the present study and in most of the other studies were large, attributable to a combination of low heritability of the fertility traits and the relatively small datasets. Royal et al. (Reference Royal, Flint and Woolliams2002) reported a correlation for CLA of −0.03 (s.e. = 0.27) and 0.49 (s.e. = 0.47) with CFS and PRFS, respectively, which are of opposite sign to the genetic correlations reported in the present study. Nonetheless, the genetic correlations in the present study were substantiated by the correlations between both the SNP effects and the genomic breeding values estimated in the bivariate genome-wide association analysis, where the covariances assumed between the traits in the analyses were estimated separately within the analyses themselves.
Considerable genetic variation was present for the different fertility traits evaluated in the present study, indicating that substantial genetic improvement can be achieved if the accuracy of identifying genetically superior animals is high. Achieving high accuracy of selection can be achieved through extensive phenotyping (of relatives and) of animals or through the use of genomics to augment the information from relatives and the animal itself.
Genome-wide associations
The model used to estimate the genome-wide associations in the present study is a model that is used to estimate genomic breeding values in genomic selection. Its performance for this purpose has been presented elsewhere for univariate applications (Calus et al., Reference Calus, Meuwissen, de Roos and Veerkamp2008), and for univariate and bivariate applications compared with other models, both for genomic prediction (Calus and Veerkamp, Reference Calus and Veerkamp2011) and QTL mapping (Calus et al., Reference Calus, Mulder and Veerkamp2011a). By estimating all SNP effects simultaneously, the rate of false positives is controlled and estimated SNP effects are less prone to overestimation compared with models where effects of only one SNP at the time are estimated. The dataset used in the present study could be considered relatively small for the detection of SNPs associated with low heritability traits. However, clear associations were identified with at least two SNPs. The association analysis for CLA was also undertaken with each country's data removed from the analysis, with the exception of the data from the Netherlands, which would have resulted in a large loss of data. The two SNPs BTA-12468-no-rs and BTA-49769-no-rs remained the most strongly associated SNPs with CLA, although the posterior probabilities were reduced to 0.007 (Swedish data removed) to 0.054 (Irish data removed) for BTA-49769-no-rs and to 0.018 (UK data removed) to 0.045 (Irish data removed) for BTA-12468-no-rs.
Although relatively large Bayes factors were obtained for CLA at two SNPs in the univariate analysis, these factors increased when CLA was included in a bivariate analysis. The degree of increase was greatest for CI and PRFS. Part of this increase may be attributable to additional information from the second trait, especially on genotyped animals that had no phenotype for CLA. Phenotypes for either of the traditional fertility traits were available on 1570 animals, whereas phenotypes for CLA were available on only 898 animals.
The BTA-49769-no-rs SNP that had a Bayes factor of 880 when estimated in the bivariate analysis with CI is a ‘decisive’ indication (Jeffreys, Reference Jeffreys1961) of a QTL at this location. The SNP is in an intergenic region on BTA 2 and a search for genes within 500 kb from this SNP revealed 16 genes at distances ranging from 9 to 462 kb; one gene was within 500 kb of BTA-12468-no-rs at a distance of 433 kb downstream. All genes were screened using Entrez Gene for functional annotation (http://www.ncbi.nlm.nih.gov/sites/entrez?db=gene) that would relate the genes to fertility traits in either cattle or humans.
The only gene near BTA-12468-no-rs was Bos taurus arrestin domain containing 4 (ARRDC4), which has a function in neurological processes. ARRDC4 was relatively distant from BTA-12468-no-rs compared with the average distance at which linkage disequilibrium (LD) decays in Holstein cattle (McKay et al., Reference McKay, Schnabel, Murdoch, Matukumalli, Aerts, Coppieters, Crews, Dias, Neto, Gill, Gao, Mannen, Stothard, Wang, Van Tassell, Williams, Taylor and Moore2007).
Of the 16 genes near BTA-49769-no-rs, three have functions that relate to fertility. The first gene, 242 kb upstream from BTA-49769-no-rs, was ‘Bos taurus inhibitor of DNA binding 3, dominant negative helix-loop-helix protein’ (ID3), which controls the growth of vascular cells through signalling that is induced by oestrogen (Felty and Porther, Reference Felty and Porther2008). However, studies that have attempted to relate mutations in the ID3 gene to fertility in cattle were not found. The second gene, located 462 kb upstream from BTA-49769-no-rs, was Bos taurus fucosidase, alpha-l-1, tissue (FUCA1), which is broadly distributed over the membrane systems of human sperm and has a role in the intimate species signature interactions between the sperm and oocyte (Venditti et al., Reference Venditti, Donigan and Bean2007). Even though FUCA1 was found to play a clear role in fertility, the evidence reported suggests a role in male rather than female fertility; although both are correlated (Berry et al., Reference Berry, Evans and Mc Parland2011). The third gene, located 199 kb upstream from BTA-49769-no-rs, was E2F transcription factor 2, which plays a crucial role in the control of cell cycle and action of tumour suppressor proteins in humans.
Several other studies have attempted to locate genomic regions associated with reproductive performance in cattle (Holmberg and Andersson-Eklund, Reference Holmberg and Andersson-Eklund2006; Huang et al., Reference Huang, Kirkpatrick, Rosa and Khatib2010; Pryce et al., Reference Pryce, Bolormaa, Chamberlain, Bowman, Savin, Goddard and Hayes2010; Sahana et al., Reference Sahana, Guldbrandtsen, Bendixen and Lund2010; Schulman et al., Reference Schulman, Sahana, Iso-Touro, McKay, Schnabel, Lund, Taylor, Virta and Vikki2011). Although not all studies have documented significant associations with different measures of reproductive performance on either BTA 2 or BTA 21 (Holmberg and Andersson-Eklund, Reference Holmberg and Andersson-Eklund2006; Pryce et al., Reference Pryce, Bolormaa, Chamberlain, Bowman, Savin, Goddard and Hayes2010), significant associations with SNPs on BTA 2 and BTA 21 have been reported in some studies (Huang et al., Reference Huang, Kirkpatrick, Rosa and Khatib2010; Sahana et al., Reference Sahana, Guldbrandtsen, Bendixen and Lund2010; Schulman et al., Reference Schulman, Sahana, Iso-Touro, McKay, Schnabel, Lund, Taylor, Virta and Vikki2011). Of particular note were the two SNPs, ARS-BFGL-NGS-36151 and BTB-00117780, reported by Huang et al. (Reference Huang, Kirkpatrick, Rosa and Khatib2010) on BTA 2 to be associated with fertilisation rate in a Holstein population independent to that used in the present study. On the basis of the UMD3.0 bovine genome assembly from the University of Maryland and the SNPs retained in the present study for inclusion in the association analysis, both SNPs reported by Huang et al. (Reference Huang, Kirkpatrick, Rosa and Khatib2010) flanked the BTA-49769-no-rs SNP shown to be associated with CLA in the present study. The ARS-BFGL-NGS-36151 SNP was the fifth SNP upstream (214 kb) from BTA-49769-no-rs, whereas BTB-00117780 was the 19th SNP downstream (1.25 Mb) from BTA-49769-no-rs. LD (r 2) between all three SNPs in the population of animals used in the present study was low (0.04 to 0.07). However, it is still possible that all three SNPs could be individually tagging the same causative mutation(s) in the region, and that the LD between associated SNP in both studies is simply low because of differences in local LD patterns between both populations.
Conclusion
The low heritability of the traditional fertility traits agree, albeit they were slightly greater in the present study, with the general consensus in the literature that traditional fertility measures are lowly heritable. Two regions of the genome with a relatively high probability of harbouring a mutation associated with CLA were identified, with the probability increasing when the genome-wide association was undertaken in a bivariate analysis; nonetheless, the proportion of the genetic variance explained by these regions was small.
Acknowledgements
This study is part of the RobustMilk project (http://www.robustmilk.eu), which is financially supported by the European Commission under the Seventh Research Framework Programme, Grant Agreement KBBE-211708. The content of this paper is the sole responsibility of the authors, and it does not necessarily represent the views of the Commission or its services.
